#include <cmath>
#include <limits>
-//#define __myDEBUG
+#define __myDEBUG
using namespace std;
enum UIndex { U_TGT = 1, U_SRC, LEN_TGT };
const double theMinSmoothCosin = 0.1;
+ const double theSmoothThickToElemSizeRatio = 0.3;
+
+ bool needSmoothing( double cosin, double tgtThick, double elemSize )
+ {
+ return cosin * tgtThick > theSmoothThickToElemSizeRatio * elemSize;
+ }
/*!
* \brief SMESH_ProxyMesh computed by _ViscousBuilder for a SOLID.
bool IsOnEdge() const { return _2neibors; }
gp_XYZ Copy( _LayerEdge& other, SMESH_MesherHelper& helper );
void SetCosin( double cosin );
+ int NbSteps() const { return _pos.size() - 1; } // nb inlation steps
};
struct _LayerEdgeCmp
{
bool CheckPrisms() const;
};
+ //--------------------------------------------------------------------------------
+ /*!
+ * \brief Layers parameters got by averaging several hypotheses
+ */
+ struct AverageHyp
+ {
+ AverageHyp( const StdMeshers_ViscousLayers* hyp = 0 )
+ :_nbLayers(0), _nbHyps(0), _thickness(0), _stretchFactor(0)
+ {
+ Add( hyp );
+ }
+ void Add( const StdMeshers_ViscousLayers* hyp )
+ {
+ if ( hyp )
+ {
+ _nbHyps++;
+ _nbLayers = hyp->GetNumberLayers();
+ //_thickness += hyp->GetTotalThickness();
+ _thickness = Max( _thickness, hyp->GetTotalThickness() );
+ _stretchFactor += hyp->GetStretchFactor();
+ }
+ }
+ double GetTotalThickness() const { return _thickness; /*_nbHyps ? _thickness / _nbHyps : 0;*/ }
+ double GetStretchFactor() const { return _nbHyps ? _stretchFactor / _nbHyps : 0; }
+ int GetNumberLayers() const { return _nbLayers; }
+ private:
+ int _nbLayers, _nbHyps;
+ double _thickness, _stretchFactor;
+ };
+
//--------------------------------------------------------------------------------
typedef map< const SMDS_MeshNode*, _LayerEdge*, TIDCompare > TNode2Edge;
*/
struct _SolidData
{
+ typedef const StdMeshers_ViscousLayers* THyp;
TopoDS_Shape _solid;
- const StdMeshers_ViscousLayers* _hyp;
- TopoDS_Shape _hypShape;
+ TGeomID _index; // SOLID id
_MeshOfSolid* _proxyMesh;
- set<TGeomID> _reversedFaceIds;
- set<TGeomID> _ignoreFaceIds; // WOL FACEs and FACEs of other SOLIDS
+ list< THyp > _hyps;
+ list< TopoDS_Shape > _hypShapes;
+ map< TGeomID, THyp > _face2hyp; // filled if _hyps.size() > 1
+ set< TGeomID > _reversedFaceIds;
+ set< TGeomID > _ignoreFaceIds; // WOL FACEs and FACEs of other SOLIDs
- double _stepSize, _stepSizeCoeff;
+ double _stepSize, _stepSizeCoeff, _geomSize;
const SMDS_MeshNode* _stepSizeNodes[2];
TNode2Edge _n2eMap; // nodes and _LayerEdge's based on them
// map to find _n2eMap of another _SolidData by a shrink shape shared by two _SolidData's
map< TGeomID, TNode2Edge* > _s2neMap;
// edges of _n2eMap. We keep same data in two containers because
- // iteration over the map is 5 time longer than over the vector
+ // iteration over the map is 5 times longer than over the vector
vector< _LayerEdge* > _edges;
// key: an id of shape (EDGE or VERTEX) shared by a FACE with
vector< int > _endEdgeOnShape;
int _nbShapesToSmooth;
- double _epsilon; // precision for SegTriaInter()
+ // data of averaged StdMeshers_ViscousLayers parameters for each shape with _LayerEdge's
+ vector< AverageHyp > _hypOnShape;
+ double _maxThickness; // of all _hyps
+ double _minThickness; // of all _hyps
- TGeomID _index; // SOLID id, for debug
+ double _epsilon; // precision for SegTriaInter()
- _SolidData(const TopoDS_Shape& s=TopoDS_Shape(),
- const StdMeshers_ViscousLayers* h=0,
- const TopoDS_Shape& hs=TopoDS_Shape(),
- _MeshOfSolid* m=0)
- :_solid(s), _hyp(h), _hypShape(hs), _proxyMesh(m) {}
+ _SolidData(const TopoDS_Shape& s=TopoDS_Shape(),
+ _MeshOfSolid* m=0)
+ :_solid(s), _proxyMesh(m) {}
~_SolidData();
Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge& E,
iEnd = _endEdgeOnShape[ end ];
}
- bool GetShapeEdges(const TGeomID shapeID, size_t& edgeEnd, int* iBeg=0, int* iEnd=0 ) const;
+ bool GetShapeEdges(const TGeomID shapeID, size_t& iEdgeEnd, int* iBeg=0, int* iEnd=0 ) const;
void AddShapesToSmooth( const set< TGeomID >& shapeIDs );
};
// does it's job
SMESH_ComputeErrorPtr Compute(SMESH_Mesh& mesh,
const TopoDS_Shape& shape);
+ // check validity of hypotheses
+ SMESH_ComputeErrorPtr CheckHypotheses( SMESH_Mesh& mesh,
+ const TopoDS_Shape& shape );
// restore event listeners used to clear an inferior dim sub-mesh modified by viscous layers
void RestoreListeners();
private:
bool findSolidsWithLayers();
- bool findFacesWithLayers();
+ bool findFacesWithLayers(const bool onlyWith=false);
+ void getIgnoreFaces(const TopoDS_Shape& solid,
+ const StdMeshers_ViscousLayers* hyp,
+ const TopoDS_Shape& hypShape,
+ set<TGeomID>& ignoreFaces);
bool makeLayer(_SolidData& data);
bool setEdgeData(_LayerEdge& edge, const set<TGeomID>& subIds,
SMESH_MesherHelper& helper, _SolidData& data);
const bool toSort = false);
void findSimplexTestEdges( _SolidData& data,
vector< vector<_LayerEdge*> >& edgesByGeom);
+ void computeGeomSize( _SolidData& data );
bool sortEdges( _SolidData& data,
vector< vector<_LayerEdge*> >& edgesByGeom);
void limitStepSizeByCurvature( _SolidData& data );
return _surface->Value( uv.X(), uv.Y() ).XYZ();
}
};
+
} // namespace VISCOUS_3D
{
// TODO
return false;
+} // --------------------------------------------------------------------------------
+SMESH_ComputeErrorPtr
+StdMeshers_ViscousLayers::CheckHypothesis(SMESH_Mesh& theMesh,
+ const TopoDS_Shape& theShape,
+ SMESH_Hypothesis::Hypothesis_Status& theStatus)
+{
+ VISCOUS_3D::_ViscousBuilder bulder;
+ SMESH_ComputeErrorPtr err = bulder.CheckHypotheses( theMesh, theShape );
+ if ( err && !err->IsOK() )
+ theStatus = SMESH_Hypothesis::HYP_INCOMPAT_HYPS;
+ else
+ theStatus = SMESH_Hypothesis::HYP_OK;
+
+ return err;
+}
+// --------------------------------------------------------------------------------
+bool StdMeshers_ViscousLayers::IsShapeWithLayers(int shapeIndex) const
+{
+ bool isIn =
+ ( std::find( _shapeIds.begin(), _shapeIds.end(), shapeIndex ) != _shapeIds.end() );
+ return IsToIgnoreShapes() ? !isIn : isIn;
}
// END StdMeshers_ViscousLayers hypothesis
//================================================================================
}
return false;
}
+ //================================================================================
+ /*!
+ * \brief Computes mimimal distance of face in-FACE nodes from an EDGE
+ * \param [in] face - the mesh face to treat
+ * \param [in] nodeOnEdge - a node on the EDGE
+ * \param [out] faceSize - the computed distance
+ * \return bool - true if faceSize computed
+ */
+ //================================================================================
+
+ bool getDistFromEdge( const SMDS_MeshElement* face,
+ const SMDS_MeshNode* nodeOnEdge,
+ double & faceSize )
+ {
+ faceSize = Precision::Infinite();
+ bool done = false;
+
+ int nbN = face->NbCornerNodes();
+ int iOnE = face->GetNodeIndex( nodeOnEdge );
+ int iNext[2] = { SMESH_MesherHelper::WrapIndex( iOnE+1, nbN ),
+ SMESH_MesherHelper::WrapIndex( iOnE-1, nbN ) };
+ const SMDS_MeshNode* nNext[2] = { face->GetNode( iNext[0] ),
+ face->GetNode( iNext[1] ) };
+ gp_XYZ segVec, segEnd = SMESH_TNodeXYZ( nodeOnEdge ); // segment on EDGE
+ double segLen = -1.;
+ // look for two neighbor not in-FACE nodes of face
+ for ( int i = 0; i < 2; ++i )
+ {
+ if ( nNext[i]->GetPosition()->GetDim() != 2 &&
+ nNext[i]->GetID() < nodeOnEdge->GetID() )
+ {
+ // look for an in-FACE node
+ for ( int iN = 0; iN < nbN; ++iN )
+ {
+ if ( iN == iOnE || iN == iNext[i] )
+ continue;
+ SMESH_TNodeXYZ pInFace = face->GetNode( iN );
+ gp_XYZ v = pInFace - segEnd;
+ if ( segLen < 0 )
+ {
+ segVec = SMESH_TNodeXYZ( nNext[i] ) - segEnd;
+ segLen = segVec.Modulus();
+ }
+ double distToSeg = v.Crossed( segVec ).Modulus() / segLen;
+ faceSize = Min( faceSize, distToSeg );
+ done = true;
+ }
+ segLen = -1;
+ }
+ }
+ return done;
+ }
+
//--------------------------------------------------------------------------------
// DEBUG. Dump intermediate node positions into a python script
// HOWTO use: run python commands written in a console to see
return _error;
}
+//================================================================================
+/*!
+ * \brief Check validity of hypotheses
+ */
+//================================================================================
+
+SMESH_ComputeErrorPtr _ViscousBuilder::CheckHypotheses( SMESH_Mesh& mesh,
+ const TopoDS_Shape& shape )
+{
+ _mesh = & mesh;
+
+ if ( _ViscousListener::GetSolidMesh( _mesh, shape, /*toCreate=*/false))
+ return SMESH_ComputeErrorPtr(); // everything already computed
+
+
+ findSolidsWithLayers();
+ bool ok = findFacesWithLayers();
+
+ // remove _MeshOfSolid's of _SolidData's
+ for ( size_t i = 0; i < _sdVec.size(); ++i )
+ _ViscousListener::RemoveSolidMesh( _mesh, _sdVec[i]._solid );
+
+ if ( !ok )
+ return _error;
+
+ return SMESH_ComputeErrorPtr();
+}
+
//================================================================================
/*!
* \brief Finds SOLIDs to compute using viscous layers. Fills _sdVec
// TODO: check if algo is hidden
const list <const SMESHDS_Hypothesis *> & allHyps =
algo->GetUsedHypothesis(*_mesh, allSolids(i), /*ignoreAuxiliary=*/false);
+ _SolidData* soData = 0;
list< const SMESHDS_Hypothesis *>::const_iterator hyp = allHyps.begin();
const StdMeshers_ViscousLayers* viscHyp = 0;
- for ( ; hyp != allHyps.end() && !viscHyp; ++hyp )
- viscHyp = dynamic_cast<const StdMeshers_ViscousLayers*>( *hyp );
- if ( viscHyp )
- {
- TopoDS_Shape hypShape;
- filter.Init( filter.Is( viscHyp ));
- _mesh->GetHypothesis( allSolids(i), filter, true, &hypShape );
+ for ( ; hyp != allHyps.end(); ++hyp )
+ if ( viscHyp = dynamic_cast<const StdMeshers_ViscousLayers*>( *hyp ))
+ {
+ TopoDS_Shape hypShape;
+ filter.Init( filter.Is( viscHyp ));
+ _mesh->GetHypothesis( allSolids(i), filter, true, &hypShape );
- _MeshOfSolid* proxyMesh = _ViscousListener::GetSolidMesh( _mesh,
- allSolids(i),
- /*toCreate=*/true);
- _sdVec.push_back( _SolidData( allSolids(i), viscHyp, hypShape, proxyMesh ));
- _sdVec.back()._index = getMeshDS()->ShapeToIndex( allSolids(i));
- }
+ if ( !soData )
+ {
+ _MeshOfSolid* proxyMesh = _ViscousListener::GetSolidMesh( _mesh,
+ allSolids(i),
+ /*toCreate=*/true);
+ _sdVec.push_back( _SolidData( allSolids(i), proxyMesh ));
+ soData = & _sdVec.back();
+ soData->_index = getMeshDS()->ShapeToIndex( allSolids(i));
+ }
+ soData->_hyps.push_back( viscHyp );
+ soData->_hypShapes.push_back( hypShape );
+ }
}
if ( _sdVec.empty() )
return error
*/
//================================================================================
-bool _ViscousBuilder::findFacesWithLayers()
+bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith)
{
SMESH_MesherHelper helper( *_mesh );
TopExp_Explorer exp;
{
solids.Add( _sdVec[i]._solid );
- vector<TGeomID> ids = _sdVec[i]._hyp->GetBndShapes();
- if ( _sdVec[i]._hyp->IsToIgnoreShapes() ) // FACEs to ignore are given
+ // get faces to ignore defined by each hyp
+ typedef const StdMeshers_ViscousLayers* THyp;
+ typedef std::pair< set<TGeomID>, THyp > TFacesOfHyp;
+ list< TFacesOfHyp > ignoreFacesOfHyps;
+ list< THyp >::iterator hyp = _sdVec[i]._hyps.begin();
+ list< TopoDS_Shape >::iterator hypShape = _sdVec[i]._hypShapes.begin();
+ for ( ; hyp != _sdVec[i]._hyps.end(); ++hyp, ++hypShape )
{
- for ( size_t ii = 0; ii < ids.size(); ++ii )
- {
- const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[ii] );
- if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
- _sdVec[i]._ignoreFaceIds.insert( ids[ii] );
- }
+ ignoreFacesOfHyps.push_back( TFacesOfHyp( set<TGeomID>(), *hyp ));
+ getIgnoreFaces( _sdVec[i]._solid, *hyp, *hypShape, ignoreFacesOfHyps.back().first );
}
- else // FACEs with layers are given
+
+ // fill _SolidData::_face2hyp and check compatibility of hypotheses
+ const int nbHyps = _sdVec[i]._hyps.size();
+ if ( nbHyps > 1 )
{
- exp.Init( _sdVec[i]._solid, TopAbs_FACE );
- for ( ; exp.More(); exp.Next() )
+ // check if two hypotheses define different parameters for the same FACE
+ list< TFacesOfHyp >::iterator igFacesOfHyp;
+ for ( exp.Init( _sdVec[i]._solid, TopAbs_FACE ); exp.More(); exp.Next() )
{
- TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
- if ( find( ids.begin(), ids.end(), faceInd ) == ids.end() )
- _sdVec[i]._ignoreFaceIds.insert( faceInd );
+ const TGeomID faceID = getMeshDS()->ShapeToIndex( exp.Current() );
+ THyp hyp = 0;
+ igFacesOfHyp = ignoreFacesOfHyps.begin();
+ for ( ; igFacesOfHyp != ignoreFacesOfHyps.end(); ++igFacesOfHyp )
+ if ( ! igFacesOfHyp->first.count( faceID ))
+ {
+ if ( hyp )
+ return error(SMESH_Comment("Several hypotheses define "
+ "Viscous Layers on the face #") << faceID );
+ hyp = igFacesOfHyp->second;
+ }
+ if ( hyp )
+ _sdVec[i]._face2hyp.insert( make_pair( faceID, hyp ));
+ else
+ _sdVec[i]._ignoreFaceIds.insert( faceID );
}
- }
- // ignore internal FACEs if inlets and outlets are specified
+ // check if two hypotheses define different number of viscous layers for
+ // adjacent faces of a solid
+ set< int > nbLayersSet;
+ igFacesOfHyp = ignoreFacesOfHyps.begin();
+ for ( ; igFacesOfHyp != ignoreFacesOfHyps.end(); ++igFacesOfHyp )
+ {
+ nbLayersSet.insert( igFacesOfHyp->second->GetNumberLayers() );
+ }
+ if ( nbLayersSet.size() > 1 )
+ {
+ for ( exp.Init( _sdVec[i]._solid, TopAbs_EDGE ); exp.More(); exp.Next() )
+ {
+ PShapeIteratorPtr fIt = helper.GetAncestors( exp.Current(), *_mesh, TopAbs_FACE );
+ THyp hyp1 = 0, hyp2 = 0;
+ while( const TopoDS_Shape* face = fIt->next() )
+ {
+ const TGeomID faceID = getMeshDS()->ShapeToIndex( *face );
+ map< TGeomID, THyp >::iterator f2h = _sdVec[i]._face2hyp.find( faceID );
+ if ( f2h != _sdVec[i]._face2hyp.end() )
+ {
+ ( hyp1 ? hyp2 : hyp1 ) = f2h->second;
+ }
+ }
+ if ( hyp1 && hyp2 &&
+ hyp1->GetNumberLayers() != hyp2->GetNumberLayers() )
+ {
+ return error("Two hypotheses define different number of "
+ "viscous layers on adjacent faces");
+ }
+ }
+ }
+ } // if ( nbHyps > 1 )
+ else
{
- TopTools_IndexedDataMapOfShapeListOfShape solidsOfFace;
- if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
- TopExp::MapShapesAndAncestors( _sdVec[i]._hypShape,
- TopAbs_FACE, TopAbs_SOLID, solidsOfFace);
+ _sdVec[i]._ignoreFaceIds.swap( ignoreFacesOfHyps.back().first );
+ }
+ // fill _SolidData::_reversedFaceIds
+ {
exp.Init( _sdVec[i]._solid.Oriented( TopAbs_FORWARD ), TopAbs_FACE );
for ( ; exp.More(); exp.Next() )
{
const TopoDS_Face& face = TopoDS::Face( exp.Current() );
- if ( helper.NbAncestors( face, *_mesh, TopAbs_SOLID ) < 2 )
- continue;
-
- const TGeomID faceInd = getMeshDS()->ShapeToIndex( face );
- if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
+ const TGeomID faceID = getMeshDS()->ShapeToIndex( face );
+ if ( //!sdVec[i]._ignoreFaceIds.count( faceID ) && ???????
+ helper.NbAncestors( face, *_mesh, TopAbs_SOLID ) > 1 &&
+ helper.IsReversedSubMesh( face ))
{
- int nbSolids = solidsOfFace.FindFromKey( face ).Extent();
- if ( nbSolids > 1 )
- _sdVec[i]._ignoreFaceIds.insert( faceInd );
- }
-
- if ( helper.IsReversedSubMesh( face ))
- {
- _sdVec[i]._reversedFaceIds.insert( faceInd );
+ _sdVec[i]._reversedFaceIds.insert( faceID );
}
}
}
- }
+ } // loop on _sdVec
+
+ if ( onlyWith ) // is called to check hypotheses compatibility only
+ return true;
// Find faces to shrink mesh on (solution 2 in issue 0020832);
TopTools_IndexedMapOfShape shapes;
return true;
}
+//================================================================================
+/*!
+ * \brief Finds FACEs w/o layers for a given SOLID by an hypothesis
+ */
+//================================================================================
+
+void _ViscousBuilder::getIgnoreFaces(const TopoDS_Shape& solid,
+ const StdMeshers_ViscousLayers* hyp,
+ const TopoDS_Shape& hypShape,
+ set<TGeomID>& ignoreFaceIds)
+{
+ TopExp_Explorer exp;
+
+ vector<TGeomID> ids = hyp->GetBndShapes();
+ if ( hyp->IsToIgnoreShapes() ) // FACEs to ignore are given
+ {
+ for ( size_t ii = 0; ii < ids.size(); ++ii )
+ {
+ const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[ii] );
+ if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
+ ignoreFaceIds.insert( ids[ii] );
+ }
+ }
+ else // FACEs with layers are given
+ {
+ exp.Init( solid, TopAbs_FACE );
+ for ( ; exp.More(); exp.Next() )
+ {
+ TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
+ if ( find( ids.begin(), ids.end(), faceInd ) == ids.end() )
+ ignoreFaceIds.insert( faceInd );
+ }
+ }
+
+ // ignore internal FACEs if inlets and outlets are specified
+ if ( hyp->IsToIgnoreShapes() )
+ {
+ TopTools_IndexedDataMapOfShapeListOfShape solidsOfFace;
+ TopExp::MapShapesAndAncestors( hypShape,
+ TopAbs_FACE, TopAbs_SOLID, solidsOfFace);
+
+ for ( exp.Init( solid, TopAbs_FACE ); exp.More(); exp.Next() )
+ {
+ const TopoDS_Face& face = TopoDS::Face( exp.Current() );
+ if ( SMESH_MesherHelper::NbAncestors( face, *_mesh, TopAbs_SOLID ) < 2 )
+ continue;
+
+ int nbSolids = solidsOfFace.FindFromKey( face ).Extent();
+ if ( nbSolids > 1 )
+ ignoreFaceIds.insert( getMeshDS()->ShapeToIndex( face ));
+ }
+ }
+}
+
//================================================================================
/*!
* \brief Create the inner surface of the viscous layer and prepare data for infation
void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
{
const int nbTestPnt = 5; // on a FACE sub-shape
- const double minCurvature = 0.9 / data._hyp->GetTotalThickness();
BRepLProp_SLProps surfProp( 2, 1e-6 );
SMESH_MesherHelper helper( *_mesh );
else
continue;
// check concavity and curvature and limit data._stepSize
+ const double minCurvature = 0.9 / data._hypOnShape[ edgesEnd ].GetTotalThickness();
int nbLEdges = iEnd - iBeg;
- int iStep = Max( 1, nbLEdges / nbTestPnt );
+ int iStep = Max( 1, nbLEdges / nbTestPnt );
for ( ; iBeg < iEnd; iBeg += iStep )
{
gp_XY uv = helper.GetNodeUV( F, data._edges[ iBeg ]->_nodes[0] );
bool _ViscousBuilder::sortEdges( _SolidData& data,
vector< vector<_LayerEdge*> >& edgesByGeom)
{
+ // define allowed thickness
+ computeGeomSize( data ); // compute data._geomSize
+
+ data._maxThickness = 0;
+ data._minThickness = 1e100;
+ list< const StdMeshers_ViscousLayers* >::iterator hyp = data._hyps.begin();
+ for ( ; hyp != data._hyps.end(); ++hyp )
+ {
+ data._maxThickness = Max( data._maxThickness, (*hyp)->GetTotalThickness() );
+ data._minThickness = Min( data._minThickness, (*hyp)->GetTotalThickness() );
+ }
+ const double tgtThick = /*Min( 0.5 * data._geomSize, */data._maxThickness;
+
// Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's
- // boundry inclined at a sharp angle to the shape
+ // boundry inclined to the shape at a sharp angle
list< TGeomID > shapesToSmooth;
TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() );
vector<_LayerEdge*>& eV = edgesByGeom[ iV ];
if ( eV.empty() ) continue;
- // double cosin = eV[0]->_cosin;
- // bool badCosin =
- // ( !eV[0]->_sWOL.IsNull() && ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE || !isShrinkEdge));
- // if ( badCosin )
- // {
- // gp_Vec dir1, dir2;
- // if ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE )
- // dir1 = getEdgeDir( TopoDS::Edge( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ));
- // else
- // dir1 = getFaceDir( TopoDS::Face( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ),
- // eV[0]->_nodes[0], helper, ok);
- // dir2 = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
- // double angle = dir1.Angle( dir2 );
- // cosin = cos( angle );
- // }
gp_Vec eDir = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
double angle = eDir.Angle( eV[0]->_normal );
double cosin = Cos( angle );
- needSmooth = ( cosin > theMinSmoothCosin );
+ if ( cosin > theMinSmoothCosin )
+ {
+ // compare tgtThick with the length of an end segment
+ SMDS_ElemIteratorPtr eIt = eV[0]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Edge);
+ while ( eIt->more() )
+ {
+ const SMDS_MeshElement* endSeg = eIt->next();
+ if ( endSeg->getshapeId() == iS )
+ {
+ double segLen =
+ SMESH_TNodeXYZ( endSeg->GetNode(0) ).Distance( endSeg->GetNode(1 ));
+ needSmooth = needSmoothing( cosin, tgtThick, segLen );
+ break;
+ }
+ }
+ }
}
break;
}
TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() );
vector<_LayerEdge*>& eE = edgesByGeom[ iE ];
if ( eE.empty() ) continue;
- if ( eE[0]->_sWOL.IsNull() )
+ // TopLoc_Location loc;
+ // Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face( S ), loc );
+ // bool isPlane = GeomLib_IsPlanarSurface( surface ).IsPlanar();
+ //if ( eE[0]->_sWOL.IsNull() )
{
+ double faceSize;
for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
- needSmooth = ( eE[i]->_cosin > theMinSmoothCosin );
- }
- else
- {
- const TopoDS_Face& F1 = TopoDS::Face( S );
- const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL );
- const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
- for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
- {
- gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok );
- gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok );
- double angle = dir1.Angle( dir2 );
- double cosin = cos( angle );
- needSmooth = ( cosin > theMinSmoothCosin );
- }
+ if ( eE[i]->_cosin > theMinSmoothCosin )
+ {
+ SMDS_ElemIteratorPtr fIt = eE[i]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() && !needSmooth )
+ {
+ const SMDS_MeshElement* face = fIt->next();
+ if ( getDistFromEdge( face, eE[i]->_nodes[0], faceSize ))
+ needSmooth = needSmoothing( eE[i]->_cosin, tgtThick, faceSize );
+ }
+ }
}
+ // else
+ // {
+ // const TopoDS_Face& F1 = TopoDS::Face( S );
+ // const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL );
+ // const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
+ // for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
+ // {
+ // gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok );
+ // gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok );
+ // double angle = dir1.Angle( );
+ // double cosin = cos( angle );
+ // needSmooth = ( cosin > theMinSmoothCosin );
+ // }
+ // }
}
break;
}
continue;
default:;
}
+
if ( needSmooth )
{
if ( S.ShapeType() == TopAbs_EDGE ) shapesToSmooth.push_front( iS );
//eVec.clear();
}
+ // compute average StdMeshers_ViscousLayers parameters for each shape
+
+ data._hypOnShape.clear();
+ if ( data._hyps.size() == 1 )
+ {
+ data._hypOnShape.resize( data._endEdgeOnShape.size(), AverageHyp( data._hyps.back() ));
+ }
+ else
+ {
+ data._hypOnShape.resize( data._endEdgeOnShape.size() );
+ map< TGeomID, const StdMeshers_ViscousLayers* >::iterator f2hyp;
+ for ( size_t i = 0; i < data._endEdgeOnShape.size(); ++i )
+ {
+ int iEnd = data._endEdgeOnShape[i];
+ _LayerEdge* LE = data._edges[ iEnd-1 ];
+ TGeomID iShape = LE->_nodes[0]->getshapeId();
+ const TopoDS_Shape& S = getMeshDS()->IndexToShape( iShape );
+ if ( S.ShapeType() == TopAbs_FACE )
+ {
+ if (( f2hyp = data._face2hyp.find( iShape )) != data._face2hyp.end() )
+ {
+ data._hypOnShape[ i ].Add( f2hyp->second );
+ }
+ }
+ else
+ {
+ PShapeIteratorPtr fIt = SMESH_MesherHelper::GetAncestors( S, *_mesh, TopAbs_FACE );
+ while ( const TopoDS_Shape* face = fIt->next() )
+ {
+ TGeomID faceID = getMeshDS()->ShapeToIndex( *face );
+ if (( f2hyp = data._face2hyp.find( faceID )) != data._face2hyp.end() )
+ {
+ data._hypOnShape[ i ].Add( f2hyp->second );
+ }
+ }
+ }
+ }
+ }
+
return ok;
}
#endif
}
+//================================================================================
+/*!
+ * \brief Find maximal _LayerEdge length (layer thickness) limited by geometry
+ */
+//================================================================================
+
+void _ViscousBuilder::computeGeomSize( _SolidData& data )
+{
+ data._geomSize = Precision::Infinite();
+ double intersecDist;
+ auto_ptr<SMESH_ElementSearcher> searcher
+ ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
+ data._proxyMesh->GetFaces( data._solid )) );
+
+ TNode2Edge::iterator n2e = data._n2eMap.begin(), n2eEnd = data._n2eMap.end();
+ for ( ; n2e != n2eEnd; ++n2e )
+ {
+ _LayerEdge* edge = n2e->second;
+ if ( edge->IsOnEdge() ) continue;
+ edge->FindIntersection( *searcher, intersecDist, data._epsilon );
+ if ( data._geomSize > intersecDist && intersecDist > 0 )
+ data._geomSize = intersecDist;
+ }
+}
+
//================================================================================
/*!
* \brief Increase length of _LayerEdge's to reach the required thickness of layers
// Limit inflation step size by geometry size found by itersecting
// normals of _LayerEdge's with mesh faces
- double geomSize = Precision::Infinite(), intersecDist;
- auto_ptr<SMESH_ElementSearcher> searcher
- ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
- data._proxyMesh->GetFaces( data._solid )) );
- for ( size_t i = 0; i < data._edges.size(); ++i )
- {
- if ( data._edges[i]->IsOnEdge() ) continue;
- data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon );
- if ( geomSize > intersecDist && intersecDist > 0 )
- geomSize = intersecDist;
- }
- if ( data._stepSize > 0.3 * geomSize )
- limitStepSize( data, 0.3 * geomSize );
+ if ( data._stepSize > 0.3 * data._geomSize )
+ limitStepSize( data, 0.3 * data._geomSize );
- const double tgtThick = data._hyp->GetTotalThickness();
- if ( data._stepSize > tgtThick )
- limitStepSize( data, tgtThick );
+ const double tgtThick = data._maxThickness;
+ if ( data._stepSize > data._minThickness )
+ limitStepSize( data, data._minThickness );
if ( data._stepSize < 1. )
data._epsilon = data._stepSize * 1e-7;
- debugMsg( "-- geomSize = " << geomSize << ", stepSize = " << data._stepSize );
+ debugMsg( "-- geomSize = " << data._geomSize << ", stepSize = " << data._stepSize );
double avgThick = 0, curThick = 0, distToIntersection = Precision::Infinite();
int nbSteps = 0, nbRepeats = 0;
- while ( 1.01 * avgThick < tgtThick )
+ int iBeg, iEnd, iS;
+ while ( avgThick < 0.99 )
{
// new target length
curThick += data._stepSize;
if ( curThick > tgtThick )
{
- curThick = tgtThick + ( tgtThick-avgThick ) * nbRepeats;
+ curThick = tgtThick + tgtThick*( 1.-avgThick ) * nbRepeats;
nbRepeats++;
}
// Elongate _LayerEdge's
dumpFunction(SMESH_Comment("inflate")<<data._index<<"_step"<<nbSteps); // debug
- for ( size_t i = 0; i < data._edges.size(); ++i )
+ for ( iBeg = 0, iS = 0; iS < data._endEdgeOnShape.size(); ++iS )
{
- data._edges[i]->SetNewLength( curThick, helper );
+ const double shapeCurThick = Min( curThick, data._hypOnShape[ iS ].GetTotalThickness() );
+ for ( iEnd = data._endEdgeOnShape[ iS ]; iBeg < iEnd; ++iBeg )
+ {
+ data._edges[iBeg]->SetNewLength( shapeCurThick, helper );
+ }
}
dumpFunctionEnd();
// Evaluate achieved thickness
avgThick = 0;
- for ( size_t i = 0; i < data._edges.size(); ++i )
- avgThick += data._edges[i]->_len;
+ for ( iBeg = 0, iS = 0; iS < data._endEdgeOnShape.size(); ++iS )
+ {
+ const double shapeTgtThick = data._hypOnShape[ iS ].GetTotalThickness();
+ for ( iEnd = data._endEdgeOnShape[ iS ]; iBeg < iEnd; ++iBeg )
+ {
+ avgThick += Min( 1., data._edges[iBeg]->_len / shapeTgtThick );
+ }
+ }
avgThick /= data._edges.size();
- debugMsg( "-- Thickness " << avgThick << " reached" );
+ debugMsg( "-- Thickness " << curThick << " ("<< avgThick*100 << "%) reached" );
- if ( distToIntersection < avgThick*1.5 )
+ if ( distToIntersection < tgtThick*avgThick*1.5 )
{
debugMsg( "-- Stop inflation since "
<< " distToIntersection( "<<distToIntersection<<" ) < avgThick( "
- << avgThick << " ) * 1.5" );
+ << tgtThick*avgThick << " ) * 1.5" );
break;
}
// new step size
data._stepSize = data._stepSizeCoeff *
SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
- } // while ( 1.01 * avgThick < tgtThick )
+ } // while ( avgThick < 0.99 )
if (nbSteps == 0 )
return error("failed at the very first inflation step", data._index);
- if ( 1.01 * avgThick < tgtThick )
+ if ( avgThick < 0.99 )
if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( data._index ))
{
SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
( new SMESH_ComputeError (COMPERR_WARNING,
SMESH_Comment("Thickness ") << tgtThick <<
" of viscous layers not reached,"
- " average reached thickness is " << avgThick ));
+ " average reached thickness is " << avgThick*tgtThick));
}
// Restore position of src nodes moved by infaltion on _noShrinkShapes
dumpFunction(SMESH_Comment("restoNoShrink_So")<<data._index); // debug
- int iBeg, iEnd = 0;
- for ( int iS = 0; iS < data._endEdgeOnShape.size(); ++iS )
+ for ( iEnd = iS = 0; iS < data._endEdgeOnShape.size(); ++iS )
{
iBeg = iEnd;
iEnd = data._endEdgeOnShape[ iS ];
iBeg = iEnd;
iEnd = data._endEdgeOnShape[ iS ];
+ // bool toSmooth = false;
+ // for ( int i = iBeg; i < iEnd; ++i )
+ // toSmooth = data._edges[ iBeg ]->NbSteps() >= nbSteps+1;
+ // if ( !toSmooth )
+ // {
+ // if ( iS+1 == data._nbShapesToSmooth )
+ // data._nbShapesToSmooth--;
+ // continue; // target length reached some steps before
+ // }
+
if ( !data._edges[ iBeg ]->_sWOL.IsNull() &&
data._edges[ iBeg ]->_sWOL.ShapeType() == TopAbs_FACE )
{
badNb = 0;
moved = false;
if ( step % 2 )
- for ( int i = iBeg; i < iEnd; ++i )
+ for ( int i = iBeg; i < iEnd; ++i ) // iterate forward
moved |= data._edges[i]->Smooth(badNb);
else
- for ( int i = iEnd-1; i >= iBeg; --i )
+ for ( int i = iEnd-1; i >= iBeg; --i ) // iterate backward
moved |= data._edges[i]->Smooth(badNb);
improved = ( badNb < oldBadNb );
//================================================================================
bool _SolidData::GetShapeEdges(const TGeomID shapeID,
- size_t & edgesEnd,
+ size_t & iEdgesEnd,
int* iBeg,
int* iEnd ) const
{
int beg = 0, end = 0;
- for ( edgesEnd = 0; edgesEnd < _endEdgeOnShape.size(); ++edgesEnd )
+ for ( iEdgesEnd = 0; iEdgesEnd < _endEdgeOnShape.size(); ++iEdgesEnd )
{
- end = _endEdgeOnShape[ edgesEnd ];
+ end = _endEdgeOnShape[ iEdgesEnd ];
TGeomID sID = _edges[ beg ]->_nodes[0]->getshapeId();
if ( sID == shapeID )
{
{
if ( _len - len > -1e-6 )
{
- _pos.push_back( _pos.back() );
+ //_pos.push_back( _pos.back() );
return;
}
// Create intermediate nodes on each _LayerEdge
+ int iS = 0, iEnd = data._endEdgeOnShape[ iS ];
+
for ( size_t i = 0; i < data._edges.size(); ++i )
{
_LayerEdge& edge = *data._edges[i];
if ( edge._nodes.size() < 2 )
continue; // on _noShrinkShapes
+ // get parameters of layers for the edge
+ if ( i == iEnd )
+ iEnd = data._endEdgeOnShape[ ++iS ];
+ const AverageHyp& hyp = data._hypOnShape[ iS ];
+
// get accumulated length of segments
vector< double > segLen( edge._pos.size() );
segLen[0] = 0.0;
const SMDS_MeshNode* tgtNode = edge._nodes.back();
if ( edge._nodes.size() == 2 )
{
- edge._nodes.resize( data._hyp->GetNumberLayers() + 1, 0 );
+ edge._nodes.resize( hyp.GetNumberLayers() + 1, 0 );
edge._nodes[1] = 0;
edge._nodes.back() = tgtNode;
}
// calculate height of the first layer
double h0;
const double T = segLen.back(); //data._hyp.GetTotalThickness();
- const double f = data._hyp->GetStretchFactor();
- const int N = data._hyp->GetNumberLayers();
+ const double f = hyp.GetStretchFactor();
+ const int N = hyp.GetNumberLayers();
const double fPowN = pow( f, N );
if ( fPowN - 1 <= numeric_limits<double>::min() )
h0 = T / N;
TopExp_Explorer exp( data._solid, TopAbs_FACE );
for ( ; exp.More(); exp.Next() )
{
- if ( data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
+ const TGeomID faceID = getMeshDS()->ShapeToIndex( exp.Current() );
+ if ( data._ignoreFaceIds.count( faceID ))
continue;
- SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() );
- SMDS_ElemIteratorPtr fIt = fSubM->GetElements();
+ const bool isReversedFace = data._reversedFaceIds.count( faceID );
+ SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() );
+ SMDS_ElemIteratorPtr fIt = fSubM->GetElements();
while ( fIt->more() )
{
const SMDS_MeshElement* face = fIt->next();
for ( int iN = 0; iN < nbNodes; ++iN )
{
const SMDS_MeshNode* n = nIt->next();
- nnVec[ iN ] = & data._n2eMap[ n ]->_nodes;
- if ( nnVec[ iN ]->size() < 2 )
+ const int i = isReversedFace ? nbNodes-1-iN : iN;
+ nnVec[ i ] = & data._n2eMap[ n ]->_nodes;
+ if ( nnVec[ i ]->size() < 2 )
degenEdgeInd.insert( iN );
else
- nbZ = nnVec[ iN ]->size();
+ nbZ = nnVec[ i ]->size();
if ( helper.HasDegeneratedEdges() )
- nnSet.insert( nnVec[ iN ]);
+ nnSet.insert( nnVec[ i ]);
}
if ( nbZ == 0 )
continue;
{
SMESH_MesherHelper helper( *_mesh );
+ vector< const SMDS_MeshNode* > faceNodes;
+
for ( size_t i = 0; i < _sdVec.size(); ++i )
{
_SolidData& data = _sdVec[i];
if ( nbSharedPyram > 1 )
continue; // not free border of the pyramid
- if ( getMeshDS()->FindFace( ledges[0]->_nodes[0], ledges[0]->_nodes[1],
- ledges[1]->_nodes[0], ledges[1]->_nodes[1]))
+ faceNodes.clear();
+ faceNodes.push_back( ledges[0]->_nodes[0] );
+ faceNodes.push_back( ledges[1]->_nodes[0] );
+ if ( ledges[0]->_nodes.size() > 1 ) faceNodes.push_back( ledges[0]->_nodes[1] );
+ if ( ledges[1]->_nodes.size() > 1 ) faceNodes.push_back( ledges[1]->_nodes[1] );
+
+ if ( getMeshDS()->FindElement( faceNodes, SMDSAbs_Face, /*noMedium=*/true))
continue; // faces already created
}
for ( ++u2n; u2n != u2nodes.end(); ++u2n )