X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ViscousLayers.cxx;h=b2f6855b01c065ec0245b23b6a550c94e1552f25;hp=a855c24a35a5a3182d9c9bcdc3f96b054ef9a581;hb=050aa87698efb9b123985f0bd69aab20ad42c728;hpb=1067ffa6e7e5c394e3a1b17219d8b355a57607cd diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index a855c24a3..b2f6855b0 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -1,9 +1,9 @@ -// Copyright (C) 2007-2012 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2014 CEA/DEN, EDF R&D, OPEN CASCADE // // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public // License as published by the Free Software Foundation; either -// version 2.1 of the License. +// version 2.1 of the License, or (at your option) any later version. // // This library is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of @@ -35,15 +35,18 @@ #include "SMESH_ControlsDef.hxx" #include "SMESH_Gen.hxx" #include "SMESH_Group.hxx" +#include "SMESH_HypoFilter.hxx" #include "SMESH_Mesh.hxx" +#include "SMESH_MeshAlgos.hxx" #include "SMESH_MesherHelper.hxx" #include "SMESH_ProxyMesh.hxx" #include "SMESH_subMesh.hxx" #include "SMESH_subMeshEventListener.hxx" - -#include "utilities.h" +#include "StdMeshers_FaceSide.hxx" #include +#include +#include #include #include #include @@ -53,16 +56,19 @@ #include #include #include +#include #include #include #include #include #include #include +#include #include #include #include #include +#include #include #include #include @@ -224,8 +230,11 @@ namespace VISCOUS_3D struct _Simplex { const SMDS_MeshNode *_nPrev, *_nNext; // nodes on a smoothed mesh surface - _Simplex(const SMDS_MeshNode* nPrev=0, const SMDS_MeshNode* nNext=0) - : _nPrev(nPrev), _nNext(nNext) {} + const SMDS_MeshNode *_nOpp; // in 2D case, a node opposite to a smoothed node in QUAD + _Simplex(const SMDS_MeshNode* nPrev=0, + const SMDS_MeshNode* nNext=0, + const SMDS_MeshNode* nOpp=0) + : _nPrev(nPrev), _nNext(nNext), _nOpp(nOpp) {} bool IsForward(const SMDS_MeshNode* nSrc, const gp_XYZ* pntTgt) const { const double M[3][3] = @@ -317,7 +326,7 @@ namespace VISCOUS_3D gp_XYZ _normal; // to solid surface vector _pos; // points computed during inflation - double _len; // length achived with the last step + double _len; // length achived with the last inflation step double _cosin; // of angle (_normal ^ surface) double _lenFactor; // to compute _len taking _cosin into account @@ -356,7 +365,7 @@ namespace VISCOUS_3D const double& epsilon) const; gp_Ax1 LastSegment(double& segLen) const; bool IsOnEdge() const { return _2neibors; } - void Copy( _LayerEdge& other, SMESH_MesherHelper& helper ); + gp_XYZ Copy( _LayerEdge& other, SMESH_MesherHelper& helper ); void SetCosin( double cosin ); }; struct _LayerEdgeCmp @@ -379,21 +388,27 @@ namespace VISCOUS_3D { TopoDS_Shape _solid; const StdMeshers_ViscousLayers* _hyp; + TopoDS_Shape _hypShape; _MeshOfSolid* _proxyMesh; set _reversedFaceIds; + set _ignoreFaceIds; double _stepSize, _stepSizeCoeff; const SMDS_MeshNode* _stepSizeNodes[2]; TNode2Edge _n2eMap; + // 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 vector< _LayerEdge* > _edges; + // edges on EDGE's with null _sWOL, whose _simplices are used to stop inflation + vector< _LayerEdge* > _simplexTestEdges; - // key: an id of shape (EDGE or VERTEX) shared by a FACE with - // layers and a FACE w/o layers + // key: an id of shape (EDGE or VERTEX) shared by a FACE with + // layers and a FACE w/o layers // value: the shape (FACE or EDGE) to shrink mesh on. - // _LayerEdge's basing on nodes on key shape are inflated along the value shape + // _LayerEdge's basing on nodes on key shape are inflated along the value shape map< TGeomID, TopoDS_Shape > _shrinkShape2Shape; // FACE's WOL, srink on which is forbiden due to algo on the adjacent SOLID @@ -411,7 +426,9 @@ namespace VISCOUS_3D _SolidData(const TopoDS_Shape& s=TopoDS_Shape(), const StdMeshers_ViscousLayers* h=0, - _MeshOfSolid* m=0) :_solid(s), _hyp(h), _proxyMesh(m) {} + const TopoDS_Shape& hs=TopoDS_Shape(), + _MeshOfSolid* m=0) + :_solid(s), _hyp(h), _hypShape(hs), _proxyMesh(m) {} ~_SolidData(); Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge& E, @@ -431,12 +448,18 @@ namespace VISCOUS_3D //vector _nodesAround; vector<_Simplex> _simplices; // for quality check + enum SmoothType { LAPLACIAN, CENTROIDAL, ANGULAR, TFI }; + bool Smooth(int& badNb, Handle(Geom_Surface)& surface, SMESH_MesherHelper& helper, const double refSign, - bool isCentroidal, + SmoothType how, bool set3D); + + gp_XY computeAngularPos(vector& uv, + const gp_XY& uvToFix, + const double refSign ); }; //-------------------------------------------------------------------------------- /*! @@ -463,6 +486,9 @@ namespace VISCOUS_3D bool makeLayer(_SolidData& data); bool setEdgeData(_LayerEdge& edge, const set& subIds, SMESH_MesherHelper& helper, _SolidData& data); + gp_XYZ getWeigthedNormal( const SMDS_MeshNode* n, + std::pair< TGeomID, gp_XYZ > fId2Normal[], + const int nbFaces ); bool findNeiborsOnEdge(const _LayerEdge* edge, const SMDS_MeshNode*& n1, const SMDS_MeshNode*& n2, @@ -471,6 +497,8 @@ namespace VISCOUS_3D const set& ingnoreShapes, const _SolidData* dataToCheckOri = 0, const bool toSort = false); + void findSimplexTestEdges( _SolidData& data, + vector< vector<_LayerEdge*> >& edgesByGeom); bool sortEdges( _SolidData& data, vector< vector<_LayerEdge*> >& edgesByGeom); void limitStepSize( _SolidData& data, @@ -491,7 +519,11 @@ namespace VISCOUS_3D bool prepareEdgeToShrink( _LayerEdge& edge, const TopoDS_Face& F, SMESH_MesherHelper& helper, const SMESHDS_SubMesh* faceSubMesh ); - void fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& helper); + void fixBadFaces(const TopoDS_Face& F, + SMESH_MesherHelper& helper, + const bool is2D, + const int step, + set * involvedNodes=NULL); bool addBoundaryElements(); bool error( const string& text, int solidID=-1 ); @@ -504,7 +536,6 @@ namespace VISCOUS_3D SMESH_ComputeErrorPtr _error; vector< _SolidData > _sdVec; - set _ignoreShapeIds; int _tmpFaceID; }; //-------------------------------------------------------------------------------- @@ -540,7 +571,7 @@ namespace VISCOUS_3D virtual vtkIdType GetVtkType() const { return -1; } virtual SMDSAbs_EntityType GetEntityType() const { return SMDSEntity_Last; } virtual SMDSAbs_GeometryType GetGeomType() const { return SMDSGeom_TRIANGLE; } -virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const + virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType) const { return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));} }; //-------------------------------------------------------------------------------- @@ -559,6 +590,44 @@ virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const _nn[3]=_le2->_nodes[0]; } }; + //-------------------------------------------------------------------------------- + /*! + * \brief Retriever of node coordinates either directly of from a surface by node UV. + * \warning Location of a surface is ignored + */ + struct NodeCoordHelper + { + SMESH_MesherHelper& _helper; + const TopoDS_Face& _face; + Handle(Geom_Surface) _surface; + gp_XYZ (NodeCoordHelper::* _fun)(const SMDS_MeshNode* n) const; + + NodeCoordHelper(const TopoDS_Face& F, SMESH_MesherHelper& helper, bool is2D) + : _helper( helper ), _face( F ) + { + if ( is2D ) + { + TopLoc_Location loc; + _surface = BRep_Tool::Surface( _face, loc ); + } + if ( _surface.IsNull() ) + _fun = & NodeCoordHelper::direct; + else + _fun = & NodeCoordHelper::byUV; + } + gp_XYZ operator()(const SMDS_MeshNode* n) const { return (this->*_fun)( n ); } + + private: + gp_XYZ direct(const SMDS_MeshNode* n) const + { + return SMESH_TNodeXYZ( n ); + } + gp_XYZ byUV (const SMDS_MeshNode* n) const + { + gp_XY uv = _helper.GetNodeUV( _face, n ); + return _surface->Value( uv.X(), uv.Y() ).XYZ(); + } + }; } // namespace VISCOUS_3D //================================================================================ @@ -566,20 +635,17 @@ virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const // StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH_Gen* gen) :SMESH_Hypothesis(hypId, studyId, gen), - _nbLayers(1), _thickness(1), _stretchFactor(1) + _isToIgnoreShapes(1), _nbLayers(1), _thickness(1), _stretchFactor(1) { _name = StdMeshers_ViscousLayers::GetHypType(); _param_algo_dim = -3; // auxiliary hyp used by 3D algos } // -------------------------------------------------------------------------------- -void StdMeshers_ViscousLayers::SetBndShapesToIgnore(const std::vector& faceIds) -{ - if ( faceIds != _ignoreBndShapeIds ) - _ignoreBndShapeIds = faceIds, NotifySubMeshesHypothesisModification(); -} // -------------------------------------------------------------------------------- -bool StdMeshers_ViscousLayers::IsIgnoredShape(const int shapeID) const +void StdMeshers_ViscousLayers::SetBndShapes(const std::vector& faceIds, bool toIgnore) { - return ( find( _ignoreBndShapeIds.begin(), _ignoreBndShapeIds.end(), shapeID ) - != _ignoreBndShapeIds.end() ); + if ( faceIds != _shapeIds ) + _shapeIds = faceIds, NotifySubMeshesHypothesisModification(); + if ( _isToIgnoreShapes != toIgnore ) + _isToIgnoreShapes = toIgnore, NotifySubMeshesHypothesisModification(); } // -------------------------------------------------------------------------------- void StdMeshers_ViscousLayers::SetTotalThickness(double thickness) { @@ -637,17 +703,22 @@ std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save) save << " " << _nbLayers << " " << _thickness << " " << _stretchFactor - << " " << _ignoreBndShapeIds.size(); - for ( unsigned i = 0; i < _ignoreBndShapeIds.size(); ++i ) - save << " " << _ignoreBndShapeIds[i]; + << " " << _shapeIds.size(); + for ( size_t i = 0; i < _shapeIds.size(); ++i ) + save << " " << _shapeIds[i]; + save << " " << !_isToIgnoreShapes; // negate to keep the behavior in old studies. return save; } // -------------------------------------------------------------------------------- std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load) { - int nbFaces, faceID; + int nbFaces, faceID, shapeToTreat; load >> _nbLayers >> _thickness >> _stretchFactor >> nbFaces; - while ( _ignoreBndShapeIds.size() < nbFaces && load >> faceID ) - _ignoreBndShapeIds.push_back( faceID ); + while ( _shapeIds.size() < nbFaces && load >> faceID ) + _shapeIds.push_back( faceID ); + if ( load >> shapeToTreat ) + _isToIgnoreShapes = !shapeToTreat; + else + _isToIgnoreShapes = true; // old behavior return load; } // -------------------------------------------------------------------------------- bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh* theMesh, @@ -735,23 +806,29 @@ namespace // get average dir of edges going fromV gp_XYZ edgeDir; //if ( edges.size() > 1 ) - for ( unsigned i = 0; i < edges.size(); ++i ) - { - edgeDir = getEdgeDir( edges[i], fromV ); - double size2 = edgeDir.SquareModulus(); - if ( size2 > numeric_limits::min() ) - edgeDir /= sqrt( size2 ); - else - ok = false; - dir += edgeDir; - } + for ( size_t i = 0; i < edges.size(); ++i ) + { + edgeDir = getEdgeDir( edges[i], fromV ); + double size2 = edgeDir.SquareModulus(); + if ( size2 > numeric_limits::min() ) + edgeDir /= sqrt( size2 ); + else + ok = false; + dir += edgeDir; + } gp_XYZ fromEdgeDir = getFaceDir( F, edges[0], node, helper, ok ); - if ( edges.size() == 1 ) - dir = fromEdgeDir; - else if ( dir.SquareModulus() < 0.1 ) // ~< 20 degrees - dir = fromEdgeDir + getFaceDir( F, edges[1], node, helper, ok ); - else if ( dir * fromEdgeDir < 0 ) - dir *= -1; + try { + if ( edges.size() == 1 ) + dir = fromEdgeDir; + else if ( dir.SquareModulus() < 0.1 ) // ~< 20 degrees + dir = fromEdgeDir.Normalized() + getFaceDir( F, edges[1], node, helper, ok ).Normalized(); + else if ( dir * fromEdgeDir < 0 ) + dir *= -1; + } + catch ( Standard_Failure ) + { + ok = false; + } if ( ok ) { //dir /= edges.size(); @@ -770,13 +847,15 @@ namespace bool isConcave( const TopoDS_Face& F, SMESH_MesherHelper& helper ) { + // if ( helper.Count( F, TopAbs_WIRE, /*useMap=*/false) > 1 ) + // return true; gp_Vec2d drv1, drv2; gp_Pnt2d p; TopExp_Explorer eExp( F.Oriented( TopAbs_FORWARD ), TopAbs_EDGE ); for ( ; eExp.More(); eExp.Next() ) { const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() ); - if ( BRep_Tool::Degenerated( E )) continue; + if ( SMESH_Algo::isDegenerated( E )) continue; // check if 2D curve is concave BRepAdaptor_Curve2d curve( E, F ); const int nbIntervals = curve.NbIntervals( GeomAbs_C2 ); @@ -791,7 +870,7 @@ namespace double cross = drv2 ^ drv1; if ( E.Orientation() == TopAbs_REVERSED ) cross = -cross; - isConvex = ( cross < 1e-9 ); + isConvex = ( cross > 0.1 ); //-1e-9 ); } // check if concavity is strong enough to care about it //const double maxAngle = 5 * Standard_PI180; @@ -821,6 +900,26 @@ namespace // } } } + // check angles at VERTEXes + TError error; + TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, *helper.GetMesh(), 0, error ); + for ( size_t iW = 0; iW < wires.size(); ++iW ) + { + const int nbEdges = wires[iW]->NbEdges(); + if ( nbEdges < 2 && SMESH_Algo::isDegenerated( wires[iW]->Edge(0))) + continue; + for ( int iE1 = 0; iE1 < nbEdges; ++iE1 ) + { + if ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE1 ))) continue; + int iE2 = ( iE1 + 1 ) % nbEdges; + while ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE2 ))) + iE2 = ( iE2 + 1 ) % nbEdges; + double angle = helper.GetAngle( wires[iW]->Edge( iE1 ), + wires[iW]->Edge( iE2 ), F ); + if ( angle < -5. * M_PI / 180. ) + return true; + } + } return false; } //-------------------------------------------------------------------------------- @@ -832,13 +931,15 @@ namespace const char* fname = "/tmp/viscous.py"; cout << "execfile('"< notSupportAlgos; notSupportAlgos.insert("Hexa_3D"); - for ( unsigned i = 0; i < _sdVec.size(); ++i ) + for ( size_t i = 0; i < _sdVec.size(); ++i ) { TopTools_MapOfShape noShrinkVertices; map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin(); @@ -1152,7 +1300,7 @@ bool _ViscousBuilder::findFacesWithLayers() SMESH_Algo* algo = _mesh->GetGen()->GetAlgo( *_mesh, *solid ); if ( !algo || !notSupportAlgos.count( algo->GetName() )) continue; notShrinkFace = true; - for ( unsigned j = 0; j < _sdVec.size(); ++j ) + for ( size_t j = 0; j < _sdVec.size(); ++j ) { if ( _sdVec[j]._solid.IsSame( *solid ) ) if ( _sdVec[j]._shrinkShape2Shape.count( edgeID )) @@ -1185,10 +1333,10 @@ bool _ViscousBuilder::findFacesWithLayers() } } } - + // Find the SHAPE along which to inflate _LayerEdge based on VERTEX - for ( unsigned i = 0; i < _sdVec.size(); ++i ) + for ( size_t i = 0; i < _sdVec.size(); ++i ) { shapes.Clear(); TopExp::MapShapes(_sdVec[i]._solid, TopAbs_VERTEX, shapes); @@ -1202,11 +1350,12 @@ bool _ViscousBuilder::findFacesWithLayers() while ( fIt->more()) { const TopoDS_Shape* f = fIt->next(); - const int fID = getMeshDS()->ShapeToIndex( *f ); if ( helper.IsSubShape( *f, _sdVec[i]._solid ) ) { totalNbFaces++; - if ( _ignoreShapeIds.count ( fID ) && ! _sdVec[i]._noShrinkFaces.count( fID )) + const int fID = getMeshDS()->ShapeToIndex( *f ); + if ( _sdVec[i]._ignoreFaceIds.count ( fID ) && + !_sdVec[i]._noShrinkFaces.count( fID )) facesWOL.push_back( *f ); } } @@ -1216,48 +1365,61 @@ bool _ViscousBuilder::findFacesWithLayers() switch ( facesWOL.size() ) { case 1: + { + helper.SetSubShape( facesWOL[0] ); + if ( helper.IsRealSeam( vInd )) // inflate along a seam edge? { - helper.SetSubShape( facesWOL[0] ); - if ( helper.IsRealSeam( vInd )) // inflate along a seam edge? + TopoDS_Shape seamEdge; + PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE); + while ( eIt->more() && seamEdge.IsNull() ) { - TopoDS_Shape seamEdge; - PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE); - while ( eIt->more() && seamEdge.IsNull() ) - { - const TopoDS_Shape* e = eIt->next(); - if ( helper.IsRealSeam( *e ) ) - seamEdge = *e; - } - if ( !seamEdge.IsNull() ) - { - _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, seamEdge )); - break; - } + const TopoDS_Shape* e = eIt->next(); + if ( helper.IsRealSeam( *e ) ) + seamEdge = *e; + } + if ( !seamEdge.IsNull() ) + { + _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, seamEdge )); + break; } - _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] )); - break; } + _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] )); + break; + } case 2: + { + // find an edge shared by 2 faces + PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE); + while ( eIt->more()) { - // find an edge shared by 2 faces - PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE); - while ( eIt->more()) + const TopoDS_Shape* e = eIt->next(); + if ( helper.IsSubShape( *e, facesWOL[0]) && + helper.IsSubShape( *e, facesWOL[1])) { - const TopoDS_Shape* e = eIt->next(); - if ( helper.IsSubShape( *e, facesWOL[0]) && - helper.IsSubShape( *e, facesWOL[1])) - { - _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, *e )); break; - } + _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, *e )); break; } - break; } + break; + } default: return error("Not yet supported case", _sdVec[i]._index); } } } + // add FACEs of other SOLIDs to _ignoreFaceIds + for ( size_t i = 0; i < _sdVec.size(); ++i ) + { + shapes.Clear(); + TopExp::MapShapes(_sdVec[i]._solid, TopAbs_FACE, shapes); + + for ( exp.Init( _mesh->GetShapeToMesh(), TopAbs_FACE ); exp.More(); exp.Next() ) + { + if ( !shapes.Contains( exp.Current() )) + _sdVec[i]._ignoreFaceIds.insert( getMeshDS()->ShapeToIndex( exp.Current() )); + } + } + return true; } @@ -1274,10 +1436,10 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) subIds = data._noShrinkFaces; TopExp_Explorer exp( data._solid, TopAbs_FACE ); for ( ; exp.More(); exp.Next() ) - if ( ! _ignoreShapeIds.count( getMeshDS()->ShapeToIndex( exp.Current() ))) { SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() ); - faceIds.insert( fSubM->GetId() ); + if ( ! data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() ))) + faceIds.insert( fSubM->GetId() ); SMESH_subMeshIteratorPtr subIt = fSubM->getDependsOnIterator(/*includeSelf=*/true, /*complexShapeFirst=*/false); while ( subIt->more() ) @@ -1285,20 +1447,19 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) } // make a map to find new nodes on sub-shapes shared with other SOLID - map< TGeomID, TNode2Edge* > s2neMap; map< TGeomID, TNode2Edge* >::iterator s2ne; map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin(); for (; s2s != data._shrinkShape2Shape.end(); ++s2s ) { TGeomID shapeInd = s2s->first; - for ( unsigned i = 0; i < _sdVec.size(); ++i ) + for ( size_t i = 0; i < _sdVec.size(); ++i ) { if ( _sdVec[i]._index == data._index ) continue; map< TGeomID, TopoDS_Shape >::iterator s2s2 = _sdVec[i]._shrinkShape2Shape.find( shapeInd ); if ( s2s2 != _sdVec[i]._shrinkShape2Shape.end() && *s2s == *s2s2 && !_sdVec[i]._n2eMap.empty() ) { - s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap )); + data._s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap )); break; } } @@ -1350,21 +1511,23 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) const int shapeID = n->getshapeId(); edgesByGeom[ shapeID ].push_back( edge ); + SMESH_TNodeXYZ xyz( n ); + // set edge data or find already refined _LayerEdge and get data from it if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE && - ( s2ne = s2neMap.find( shapeID )) != s2neMap.end() && + ( s2ne = data._s2neMap.find( shapeID )) != data._s2neMap.end() && ( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end()) { _LayerEdge* foundEdge = (*n2e2).second; - edge->Copy( *foundEdge, helper ); - // location of the last node is modified but we can restore - // it by node position on _sWOL stored by the node + gp_XYZ lastPos = edge->Copy( *foundEdge, helper ); + foundEdge->_pos.push_back( lastPos ); + // location of the last node is modified and we restore it by foundEdge->_pos.back() const_cast< SMDS_MeshNode* > - ( edge->_nodes.back() )->setXYZ( n->X(), n->Y(), n->Z() ); + ( edge->_nodes.back() )->setXYZ( xyz.X(), xyz.Y(), xyz.Z() ); } else { - edge->_nodes.push_back( helper.AddNode( n->X(), n->Y(), n->Z() )); + edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() )); if ( !setEdgeData( *edge, subIds, helper, data )) return false; } @@ -1391,14 +1554,16 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) if ( data._stepSize < 1. ) data._epsilon *= data._stepSize; - // Put _LayerEdge's into the vector data._edges + // fill data._simplexTestEdges + findSimplexTestEdges( data, edgesByGeom ); + // Put _LayerEdge's into the vector data._edges if ( !sortEdges( data, edgesByGeom )) return false; - // Set target nodes into _Simplex and _2NearEdges + // Set target nodes into _Simplex and _2NearEdges of _LayerEdge's TNode2Edge::iterator n2e; - for ( unsigned i = 0; i < data._edges.size(); ++i ) + for ( size_t i = 0; i < data._edges.size(); ++i ) { if ( data._edges[i]->IsOnEdge()) for ( int j = 0; j < 2; ++j ) @@ -1411,13 +1576,13 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) n = (*n2e).second->_nodes.back(); data._edges[i]->_2neibors->_edges[j] = n2e->second; } - else - for ( unsigned j = 0; j < data._edges[i]->_simplices.size(); ++j ) - { - _Simplex& s = data._edges[i]->_simplices[j]; - s._nNext = data._n2eMap[ s._nNext ]->_nodes.back(); - s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back(); - } + //else + for ( size_t j = 0; j < data._edges[i]->_simplices.size(); ++j ) + { + _Simplex& s = data._edges[i]->_simplices[j]; + s._nNext = data._n2eMap[ s._nNext ]->_nodes.back(); + s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back(); + } } dumpFunctionEnd(); @@ -1465,7 +1630,7 @@ void _ViscousBuilder::limitStepSize( _SolidData& data, */ //================================================================================ -void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize) +void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize ) { if ( minSize < data._stepSize ) { @@ -1479,6 +1644,96 @@ void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize) } } +//================================================================================ +/*! + * Fill data._simplexTestEdges. These _LayerEdge's are used to stop inflation + * in the case where there are no _LayerEdge's on a curved convex FACE, + * as e.g. on a fillet surface with no internal nodes - issue 22580, + * so that collision of viscous internal faces is not detected by check of + * intersection of _LayerEdge's with the viscous internal faces. + */ +//================================================================================ + +void _ViscousBuilder::findSimplexTestEdges( _SolidData& data, + vector< vector<_LayerEdge*> >& edgesByGeom) +{ + data._simplexTestEdges.clear(); + + SMESH_MesherHelper helper( *_mesh ); + + vector< vector<_LayerEdge*> * > ledgesOnEdges; + set< const SMDS_MeshNode* > usedNodes; + + const double minCurvature = 1. / data._hyp->GetTotalThickness(); + + for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS ) + { + // look for a FACE with layers and w/o _LayerEdge's + const vector<_LayerEdge*>& eS = edgesByGeom[iS]; + if ( !eS.empty() ) continue; + const TopoDS_Shape& S = getMeshDS()->IndexToShape( iS ); + if ( S.IsNull() || S.ShapeType() != TopAbs_FACE ) continue; + if ( data._ignoreFaceIds.count( iS )) continue; + + const TopoDS_Face& F = TopoDS::Face( S ); + + // look for _LayerEdge's on EDGEs with null _sWOL + ledgesOnEdges.clear(); + TopExp_Explorer eExp( F, TopAbs_EDGE ); + for ( ; eExp.More(); eExp.Next() ) + { + TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() ); + vector<_LayerEdge*>& eE = edgesByGeom[iE]; + if ( !eE.empty() && eE[0]->_sWOL.IsNull() ) + ledgesOnEdges.push_back( & eE ); + } + if ( ledgesOnEdges.empty() ) continue; + + // check FACE convexity + const _LayerEdge* le = ledgesOnEdges[0]->back(); + gp_XY uv = helper.GetNodeUV( F, le->_nodes[0] ); + BRepAdaptor_Surface surf( F ); + BRepLProp_SLProps surfProp( surf, uv.X(), uv.Y(), 2, 1e-6 ); + if ( !surfProp.IsCurvatureDefined() ) + continue; + double surfCurvature = Max( Abs( surfProp.MaxCurvature() ), + Abs( surfProp.MinCurvature() )); + if ( surfCurvature < minCurvature ) + continue; + gp_Dir minDir, maxDir; + surfProp.CurvatureDirections( maxDir, minDir ); + if ( F.Orientation() == TopAbs_REVERSED ) { + maxDir.Reverse(); minDir.Reverse(); + } + const gp_XYZ& inDir = le->_normal; + if ( inDir * maxDir.XYZ() < 0 && + inDir * minDir.XYZ() < 0 ) + continue; + + limitStepSize( data, 0.9 / surfCurvature ); + + // add _simplices to the _LayerEdge's + for ( size_t iE = 0; iE < ledgesOnEdges.size(); ++iE ) + { + const vector<_LayerEdge*>& ledges = *ledgesOnEdges[iE]; + for ( size_t iLE = 0; iLE < ledges.size(); ++iLE ) + { + _LayerEdge* ledge = ledges[iLE]; + const SMDS_MeshNode* srcNode = ledge->_nodes[0]; + if ( !usedNodes.insert( srcNode ).second ) continue; + + getSimplices( srcNode, ledge->_simplices, data._ignoreFaceIds, &data ); + for ( size_t i = 0; i < ledge->_simplices.size(); ++i ) + { + usedNodes.insert( ledge->_simplices[i]._nPrev ); + usedNodes.insert( ledge->_simplices[i]._nNext ); + } + data._simplexTestEdges.push_back( ledge ); + } + } + } +} + //================================================================================ /*! * \brief Separate shapes (and _LayerEdge's on them) to smooth from the rest ones @@ -1496,11 +1751,11 @@ bool _ViscousBuilder::sortEdges( _SolidData& data, SMESH_MesherHelper helper( *_mesh ); bool ok = true; - for ( unsigned iS = 0; iS < edgesByGeom.size(); ++iS ) + for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS ) { vector<_LayerEdge*>& eS = edgesByGeom[iS]; if ( eS.empty() ) continue; - TopoDS_Shape S = getMeshDS()->IndexToShape( iS ); + const TopoDS_Shape& S = getMeshDS()->IndexToShape( iS ); bool needSmooth = false; switch ( S.ShapeType() ) { @@ -1540,7 +1795,7 @@ bool _ViscousBuilder::sortEdges( _SolidData& data, if ( eE.empty() ) continue; if ( eE[0]->_sWOL.IsNull() ) { - for ( unsigned i = 0; i < eE.size() && !needSmooth; ++i ) + for ( size_t i = 0; i < eE.size() && !needSmooth; ++i ) needSmooth = ( eE[i]->_cosin > 0.1 ); } else @@ -1548,7 +1803,7 @@ bool _ViscousBuilder::sortEdges( _SolidData& data, 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 ( unsigned i = 0; i < eE.size() && !needSmooth; ++i ) + 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 ); @@ -1587,7 +1842,7 @@ bool _ViscousBuilder::sortEdges( _SolidData& data, } // then the rest _LayerEdge's - for ( unsigned iS = 0; iS < edgesByGeom.size(); ++iS ) + for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS ) { vector<_LayerEdge*>& eVec = edgesByGeom[iS]; data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() ); @@ -1673,32 +1928,62 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, set::iterator id = faceIds.begin(); TopoDS_Face F; + std::pair< TGeomID, gp_XYZ > id2Norm[20]; for ( ; id != faceIds.end(); ++id ) { const TopoDS_Shape& s = getMeshDS()->IndexToShape( *id ); if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id )) continue; - totalNbFaces++; - //nbLayerFaces += subIds.count( *id ); F = TopoDS::Face( s ); gp_XY uv = helper.GetNodeUV( F, node, 0, &normOK ); Handle(Geom_Surface) surface = BRep_Tool::Surface( F ); - surface->D1( uv.X(),uv.Y(), p, du,dv ); - geomNorm = du ^ dv; - double size2 = geomNorm.SquareMagnitude(); - if ( size2 > numeric_limits::min() ) - geomNorm /= sqrt( size2 ); - else - normOK = false; + { + gp_Dir normal; + if ( GeomLib::NormEstim( surface, uv, 1e-10, normal ) < 3 ) + { + geomNorm = normal; + normOK = true; + } + else // hard singularity + { + SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face); + while ( fIt->more() ) + { + const SMDS_MeshElement* f = fIt->next(); + if ( editor.FindShape( f ) == *id ) + { + SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) geomNorm.XYZ(), /*normalized=*/false ); + if ( helper.IsReversedSubMesh( F )) + geomNorm.Reverse(); + break; + } + } + double size2 = geomNorm.SquareMagnitude(); + if ( size2 > numeric_limits::min() ) + geomNorm /= sqrt( size2 ); + else + normOK = false; + } + } if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED ) geomNorm.Reverse(); + id2Norm[ totalNbFaces ].first = *id; + id2Norm[ totalNbFaces ].second = geomNorm.XYZ(); + totalNbFaces++; edge._normal += geomNorm.XYZ(); } if ( totalNbFaces == 0 ) return error(SMESH_Comment("Can't get normal to node ") << node->GetID(), data._index); - edge._normal /= totalNbFaces; + if ( totalNbFaces < 3 ) + { + //edge._normal /= totalNbFaces; + } + else + { + edge._normal = getWeigthedNormal( node, id2Norm, totalNbFaces ); + } switch ( posType ) { @@ -1764,9 +2049,9 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, if ( posType == SMDS_TOP_FACE ) { - getSimplices( node, edge._simplices, _ignoreShapeIds, &data ); + getSimplices( node, edge._simplices, data._ignoreFaceIds, &data ); double avgNormProj = 0, avgLen = 0; - for ( unsigned i = 0; i < edge._simplices.size(); ++i ) + for ( size_t i = 0; i < edge._simplices.size(); ++i ) { gp_XYZ vec = edge._pos.back() - SMESH_TNodeXYZ( edge._simplices[i]._nPrev ); avgNormProj += edge._normal * vec; @@ -1800,6 +2085,87 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, return true; } +//================================================================================ +/*! + * \brief Return normal at a node weighted with angles taken by FACEs + * \param [in] n - the node + * \param [in] fId2Normal - FACE ids and normals + * \param [in] nbFaces - nb of FACEs meeting at the node + * \return gp_XYZ - computed normal + */ +//================================================================================ + +gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode* n, + std::pair< TGeomID, gp_XYZ > fId2Normal[], + const int nbFaces ) +{ + gp_XYZ resNorm(0,0,0); + TopoDS_Shape V = SMESH_MesherHelper::GetSubShapeByNode( n, getMeshDS() ); + if ( V.ShapeType() != TopAbs_VERTEX ) + { + for ( int i = 0; i < nbFaces; ++i ) + resNorm += fId2Normal[i].second / nbFaces ; + return resNorm; + } + + double angles[30]; + for ( int i = 0; i < nbFaces; ++i ) + { + const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( fId2Normal[i].first )); + + // look for two EDGEs shared by F and other FACEs within fId2Normal + TopoDS_Edge ee[2]; + int nbE = 0; + PShapeIteratorPtr eIt = SMESH_MesherHelper::GetAncestors( V, *_mesh, TopAbs_EDGE ); + while ( const TopoDS_Shape* E = eIt->next() ) + { + if ( !SMESH_MesherHelper::IsSubShape( *E, F )) + continue; + bool isSharedEdge = false; + for ( int j = 0; j < nbFaces && !isSharedEdge; ++j ) + { + if ( i == j ) continue; + const TopoDS_Shape& otherF = getMeshDS()->IndexToShape( fId2Normal[j].first ); + isSharedEdge = SMESH_MesherHelper::IsSubShape( *E, otherF ); + } + if ( !isSharedEdge ) + continue; + ee[ nbE ] = TopoDS::Edge( *E ); + ee[ nbE ].Orientation( SMESH_MesherHelper::GetSubShapeOri( F, *E )); + if ( ++nbE == 2 ) + break; + } + + // get an angle between the two EDGEs + angles[i] = 0; + if ( nbE < 1 ) continue; + if ( nbE == 1 ) + { + ee[ 1 ] == ee[ 0 ]; + } + else + { + TopoDS_Vertex v10 = SMESH_MesherHelper::IthVertex( 1, ee[ 0 ]); + TopoDS_Vertex v01 = SMESH_MesherHelper::IthVertex( 0, ee[ 1 ]); + if ( !v10.IsSame( v01 )) + std::swap( ee[0], ee[1] ); + } + angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F ); + } + + // compute a weighted normal + double sumAngle = 0; + for ( int i = 0; i < nbFaces; ++i ) + { + angles[i] = ( angles[i] > 2*M_PI ) ? 0 : M_PI - angles[i]; + sumAngle += angles[i]; + } + for ( int i = 0; i < nbFaces; ++i ) + resNorm += angles[i] / sumAngle * fId2Normal[i].second; + + return resNorm; +} + //================================================================================ /*! * \brief Find 2 neigbor nodes of a node on EDGE @@ -1903,7 +2269,7 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1, */ //================================================================================ -void _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper ) +gp_XYZ _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper ) { _nodes = other._nodes; _normal = other._normal; @@ -1915,16 +2281,25 @@ void _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper ) _curvature = 0; std::swap( _curvature, other._curvature ); _2neibors = 0; std::swap( _2neibors, other._2neibors ); + gp_XYZ lastPos( 0,0,0 ); if ( _sWOL.ShapeType() == TopAbs_EDGE ) { double u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes[0] ); _pos.push_back( gp_XYZ( u, 0, 0)); + + u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes.back() ); + lastPos.SetX( u ); } else // TopAbs_FACE { gp_XY uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes[0]); _pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0)); + + uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes.back() ); + lastPos.SetX( uv.X() ); + lastPos.SetY( uv.Y() ); } + return lastPos; } //================================================================================ @@ -1936,7 +2311,7 @@ void _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper ) void _LayerEdge::SetCosin( double cosin ) { _cosin = cosin; - _lenFactor = ( _cosin > 0.1 ) ? 1./sqrt(1-_cosin*_cosin) : 1.0; + _lenFactor = ( Abs( _cosin ) > 0.1 ) ? 1./sqrt(1-_cosin*_cosin) : 1.0; } //================================================================================ @@ -1951,19 +2326,21 @@ void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node, const _SolidData* dataToCheckOri, const bool toSort) { + simplices.clear(); SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face); while ( fIt->more() ) { const SMDS_MeshElement* f = fIt->next(); - const TGeomID shapeInd = f->getshapeId(); + const TGeomID shapeInd = f->getshapeId(); if ( ingnoreShapes.count( shapeInd )) continue; const int nbNodes = f->NbCornerNodes(); - int srcInd = f->GetNodeIndex( node ); + const int srcInd = f->GetNodeIndex( node ); const SMDS_MeshNode* nPrev = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd-1, nbNodes )); const SMDS_MeshNode* nNext = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+1, nbNodes )); + const SMDS_MeshNode* nOpp = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+2, nbNodes )); if ( dataToCheckOri && dataToCheckOri->_reversedFaceIds.count( shapeInd )) std::swap( nPrev, nNext ); - simplices.push_back( _Simplex( nPrev, nNext )); + simplices.push_back( _Simplex( nPrev, nNext, nOpp )); } if ( toSort ) @@ -1995,7 +2372,7 @@ void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node, void _ViscousBuilder::makeGroupOfLE() { #ifdef _DEBUG_ - for ( unsigned i = 0 ; i < _sdVec.size(); ++i ) + for ( size_t i = 0 ; i < _sdVec.size(); ++i ) { if ( _sdVec[i]._edges.empty() ) continue; // string name = SMESH_Comment("_LayerEdge's_") << i; @@ -2005,10 +2382,10 @@ void _ViscousBuilder::makeGroupOfLE() // SMESHDS_Mesh* mDS = _mesh->GetMeshDS(); dumpFunction( SMESH_Comment("make_LayerEdge_") << i ); - for ( unsigned j = 0 ; j < _sdVec[i]._edges.size(); ++j ) + for ( size_t j = 0 ; j < _sdVec[i]._edges.size(); ++j ) { _LayerEdge* le = _sdVec[i]._edges[j]; - for ( unsigned iN = 1; iN < le->_nodes.size(); ++iN ) + for ( size_t iN = 1; iN < le->_nodes.size(); ++iN ) dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <_nodes[iN-1]->GetID() << ", " << le->_nodes[iN]->GetID() <<"])"); //gDS->SMDSGroup().Add( mDS->AddEdge( le->_nodes[iN-1], le->_nodes[iN])); @@ -2016,7 +2393,7 @@ void _ViscousBuilder::makeGroupOfLE() dumpFunctionEnd(); dumpFunction( SMESH_Comment("makeNormals") << i ); - for ( unsigned j = 0 ; j < _sdVec[i]._edges.size(); ++j ) + for ( size_t j = 0 ; j < _sdVec[i]._edges.size(); ++j ) { _LayerEdge& edge = *_sdVec[i]._edges[j]; SMESH_TNodeXYZ nXYZ( edge._nodes[0] ); @@ -2067,14 +2444,14 @@ bool _ViscousBuilder::inflate(_SolidData& data) // Limit inflation step size by geometry size found by itersecting // normals of _LayerEdge's with mesh faces double geomSize = Precision::Infinite(), intersecDist; - SMESH_MeshEditor editor( _mesh ); auto_ptr searcher - ( editor.GetElementSearcher( data._proxyMesh->GetFaces( data._solid )) ); - for ( unsigned i = 0; i < data._edges.size(); ++i ) + ( 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 ) + if ( geomSize > intersecDist && intersecDist > 0 ) geomSize = intersecDist; } if ( data._stepSize > 0.3 * geomSize ) @@ -2105,7 +2482,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) // Elongate _LayerEdge's dumpFunction(SMESH_Comment("inflate")<SetNewLength( curThick, helper ); } @@ -2121,7 +2498,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) if ( nbSteps > 0 ) { dumpFunction(SMESH_Comment("invalidate")<InvalidateStep( nbSteps+1 ); } @@ -2133,7 +2510,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) // Evaluate achieved thickness avgThick = 0; - for ( unsigned i = 0; i < data._edges.size(); ++i ) + for ( size_t i = 0; i < data._edges.size(); ++i ) avgThick += data._edges[i]->_len; avgThick /= data._edges.size(); #ifdef __myDEBUG @@ -2181,7 +2558,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, TopoDS_Face F; int iBeg, iEnd = 0; - for ( unsigned iS = 0; iS < data._endEdgeToSmooth.size(); ++iS ) + for ( size_t iS = 0; iS < data._endEdgeToSmooth.size(); ++iS ) { iBeg = iEnd; iEnd = data._endEdgeToSmooth[ iS ]; @@ -2226,8 +2603,8 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, else { // smooth on FACE's - int step = 0, badNb = 0; moved = true; - while (( ++step <= 5 && moved ) || improved ) + int step = 0, stepLimit = 5, badNb = 0; moved = true; + while (( ++step <= stepLimit && moved ) || improved ) { dumpFunction(SMESH_Comment("smooth")<Smooth(badNb); improved = ( badNb < oldBadNb ); + // issue 22576. no bad faces but still there are intersections to fix + if ( improved && badNb == 0 ) + stepLimit = step + 3; + dumpFunctionEnd(); } if ( badNb > 0 ) @@ -2247,7 +2628,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, { _LayerEdge* edge = data._edges[i]; SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() ); - for ( unsigned j = 0; j < edge->_simplices.size(); ++j ) + for ( size_t j = 0; j < edge->_simplices.size(); ++j ) if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ )) { cout << "Bad simplex ( " << edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID() @@ -2262,12 +2643,30 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, } } // loop on shapes to smooth + // Check orientation of simplices of _simplexTestEdges + for ( size_t i = 0; i < data._simplexTestEdges.size(); ++i ) + { + const _LayerEdge* edge = data._simplexTestEdges[i]; + SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() ); + for ( size_t j = 0; j < edge->_simplices.size(); ++j ) + if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ )) + { +#ifdef __myDEBUG + cout << "Bad simplex of _simplexTestEdges (" + << " "<< edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID() + << " "<< edge->_simplices[j]._nPrev->GetID() + << " "<< edge->_simplices[j]._nNext->GetID() << " )" << endl; +#endif + return false; + } + } + // Check if the last segments of _LayerEdge intersects 2D elements; // checked elements are either temporary faces or faces on surfaces w/o the layers - SMESH_MeshEditor editor( _mesh ); auto_ptr searcher - ( editor.GetElementSearcher( data._proxyMesh->GetFaces( data._solid )) ); + ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), + data._proxyMesh->GetFaces( data._solid )) ); distToIntersection = Precision::Infinite(); double dist; @@ -2276,7 +2675,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, const SMDS_MeshElement* closestFace = 0; int iLE = 0; #endif - for ( unsigned i = 0; i < data._edges.size(); ++i ) + for ( size_t i = 0; i < data._edges.size(); ++i ) { if ( data._edges[i]->FindIntersection( *searcher, dist, data._epsilon, &intFace )) return false; @@ -2522,6 +2921,10 @@ bool _ViscousBuilder::smoothAnalyticEdge( _SolidData& data, if ( F.IsNull() ) // 3D { + if ( data._edges[iFrom]->_2neibors->_nodes[0] == + data._edges[iTo-1]->_2neibors->_nodes[1] ) + return true; // closed EDGE - nothing to do + return false; // TODO ??? } else // 2D @@ -2583,7 +2986,7 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, vector< const SMDS_MeshNode*> nodes(4); // of a tmp mesh face dumpFunction(SMESH_Comment("makeTmpFacesOnEdges")<IsOnEdge() || !edge->_sWOL.IsNull() ) continue; @@ -2600,7 +3003,7 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, } // look for a _LayerEdge containg tgt2 // _LayerEdge* neiborEdge = 0; -// unsigned di = 0; // check _edges[i+di] and _edges[i-di] +// size_t di = 0; // check _edges[i+di] and _edges[i-di] // while ( !neiborEdge && ++di <= data._edges.size() ) // { // if ( i+di < data._edges.size() && data._edges[i+di]->_nodes.back() == tgt2 ) @@ -2627,10 +3030,10 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, // 1) to find and fix intersection // 2) to check that no new intersection appears as result of 1) - SMESH_MeshEditor editor( _mesh ); SMDS_ElemIteratorPtr fIt( new SMDS_ElementVectorIterator( tmpFaces.begin(), tmpFaces.end())); - auto_ptr searcher ( editor.GetElementSearcher( fIt )); + auto_ptr searcher + ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), fIt )); // 1) Find intersections double dist; @@ -2639,7 +3042,7 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, TLEdge2LEdgeSet edge2CloseEdge; const double eps = data._epsilon * data._epsilon; - for ( unsigned i = 0; i < data._edges.size(); ++i ) + for ( size_t i = 0; i < data._edges.size(); ++i ) { _LayerEdge* edge = data._edges[i]; if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue; @@ -2710,19 +3113,19 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, if ( FF1[0].IsNull() || FF2[0].IsNull() ) continue; -// // get a new normal for edge1 + // get a new normal for edge1 bool ok; gp_Vec dir1 = edge1->_normal, dir2 = edge2->_normal; if ( edge1->_cosin < 0 ) dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ).Normalized(); if ( edge2->_cosin < 0 ) dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok ).Normalized(); - // gp_Vec dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ); -// gp_Vec dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok2 ); -// double wgt1 = ( edge1->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 ); -// double wgt2 = ( edge2->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 ); -// gp_Vec newNorm = wgt1 * dir1 + wgt2 * dir2; -// newNorm.Normalize(); + // gp_Vec dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ); + // gp_Vec dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok2 ); + // double wgt1 = ( edge1->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 ); + // double wgt2 = ( edge2->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 ); + // gp_Vec newNorm = wgt1 * dir1 + wgt2 * dir2; + // newNorm.Normalize(); double wgt1 = ( edge1->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 ); double wgt2 = ( edge2->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 ); @@ -2736,7 +3139,7 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, n1 = edge1->_2neibors->_edges[0]->_nodes[0]; n2 = edge1->_2neibors->_edges[1]->_nodes[0]; //if ( !findNeiborsOnEdge( edge1, n1, n2, data )) - //continue; + // continue; edge1->SetDataByNeighbors( n1, n2, helper ); gp_Vec dirInFace; if ( edge1->_cosin < 0 ) @@ -2810,7 +3213,7 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, // 2) Check absence of intersections // TODO? - for ( unsigned i = 0 ; i < tmpFaces.size(); ++i ) + for ( size_t i = 0 ; i < tmpFaces.size(); ++i ) delete tmpFaces[i]; return true; @@ -2836,7 +3239,7 @@ bool _LayerEdge::FindIntersection( SMESH_ElementSearcher& searcher, bool segmentIntersected = false; distance = Precision::Infinite(); int iFace = -1; // intersected face - for ( unsigned j = 0 ; j < suspectFaces.size() && !segmentIntersected; ++j ) + for ( size_t j = 0 ; j < suspectFaces.size() && !segmentIntersected; ++j ) { const SMDS_MeshElement* face = suspectFaces[j]; if ( face->GetNodeIndex( _nodes.back() ) >= 0 || @@ -3124,7 +3527,7 @@ bool _LayerEdge::Smooth(int& badNb) // compute new position for the last _pos gp_XYZ newPos (0,0,0); - for ( unsigned i = 0; i < _simplices.size(); ++i ) + for ( size_t i = 0; i < _simplices.size(); ++i ) newPos += SMESH_TNodeXYZ( _simplices[i]._nPrev ); newPos /= _simplices.size(); @@ -3147,11 +3550,11 @@ bool _LayerEdge::Smooth(int& badNb) // count quality metrics (orientation) of tetras around _tgtNode int nbOkBefore = 0; SMESH_TNodeXYZ tgtXYZ( _nodes.back() ); - for ( unsigned i = 0; i < _simplices.size(); ++i ) + for ( size_t i = 0; i < _simplices.size(); ++i ) nbOkBefore += _simplices[i].IsForward( _nodes[0], &tgtXYZ ); int nbOkAfter = 0; - for ( unsigned i = 0; i < _simplices.size(); ++i ) + for ( size_t i = 0; i < _simplices.size(); ++i ) nbOkAfter += _simplices[i].IsForward( _nodes[0], &newPos ); if ( nbOkAfter < nbOkBefore ) @@ -3272,19 +3675,23 @@ bool _ViscousBuilder::refine(_SolidData& data) Handle(Geom_Surface) surface; TopoDS_Edge geomEdge; TopoDS_Face geomFace; + TopoDS_Shape prevSWOL; TopLoc_Location loc; - double f,l, u/*, distXYZ[4]*/; + double f,l, u; gp_XY uv; bool isOnEdge; + TGeomID prevBaseId = -1; + TNode2Edge* n2eMap = 0; + TNode2Edge::iterator n2e; - for ( unsigned i = 0; i < data._edges.size(); ++i ) + for ( size_t i = 0; i < data._edges.size(); ++i ) { _LayerEdge& edge = *data._edges[i]; // get accumulated length of segments vector< double > segLen( edge._pos.size() ); segLen[0] = 0.0; - for ( unsigned j = 1; j < edge._pos.size(); ++j ) + for ( size_t j = 1; j < edge._pos.size(); ++j ) segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus(); // allocate memory for new nodes if it is not yet refined @@ -3295,27 +3702,46 @@ bool _ViscousBuilder::refine(_SolidData& data) edge._nodes[1] = 0; edge._nodes.back() = tgtNode; } - if ( !edge._sWOL.IsNull() ) + // get data of a shrink shape + if ( !edge._sWOL.IsNull() && edge._sWOL != prevSWOL ) { isOnEdge = ( edge._sWOL.ShapeType() == TopAbs_EDGE ); - // restore position of the last node -// gp_Pnt p; if ( isOnEdge ) { geomEdge = TopoDS::Edge( edge._sWOL ); - curve = BRep_Tool::Curve( geomEdge, loc, f,l); -// double u = helper.GetNodeU( tgtNode ); -// p = curve->Value( u ); + curve = BRep_Tool::Curve( geomEdge, loc, f,l); } else { geomFace = TopoDS::Face( edge._sWOL ); - surface = BRep_Tool::Surface( geomFace, loc ); -// gp_XY uv = helper.GetNodeUV( tgtNode ); -// p = surface->Value( uv.X(), uv.Y() ); + surface = BRep_Tool::Surface( geomFace, loc ); + } + prevSWOL = edge._sWOL; + } + // restore shapePos of the last node by already treated _LayerEdge of another _SolidData + const TGeomID baseShapeId = edge._nodes[0]->getshapeId(); + if ( baseShapeId != prevBaseId ) + { + map< TGeomID, TNode2Edge* >::iterator s2ne = data._s2neMap.find( baseShapeId ); + n2eMap = ( s2ne == data._s2neMap.end() ) ? 0 : n2eMap = s2ne->second; + prevBaseId = baseShapeId; + } + if ( n2eMap && (( n2e = n2eMap->find( edge._nodes[0] )) != n2eMap->end() )) + { + _LayerEdge* foundEdge = n2e->second; + const gp_XYZ& foundPos = foundEdge->_pos.back(); + SMDS_PositionPtr lastPos = tgtNode->GetPosition(); + if ( isOnEdge ) + { + SMDS_EdgePosition* epos = static_cast( lastPos ); + epos->SetUParameter( foundPos.X() ); + } + else + { + SMDS_FacePosition* fpos = static_cast( lastPos ); + fpos->SetUParameter( foundPos.X() ); + fpos->SetVParameter( foundPos.Y() ); } -// p.Transform( loc ); -// const_cast< SMDS_MeshNode* >( tgtNode )->setXYZ( p.X(), p.Y(), p.Z() ); } // calculate height of the first layer double h0; @@ -3332,8 +3758,8 @@ bool _ViscousBuilder::refine(_SolidData& data) // create intermediate nodes double hSum = 0, hi = h0/f; - unsigned iSeg = 1; - for ( unsigned iStep = 1; iStep < edge._nodes.size(); ++iStep ) + size_t iSeg = 1; + for ( size_t iStep = 1; iStep < edge._nodes.size(); ++iStep ) { // compute an intermediate position hi *= f; @@ -3353,12 +3779,14 @@ bool _ViscousBuilder::refine(_SolidData& data) if ( isOnEdge ) { u = pos.X(); - pos = curve->Value( u ).Transformed(loc); + if ( !node ) + pos = curve->Value( u ).Transformed(loc); } else { uv.SetCoord( pos.X(), pos.Y() ); - pos = surface->Value( pos.X(), pos.Y() ).Transformed(loc); + if ( !node ) + pos = surface->Value( pos.X(), pos.Y() ).Transformed(loc); } } // create or update the node @@ -3386,11 +3814,18 @@ bool _ViscousBuilder::refine(_SolidData& data) { u = 0.5 * ( u + helper.GetNodeU( geomEdge, node )); pos = curve->Value( u ).Transformed(loc); + + SMDS_EdgePosition* epos = static_cast( node->GetPosition() ); + epos->SetUParameter( u ); } else { uv = 0.5 * ( uv + helper.GetNodeUV( geomFace, node )); pos = surface->Value( uv.X(), uv.Y()).Transformed(loc); + + SMDS_FacePosition* fpos = static_cast( node->GetPosition() ); + fpos->SetUParameter( uv.X() ); + fpos->SetVParameter( uv.Y() ); } } node->setXYZ( pos.X(), pos.Y(), pos.Z() ); @@ -3400,7 +3835,7 @@ bool _ViscousBuilder::refine(_SolidData& data) if ( !getMeshDS()->IsEmbeddedMode() ) // Log node movement - for ( unsigned i = 0; i < data._edges.size(); ++i ) + for ( size_t i = 0; i < data._edges.size(); ++i ) { _LayerEdge& edge = *data._edges[i]; SMESH_TNodeXYZ p ( edge._nodes.back() ); @@ -3414,7 +3849,7 @@ bool _ViscousBuilder::refine(_SolidData& data) TopExp_Explorer exp( data._solid, TopAbs_FACE ); for ( ; exp.More(); exp.Next() ) { - if ( _ignoreShapeIds.count( getMeshDS()->ShapeToIndex( exp.Current() ))) + if ( data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() ))) continue; SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() ); SMDS_ElemIteratorPtr fIt = fSubM->GetElements(); @@ -3465,7 +3900,7 @@ bool _ViscousBuilder::shrink() // make map of (ids of FACEs to shrink mesh on) to (_SolidData containing _LayerEdge's // inflated along FACE or EDGE) map< TGeomID, _SolidData* > f2sdMap; - for ( unsigned i = 0 ; i < _sdVec.size(); ++i ) + for ( size_t i = 0 ; i < _sdVec.size(); ++i ) { _SolidData& data = _sdVec[i]; TopTools_MapOfShape FFMap; @@ -3550,7 +3985,7 @@ bool _ViscousBuilder::shrink() sm->getDependsOnIterator(/*includeSelf=*/false, /*complexShapeFirst=*/false); while ( subIt->more() ) { - SMESH_subMesh* sub = subIt->next(); + SMESH_subMesh* sub = subIt->next(); SMESHDS_SubMesh* subDS = sub->GetSubMeshDS(); if ( subDS->NbNodes() == 0 || !n2eMap.count( subDS->GetNodes()->next() )) continue; @@ -3564,9 +3999,15 @@ bool _ViscousBuilder::shrink() } } + dumpFunction(SMESH_Comment("beforeShrinkFace")<first); // debug + SMDS_ElemIteratorPtr fIt = smDS->GetElements(); + while ( fIt->more() ) + if ( const SMDS_MeshElement* f = fIt->next() ) + dumpChangeNodes( f ); + // Replace source nodes by target nodes in mesh faces to shrink const SMDS_MeshNode* nodes[20]; - for ( unsigned i = 0; i < lEdges.size(); ++i ) + for ( size_t i = 0; i < lEdges.size(); ++i ) { _LayerEdge& edge = *lEdges[i]; const SMDS_MeshNode* srcNode = edge._nodes[0]; @@ -3577,10 +4018,10 @@ bool _ViscousBuilder::shrink() const SMDS_MeshElement* f = fIt->next(); if ( !smDS->Contains( f )) continue; - SMDS_ElemIteratorPtr nIt = f->nodesIterator(); - for ( int iN = 0; iN < f->NbNodes(); ++iN ) + SMDS_NodeIteratorPtr nIt = f->nodeIterator(); + for ( int iN = 0; nIt->more(); ++iN ) { - const SMDS_MeshNode* n = static_cast( nIt->next() ); + const SMDS_MeshNode* n = nIt->next(); nodes[iN] = ( n == srcNode ? tgtNode : n ); } helper.GetMeshDS()->ChangeElementNodes( f, nodes, f->NbNodes() ); @@ -3593,25 +4034,24 @@ bool _ViscousBuilder::shrink() // Create _SmoothNode's on face F vector< _SmoothNode > nodesToSmooth( smoothNodes.size() ); { - dumpFunction(SMESH_Comment("beforeShrinkFace")<first); // debug - for ( unsigned i = 0; i < smoothNodes.size(); ++i ) + const bool sortSimplices = isConcaveFace; + for ( size_t i = 0; i < smoothNodes.size(); ++i ) { const SMDS_MeshNode* n = smoothNodes[i]; nodesToSmooth[ i ]._node = n; // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices - getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, NULL, isConcaveFace ); + getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, NULL, sortSimplices ); // fix up incorrect uv of nodes on the FACE helper.GetNodeUV( F, n, 0, &isOkUV); dumpMove( n ); } - dumpFunctionEnd(); } //if ( nodesToSmooth.empty() ) continue; - // Find EDGE's to shrink + // Find EDGE's to shrink and set simpices to LayerEdge's set< _Shrinker1D* > eShri1D; { - for ( unsigned i = 0; i < lEdges.size(); ++i ) + for ( size_t i = 0; i < lEdges.size(); ++i ) { _LayerEdge* edge = lEdges[i]; if ( edge->_sWOL.ShapeType() == TopAbs_EDGE ) @@ -3625,6 +4065,23 @@ bool _ViscousBuilder::shrink() // srinked while srinking another FACE srinker.RestoreParams(); } + getSimplices( /*tgtNode=*/edge->_nodes.back(), edge->_simplices, ignoreShapes ); + } + } + + bool toFixTria = false; // to improve quality of trias by diagonal swap + if ( isConcaveFace ) + { + const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles(); + if ( hasTria != hasQuad ) { + toFixTria = hasTria; + } + else { + set nbNodesSet; + SMDS_ElemIteratorPtr fIt = smDS->GetElements(); + while ( fIt->more() && nbNodesSet.size() < 2 ) + nbNodesSet.insert( fIt->next()->NbCornerNodes() ); + toFixTria = ( *nbNodesSet.begin() == 3 ); } } @@ -3634,19 +4091,23 @@ bool _ViscousBuilder::shrink() bool shrinked = true; int badNb, shriStep=0, smooStep=0; + _SmoothNode::SmoothType smoothType + = isConcaveFace ? _SmoothNode::ANGULAR : _SmoothNode::LAPLACIAN; while ( shrinked ) { + shriStep++; // Move boundary nodes (actually just set new UV) // ----------------------------------------------- - dumpFunction(SMESH_Comment("moveBoundaryOnF")<first<<"_st"<first<<"_st"<SetNewLength2d( surface,F,helper ); } dumpFunctionEnd(); // Move nodes on EDGE's + // (XYZ is set as soon as a needed length reached in SetNewLength2d()) set< _Shrinker1D* >::iterator shr = eShri1D.begin(); for ( ; shr != eShri1D.end(); ++shr ) (*shr)->Compute( /*set3D=*/false, helper ); @@ -3654,7 +4115,7 @@ bool _ViscousBuilder::shrink() // Smoothing in 2D // ----------------- int nbNoImpSteps = 0; - bool moved = true; + bool moved = true; badNb = 1; while (( nbNoImpSteps < 5 && badNb > 0) && moved) { @@ -3663,11 +4124,10 @@ bool _ViscousBuilder::shrink() int oldBadNb = badNb; badNb = 0; moved = false; - for ( unsigned i = 0; i < nodesToSmooth.size(); ++i ) + for ( size_t i = 0; i < nodesToSmooth.size(); ++i ) { moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign, - /*isCentroidal=*/isConcaveFace, - /*set3D=*/isConcaveFace); + smoothType, /*set3D=*/isConcaveFace); } if ( badNb < oldBadNb ) nbNoImpSteps = 0; @@ -3678,39 +4138,75 @@ bool _ViscousBuilder::shrink() } if ( badNb > 0 ) return error(SMESH_Comment("Can't shrink 2D mesh on face ") << f2sd->first ); - } + if ( shriStep > 200 ) + return error(SMESH_Comment("Infinite loop at shrinking 2D mesh on face ") << f2sd->first ); + + // Fix narrow triangles by swapping diagonals + // --------------------------------------- + if ( toFixTria ) + { + set usedNodes; + fixBadFaces( F, helper, /*is2D=*/true, shriStep, & usedNodes); // swap diagonals + + // update working data + set::iterator n; + for ( size_t i = 0; i < nodesToSmooth.size() && !usedNodes.empty(); ++i ) + { + n = usedNodes.find( nodesToSmooth[ i ]._node ); + if ( n != usedNodes.end()) + { + getSimplices( nodesToSmooth[ i ]._node, + nodesToSmooth[ i ]._simplices, + ignoreShapes, NULL, + /*sortSimplices=*/ smoothType == _SmoothNode::ANGULAR ); + usedNodes.erase( n ); + } + } + for ( size_t i = 0; i < lEdges.size() && !usedNodes.empty(); ++i ) + { + n = usedNodes.find( /*tgtNode=*/ lEdges[i]->_nodes.back() ); + if ( n != usedNodes.end()) + { + getSimplices( lEdges[i]->_nodes.back(), + lEdges[i]->_simplices, + ignoreShapes ); + usedNodes.erase( n ); + } + } + } + } // while ( shrinked ) + // No wrongly shaped faces remain; final smooth. Set node XYZ. - // First, find out a needed quality of smoothing (high for quadrangles only) - bool highQuality; + bool isStructuredFixed = false; + if ( SMESH_2D_Algo* algo = dynamic_cast( sm->GetAlgo() )) + isStructuredFixed = algo->FixInternalNodes( *data._proxyMesh, F ); + if ( !isStructuredFixed ) { - const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles(); - if ( hasTria != hasQuad ) { - highQuality = hasQuad; - } - else { - set nbNodesSet; - SMDS_ElemIteratorPtr fIt = smDS->GetElements(); - while ( fIt->more() && nbNodesSet.size() < 2 ) - nbNodesSet.insert( fIt->next()->NbCornerNodes() ); - highQuality = ( *nbNodesSet.begin() == 4 ); + if ( isConcaveFace ) // fix narrow faces by swapping diagonals + fixBadFaces( F, helper, /*is2D=*/false, ++shriStep ); + + for ( int st = 3; st; --st ) + { + switch( st ) { + case 1: smoothType = _SmoothNode::LAPLACIAN; break; + case 2: smoothType = _SmoothNode::LAPLACIAN; break; + case 3: smoothType = _SmoothNode::ANGULAR; break; + } + dumpFunction(SMESH_Comment("shrinkFace")<first<<"_st"<<++smooStep); // debug + for ( size_t i = 0; i < nodesToSmooth.size(); ++i ) + { + nodesToSmooth[i].Smooth( badNb,surface,helper,refSign, + smoothType,/*set3D=*/st==1 ); + } + dumpFunctionEnd(); } } - if ( !highQuality && isConcaveFace ) - fixBadFaces( F, helper ); // fix narrow faces by swaping diagonals - for ( int st = highQuality ? 10 : 3; st; --st ) - { - dumpFunction(SMESH_Comment("shrinkFace")<first<<"_st"<<++smooStep); // debug - for ( unsigned i = 0; i < nodesToSmooth.size(); ++i ) - nodesToSmooth[i].Smooth( badNb,surface,helper,refSign, - /*isCentroidal=*/isConcaveFace,/*set3D=*/st==1 ); - dumpFunctionEnd(); - } // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh VISCOUS_3D::ToClearSubWithMain( sm, data._solid ); if ( !getMeshDS()->IsEmbeddedMode() ) // Log node movement - for ( unsigned i = 0; i < nodesToSmooth.size(); ++i ) + for ( size_t i = 0; i < nodesToSmooth.size(); ++i ) { SMESH_TNodeXYZ p ( nodesToSmooth[i]._node ); getMeshDS()->MoveNode( nodesToSmooth[i]._node, p.X(), p.Y(), p.Z() ); @@ -3752,55 +4248,56 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, double uvLen = uvDir.Magnitude(); uvDir /= uvLen; edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0); - - // IMPORTANT to have src nodes NOT yet REPLACED by tgt nodes in shrinked faces - vector faces; - multimap< double, const SMDS_MeshNode* > proj2node; - SMDS_ElemIteratorPtr fIt = srcNode->GetInverseElementIterator(SMDSAbs_Face); - while ( fIt->more() ) - { - const SMDS_MeshElement* f = fIt->next(); - if ( faceSubMesh->Contains( f )) - faces.push_back( f ); - } - for ( unsigned i = 0; i < faces.size(); ++i ) - { - const int nbNodes = faces[i]->NbCornerNodes(); - for ( int j = 0; j < nbNodes; ++j ) - { - const SMDS_MeshNode* n = faces[i]->GetNode(j); - if ( n == srcNode ) continue; - if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE && - ( faces.size() > 1 || nbNodes > 3 )) - continue; - gp_Pnt2d uv = helper.GetNodeUV( F, n ); - gp_Vec2d uvDirN( srcUV, uv ); - double proj = uvDirN * uvDir; - proj2node.insert( make_pair( proj, n )); - } - } - - multimap< double, const SMDS_MeshNode* >::iterator p2n = proj2node.begin(), p2nEnd; - const double minProj = p2n->first; - const double projThreshold = 1.1 * uvLen; - if ( minProj > projThreshold ) - { - // tgtNode is located so that it does not make faces with wrong orientation - return true; - } + edge._len = uvLen; + + // // IMPORTANT to have src nodes NOT yet REPLACED by tgt nodes in shrinked faces + // vector faces; + // multimap< double, const SMDS_MeshNode* > proj2node; + // SMDS_ElemIteratorPtr fIt = srcNode->GetInverseElementIterator(SMDSAbs_Face); + // while ( fIt->more() ) + // { + // const SMDS_MeshElement* f = fIt->next(); + // if ( faceSubMesh->Contains( f )) + // faces.push_back( f ); + // } + // for ( size_t i = 0; i < faces.size(); ++i ) + // { + // const int nbNodes = faces[i]->NbCornerNodes(); + // for ( int j = 0; j < nbNodes; ++j ) + // { + // const SMDS_MeshNode* n = faces[i]->GetNode(j); + // if ( n == srcNode ) continue; + // if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE && + // ( faces.size() > 1 || nbNodes > 3 )) + // continue; + // gp_Pnt2d uv = helper.GetNodeUV( F, n ); + // gp_Vec2d uvDirN( srcUV, uv ); + // double proj = uvDirN * uvDir; + // proj2node.insert( make_pair( proj, n )); + // } + // } + + // multimap< double, const SMDS_MeshNode* >::iterator p2n = proj2node.begin(), p2nEnd; + // const double minProj = p2n->first; + // const double projThreshold = 1.1 * uvLen; + // if ( minProj > projThreshold ) + // { + // // tgtNode is located so that it does not make faces with wrong orientation + // return true; + // } edge._pos.resize(1); edge._pos[0].SetCoord( tgtUV.X(), tgtUV.Y(), 0 ); // store most risky nodes in _simplices - p2nEnd = proj2node.lower_bound( projThreshold ); - int nbSimpl = ( std::distance( p2n, p2nEnd ) + 1) / 2; - edge._simplices.resize( nbSimpl ); - for ( int i = 0; i < nbSimpl; ++i ) - { - edge._simplices[i]._nPrev = p2n->second; - if ( ++p2n != p2nEnd ) - edge._simplices[i]._nNext = p2n->second; - } + // p2nEnd = proj2node.lower_bound( projThreshold ); + // int nbSimpl = ( std::distance( p2n, p2nEnd ) + 1) / 2; + // edge._simplices.resize( nbSimpl ); + // for ( int i = 0; i < nbSimpl; ++i ) + // { + // edge._simplices[i]._nPrev = p2n->second; + // if ( ++p2n != p2nEnd ) + // edge._simplices[i]._nNext = p2n->second; + // } // set UV of source node to target node SMDS_FacePosition* pos = static_cast( tgtNode->GetPosition() ); pos->SetUParameter( srcUV.X() ); @@ -3908,23 +4405,28 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, */ //================================================================================ -void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& helper) +void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, + SMESH_MesherHelper& helper, + const bool is2D, + const int step, + set * involvedNodes) { SMESH::Controls::AspectRatio qualifier; SMESH::Controls::TSequenceOfXYZ points(3), points1(3), points2(3); - const double maxAspectRatio = 4.; + const double maxAspectRatio = is2D ? 4. : 2; + NodeCoordHelper xyz( F, helper, is2D ); // find bad triangles vector< const SMDS_MeshElement* > badTrias; vector< double > badAspects; - SMESHDS_SubMesh* sm = helper.GetMeshDS()->MeshElements( F ); + SMESHDS_SubMesh* sm = helper.GetMeshDS()->MeshElements( F ); SMDS_ElemIteratorPtr fIt = sm->GetElements(); while ( fIt->more() ) { const SMDS_MeshElement * f = fIt->next(); if ( f->NbCornerNodes() != 3 ) continue; - for ( int iP = 0; iP < 3; ++iP ) points(iP+1) = SMESH_TNodeXYZ( f->GetNode(iP)); + for ( int iP = 0; iP < 3; ++iP ) points(iP+1) = xyz( f->GetNode(iP)); double aspect = qualifier.GetValue( points ); if ( aspect > maxAspectRatio ) { @@ -3932,6 +4434,18 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help badAspects.push_back( aspect ); } } + if ( step == 1 ) + { + dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<GetElements(); + while ( fIt->more() ) + { + const SMDS_MeshElement * f = fIt->next(); + if ( f->NbCornerNodes() == 3 ) + dumpChangeNodes( f ); + } + dumpFunctionEnd(); + } if ( badTrias.empty() ) return; @@ -3947,9 +4461,10 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help double aspRatio [3]; int i1, i2, i3; - involvedFaces.insert( badTrias[iTia] ); + if ( !involvedFaces.insert( badTrias[iTia] ).second ) + continue; for ( int iP = 0; iP < 3; ++iP ) - points(iP+1) = SMESH_TNodeXYZ( badTrias[iTia]->GetNode(iP)); + points(iP+1) = xyz( badTrias[iTia]->GetNode(iP)); // find triangles adjacent to badTrias[iTia] with better aspect ratio after diag-swaping int bestCouple = -1; @@ -3958,14 +4473,16 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help const SMDS_MeshNode* n1 = badTrias[iTia]->GetNode( iSide ); const SMDS_MeshNode* n2 = badTrias[iTia]->GetNode(( iSide+1 ) % 3 ); trias [iSide].first = badTrias[iTia]; - trias [iSide].second = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, involvedFaces, - & i1, & i2 ); - if ( ! trias[iSide].second || trias[iSide].second->NbCornerNodes() != 3 ) + trias [iSide].second = SMESH_MeshAlgos::FindFaceInSet( n1, n2, emptySet, involvedFaces, + & i1, & i2 ); + if (( ! trias[iSide].second ) || + ( trias[iSide].second->NbCornerNodes() != 3 ) || + ( ! sm->Contains( trias[iSide].second ))) continue; // aspect ratio of an adjacent tria for ( int iP = 0; iP < 3; ++iP ) - points2(iP+1) = SMESH_TNodeXYZ( trias[iSide].second->GetNode(iP)); + points2(iP+1) = xyz( trias[iSide].second->GetNode(iP)); double aspectInit = qualifier.GetValue( points2 ); // arrange nodes as after diag-swaping @@ -3982,6 +4499,12 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help if ( aspRatio[ iSide ] > aspectInit + badAspects[ iTia ] ) continue; + // prevent inversion of a triangle + gp_Vec norm1 = gp_Vec( points1(1), points1(3) ) ^ gp_Vec( points1(1), points1(2) ); + gp_Vec norm2 = gp_Vec( points2(1), points2(3) ) ^ gp_Vec( points2(1), points2(2) ); + if ( norm1 * norm2 < 0. && norm1.Angle( norm2 ) > 70./180.*M_PI ) + continue; + if ( bestCouple < 0 || aspRatio[ bestCouple ] > aspRatio[ iSide ] ) bestCouple = iSide; } @@ -4002,17 +4525,25 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help // swap diagonals SMESH_MeshEditor editor( helper.GetMesh() ); - dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<insert( triaCouples[i].first->begin_nodes(), + triaCouples[i].first->end_nodes() ); + involvedNodes->insert( triaCouples[i].second->begin_nodes(), + triaCouples[i].second->end_nodes() ); + } // just for debug dump resulting triangles - dumpFunction(SMESH_Comment("swapDiagonals_F")< minStepSize ) - stepSize = proj; - } + // find intersection of 2 lines: curUV-tgtUV and that connecting simplex nodes + gp_XY uvN1 = helper.GetNodeUV( F, _simplices[i]._nPrev ); + gp_XY uvN2 = helper.GetNodeUV( F, _simplices[i]._nNext ); + gp_XY dirN = uvN2 - uvN1; + double det = uvDir.Crossed( dirN ); + if ( Abs( det ) < std::numeric_limits::min() ) continue; + gp_XY dirN2Cur = curUV - uvN1; + double step = dirN.Crossed( dirN2Cur ) / det; + if ( step > 0 ) + stepSize = Min( step, stepSize ); } - gp_Pnt2d newUV; - if ( stepSize == uvLen ) + if ( uvLen - stepSize < _len / 200. ) { newUV = tgtUV; _pos.clear(); } + else if ( stepSize > 0 ) + { + newUV = curUV + uvDir.XY() * stepSize * kSafe; + } else { - newUV = curUV + uvDir.XY() * stepSize; + return true; } - SMDS_FacePosition* pos = static_cast( tgtNode->GetPosition() ); pos->SetUParameter( newUV.X() ); pos->SetVParameter( newUV.Y() ); @@ -4085,22 +4617,22 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface, { TopoDS_Edge E = TopoDS::Edge( _sWOL ); const SMDS_MeshNode* n2 = _simplices[0]._nPrev; + SMDS_EdgePosition* tgtPos = static_cast( tgtNode->GetPosition() ); const double u2 = helper.GetNodeU( E, n2, tgtNode ); const double uSrc = _pos[0].Coord( U_SRC ); const double lenTgt = _pos[0].Coord( LEN_TGT ); double newU = _pos[0].Coord( U_TGT ); - if ( lenTgt < 0.99 * fabs( uSrc-u2 )) + if ( lenTgt < 0.99 * fabs( uSrc-u2 )) // n2 got out of src-tgt range { _pos.clear(); } else { - newU = 0.1 * uSrc + 0.9 * u2; + newU = 0.1 * tgtPos->GetUParameter() + 0.9 * u2; } - SMDS_EdgePosition* pos = static_cast( tgtNode->GetPosition() ); - pos->SetUParameter( newU ); + tgtPos->SetUParameter( newU ); #ifdef __myDEBUG gp_XY newUV = helper.GetNodeUV( F, tgtNode, _nodes[0]); gp_Pnt p = surface->Value( newUV.X(), newUV.Y() ); @@ -4122,7 +4654,7 @@ bool _SmoothNode::Smooth(int& badNb, Handle(Geom_Surface)& surface, SMESH_MesherHelper& helper, const double refSign, - bool isCentroidal, + SmoothType how, bool set3D) { const TopoDS_Face& face = TopoDS::Face( helper.GetSubShape() ); @@ -4134,7 +4666,24 @@ bool _SmoothNode::Smooth(int& badNb, // compute new UV for the node gp_XY newPos (0,0); - if ( isCentroidal && _simplices.size() > 3 ) + if ( how == TFI && _simplices.size() == 4 ) + { + gp_XY corners[4]; + for ( size_t i = 0; i < _simplices.size(); ++i ) + if ( _simplices[i]._nOpp ) + corners[i] = helper.GetNodeUV( face, _simplices[i]._nOpp, _node ); + else + throw SALOME_Exception(LOCALIZED("TFI smoothing: _Simplex::_nOpp not set!")); + + newPos = helper.calcTFI ( 0.5, 0.5, + corners[0], corners[1], corners[2], corners[3], + uv[1], uv[2], uv[3], uv[0] ); + } + else if ( how == ANGULAR ) + { + newPos = computeAngularPos( uv, helper.GetNodeUV( face, _node ), refSign ); + } + else if ( how == CENTROIDAL && _simplices.size() > 3 ) { // average centers of diagonals wieghted with their reciprocal lengths if ( _simplices.size() == 4 ) @@ -4159,13 +4708,12 @@ bool _SmoothNode::Smooth(int& badNb, newPos += w * ( uv[i]+uv[i2] ); } } - newPos /= 2 * sumWeight; + newPos /= 2 * sumWeight; // 2 is to get a middle between uv's } } else { // Laplacian smooth - isCentroidal = false; for ( size_t i = 0; i < _simplices.size(); ++i ) newPos += uv[i]; newPos /= _simplices.size(); @@ -4174,17 +4722,15 @@ bool _SmoothNode::Smooth(int& badNb, // count quality metrics (orientation) of triangles around the node int nbOkBefore = 0; gp_XY tgtUV = helper.GetNodeUV( face, _node ); - for ( unsigned i = 0; i < _simplices.size(); ++i ) + for ( size_t i = 0; i < _simplices.size(); ++i ) nbOkBefore += _simplices[i].IsForward( tgtUV, _node, face, helper, refSign ); int nbOkAfter = 0; - for ( unsigned i = 0; i < _simplices.size(); ++i ) + for ( size_t i = 0; i < _simplices.size(); ++i ) nbOkAfter += _simplices[i].IsForward( newPos, _node, face, helper, refSign ); if ( nbOkAfter < nbOkBefore ) { - // if ( isCentroidal ) - // return Smooth( badNb, surface, helper, refSign, !isCentroidal, set3D ); badNb += _simplices.size() - nbOkBefore; return false; } @@ -4207,6 +4753,66 @@ bool _SmoothNode::Smooth(int& badNb, return ( (tgtUV-newPos).SquareModulus() > 1e-10 ); } +//================================================================================ +/*! + * \brief Computes new UV using angle based smoothing technic + */ +//================================================================================ + +gp_XY _SmoothNode::computeAngularPos(vector& uv, + const gp_XY& uvToFix, + const double refSign) +{ + uv.push_back( uv.front() ); + + vector< gp_XY > edgeDir ( uv.size() ); + vector< double > edgeSize( uv.size() ); + for ( size_t i = 1; i < edgeDir.size(); ++i ) + { + edgeDir [i-1] = uv[i] - uv[i-1]; + edgeSize[i-1] = edgeDir[i-1].Modulus(); + if ( edgeSize[i-1] < numeric_limits::min() ) + edgeDir[i-1].SetX( 100 ); + else + edgeDir[i-1] /= edgeSize[i-1] * refSign; + } + edgeDir.back() = edgeDir.front(); + edgeSize.back() = edgeSize.front(); + + gp_XY newPos(0,0); + int nbEdges = 0; + double sumSize = 0; + for ( size_t i = 1; i < edgeDir.size(); ++i ) + { + if ( edgeDir[i-1].X() > 1. ) continue; + int i1 = i-1; + while ( edgeDir[i].X() > 1. && ++i < edgeDir.size() ); + if ( i == edgeDir.size() ) break; + gp_XY p = uv[i]; + gp_XY norm1( -edgeDir[i1].Y(), edgeDir[i1].X() ); + gp_XY norm2( -edgeDir[i].Y(), edgeDir[i].X() ); + gp_XY bisec = norm1 + norm2; + double bisecSize = bisec.Modulus(); + if ( bisecSize < numeric_limits::min() ) + { + bisec = -edgeDir[i1] + edgeDir[i]; + bisecSize = bisec.Modulus(); + } + bisec /= bisecSize; + + gp_XY dirToN = uvToFix - p; + double distToN = dirToN.Modulus(); + if ( bisec * dirToN < 0 ) + distToN = -distToN; + + newPos += ( p + bisec * distToN ) * ( edgeSize[i1] + edgeSize[i] ); + ++nbEdges; + sumSize += edgeSize[i1] + edgeSize[i]; + } + newPos /= /*nbEdges * */sumSize; + return newPos; +} + //================================================================================ /*! * \brief Delete _SolidData @@ -4215,7 +4821,7 @@ bool _SmoothNode::Smooth(int& badNb, _SolidData::~_SolidData() { - for ( unsigned i = 0; i < _edges.size(); ++i ) + for ( size_t i = 0; i < _edges.size(); ++i ) { if ( _edges[i] && _edges[i]->_2neibors ) delete _edges[i]->_2neibors; @@ -4288,7 +4894,7 @@ void _Shrinker1D::AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper ) { // remove target node of the _LayerEdge from _nodes int nbFound = 0; - for ( unsigned i = 0; i < _nodes.size(); ++i ) + for ( size_t i = 0; i < _nodes.size(); ++i ) if ( !_nodes[i] || _nodes[i] == tgtNode0 || _nodes[i] == tgtNode1 ) _nodes[i] = 0, nbFound++; if ( nbFound == _nodes.size() ) @@ -4326,7 +4932,7 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper) l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() ); double totLen = GCPnts_AbscissaPoint::Length( aCurve, f, l ); - for ( unsigned i = 0; i < _nodes.size(); ++i ) + for ( size_t i = 0; i < _nodes.size(); ++i ) { if ( !_nodes[i] ) continue; double len = totLen * _normPar[i]; @@ -4348,7 +4954,7 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper) if ( _edges[1] ) l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() ); - for ( unsigned i = 0; i < _nodes.size(); ++i ) + for ( size_t i = 0; i < _nodes.size(); ++i ) { if ( !_nodes[i] ) continue; double u = f * ( 1-_normPar[i] ) + l * _normPar[i]; @@ -4367,7 +4973,7 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper) void _Shrinker1D::RestoreParams() { if ( _done ) - for ( unsigned i = 0; i < _nodes.size(); ++i ) + for ( size_t i = 0; i < _nodes.size(); ++i ) { if ( !_nodes[i] ) continue; SMDS_EdgePosition* pos = static_cast( _nodes[i]->GetPosition() ); @@ -4420,7 +5026,7 @@ bool _ViscousBuilder::addBoundaryElements() { SMESH_MesherHelper helper( *_mesh ); - for ( unsigned i = 0; i < _sdVec.size(); ++i ) + for ( size_t i = 0; i < _sdVec.size(); ++i ) { _SolidData& data = _sdVec[i]; TopTools_IndexedMapOfShape geomEdges; @@ -4480,7 +5086,7 @@ bool _ViscousBuilder::addBoundaryElements() reverse = ( helper.GetSubShapeOri( F, E ) == TopAbs_REVERSED ); if ( helper.GetSubShapeOri( data._solid, F ) == TopAbs_REVERSED ) reverse = !reverse, F.Reverse(); - if ( SMESH_Algo::IsReversedSubMesh( TopoDS::Face(F), getMeshDS() )) + if ( helper.IsReversedSubMesh( TopoDS::Face(F) )) reverse = !reverse; } else @@ -4491,7 +5097,7 @@ bool _ViscousBuilder::addBoundaryElements() { const TopoDS_Shape* pF = fIt->next(); if ( helper.IsSubShape( *pF, data._solid) && - !_ignoreShapeIds.count( e2f->first )) + !data._ignoreFaceIds.count( e2f->first )) F = *pF; } } @@ -4507,7 +5113,7 @@ bool _ViscousBuilder::addBoundaryElements() // Make faces const int dj1 = reverse ? 0 : 1; const int dj2 = reverse ? 1 : 0; - for ( unsigned j = 1; j < ledges.size(); ++j ) + for ( size_t j = 1; j < ledges.size(); ++j ) { vector< const SMDS_MeshNode*>& nn1 = ledges[j-dj1]->_nodes; vector< const SMDS_MeshNode*>& nn2 = ledges[j-dj2]->_nodes;