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=67c79fef737d2928a6b25a49458d7c07dbdbf77e;hb=050aa87698efb9b123985f0bd69aab20ad42c728;hpb=a27b03972c8d020b1c78b4fd85a6e0430bac702b diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index 67c79fef7..b2f6855b0 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -1,9 +1,9 @@ -// Copyright (C) 2007-2013 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,14 +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 "StdMeshers_FaceSide.hxx" #include +#include +#include #include #include #include @@ -52,16 +56,19 @@ #include #include #include +#include #include #include #include #include #include #include +#include #include #include #include #include +#include #include #include #include @@ -319,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 @@ -358,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 @@ -381,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 @@ -413,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, @@ -433,7 +448,7 @@ namespace VISCOUS_3D //vector _nodesAround; vector<_Simplex> _simplices; // for quality check - enum SmoothType { LAPLACIAN, CENTROIDAL, ANGULAR }; + enum SmoothType { LAPLACIAN, CENTROIDAL, ANGULAR, TFI }; bool Smooth(int& badNb, Handle(Geom_Surface)& surface, @@ -471,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, @@ -479,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, @@ -499,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 ); @@ -512,7 +536,6 @@ namespace VISCOUS_3D SMESH_ComputeErrorPtr _error; vector< _SolidData > _sdVec; - set _ignoreShapeIds; int _tmpFaceID; }; //-------------------------------------------------------------------------------- @@ -548,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()));} }; //-------------------------------------------------------------------------------- @@ -567,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 //================================================================================ @@ -574,7 +635,7 @@ virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const // StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH_Gen* gen) :SMESH_Hypothesis(hypId, studyId, gen), - _isToIgnoreShapes(18), _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 @@ -643,7 +704,7 @@ std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save) << " " << _thickness << " " << _stretchFactor << " " << _shapeIds.size(); - for ( unsigned i = 0; i < _shapeIds.size(); ++i ) + for ( size_t i = 0; i < _shapeIds.size(); ++i ) save << " " << _shapeIds[i]; save << " " << !_isToIgnoreShapes; // negate to keep the behavior in old studies. return save; @@ -745,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(); @@ -780,6 +847,8 @@ 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 ); @@ -801,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; @@ -831,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; } //-------------------------------------------------------------------------------- @@ -842,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(); @@ -1162,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 )) @@ -1195,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); @@ -1212,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 ); } } @@ -1226,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; } @@ -1284,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() ) @@ -1295,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; } } @@ -1360,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; } @@ -1401,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 ) @@ -1421,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(); @@ -1475,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 ) { @@ -1489,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 @@ -1506,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() ) { @@ -1550,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 @@ -1558,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 ); @@ -1597,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() ); @@ -1683,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 ) { @@ -1774,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; @@ -1810,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 @@ -1913,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; @@ -1925,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; } //================================================================================ @@ -1946,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; } //================================================================================ @@ -1961,14 +2326,15 @@ 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 )); @@ -2006,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; @@ -2016,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])); @@ -2027,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] ); @@ -2081,7 +2447,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) auto_ptr searcher ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), data._proxyMesh->GetFaces( data._solid )) ); - for ( unsigned i = 0; i < data._edges.size(); ++i ) + for ( size_t i = 0; i < data._edges.size(); ++i ) { if ( data._edges[i]->IsOnEdge() ) continue; data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon ); @@ -2116,7 +2482,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) // Elongate _LayerEdge's dumpFunction(SMESH_Comment("inflate")<SetNewLength( curThick, helper ); } @@ -2132,7 +2498,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) if ( nbSteps > 0 ) { dumpFunction(SMESH_Comment("invalidate")<InvalidateStep( nbSteps+1 ); } @@ -2144,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 @@ -2192,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 ]; @@ -2237,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 ) @@ -2258,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() @@ -2273,6 +2643,24 @@ 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 @@ -2287,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; @@ -2598,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; @@ -2615,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 ) @@ -2654,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; @@ -2725,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 ); @@ -2751,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 ) @@ -2825,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; @@ -2851,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 || @@ -3139,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(); @@ -3162,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 ) @@ -3287,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 @@ -3310,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; @@ -3347,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; @@ -3368,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 @@ -3401,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() ); @@ -3415,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() ); @@ -3429,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(); @@ -3480,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; @@ -3579,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]; @@ -3592,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() ); @@ -3608,9 +4034,8 @@ bool _ViscousBuilder::shrink() // Create _SmoothNode's on face F vector< _SmoothNode > nodesToSmooth( smoothNodes.size() ); { - dumpFunction(SMESH_Comment("beforeShrinkFace")<first); // debug const bool sortSimplices = isConcaveFace; - for ( unsigned i = 0; i < smoothNodes.size(); ++i ) + for ( size_t i = 0; i < smoothNodes.size(); ++i ) { const SMDS_MeshNode* n = smoothNodes[i]; nodesToSmooth[ i ]._node = n; @@ -3620,14 +4045,13 @@ bool _ViscousBuilder::shrink() 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 ) @@ -3641,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 ); } } @@ -3651,14 +4092,15 @@ bool _ViscousBuilder::shrink() bool shrinked = true; int badNb, shriStep=0, smooStep=0; _SmoothNode::SmoothType smoothType - = isConcaveFace ? _SmoothNode::CENTROIDAL : _SmoothNode::LAPLACIAN; + = 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 ); } @@ -3673,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) { @@ -3682,7 +4124,7 @@ 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, smoothType, /*set3D=*/isConcaveFace); @@ -3696,7 +4138,43 @@ 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. bool isStructuredFixed = false; @@ -3704,12 +4182,18 @@ bool _ViscousBuilder::shrink() isStructuredFixed = algo->FixInternalNodes( *data._proxyMesh, F ); if ( !isStructuredFixed ) { - if ( isConcaveFace ) - fixBadFaces( F, helper ); // fix narrow faces by swapping diagonals - for ( int st = /*highQuality ? 10 :*/ 3; st; --st ) + 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 ( unsigned i = 0; i < nodesToSmooth.size(); ++i ) + for ( size_t i = 0; i < nodesToSmooth.size(); ++i ) { nodesToSmooth[i].Smooth( badNb,surface,helper,refSign, smoothType,/*set3D=*/st==1 ); @@ -3722,7 +4206,7 @@ bool _ViscousBuilder::shrink() 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() ); @@ -3766,54 +4250,54 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0); 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 ( 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 )); - } - } + // // 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; - } + // 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() ); @@ -3921,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 ) { @@ -3945,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; @@ -3960,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; @@ -3973,12 +4475,14 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help trias [iSide].first = badTrias[iTia]; trias [iSide].second = SMESH_MeshAlgos::FindFaceInSet( n1, n2, emptySet, involvedFaces, & i1, & i2 ); - if ( ! trias[iSide].second || trias[iSide].second->NbCornerNodes() != 3 ) + 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 @@ -3995,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; } @@ -4015,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; - else if ( proj < minStepSize ) - stepSize = minStepSize; - } + // 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 ( uvLen - stepSize < _len / 20. ) + 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() ); @@ -4149,30 +4666,24 @@ bool _SmoothNode::Smooth(int& badNb, // compute new UV for the node gp_XY newPos (0,0); -/* if ( how == ANGULAR && _simplices.size() == 4 ) + if ( how == TFI && _simplices.size() == 4 ) { - vector corners; corners.reserve(4); + gp_XY corners[4]; for ( size_t i = 0; i < _simplices.size(); ++i ) if ( _simplices[i]._nOpp ) - corners.push_back( helper.GetNodeUV( face, _simplices[i]._nOpp, _node )); - if ( corners.size() == 4 ) - { - newPos = helper.calcTFI - ( 0.5, 0.5, - corners[0], corners[1], corners[2], corners[3], - uv[1], uv[2], uv[3], uv[0] ); - } - // vector p( _simplices.size() * 2 + 1 ); - // p.clear(); - // for ( size_t i = 0; i < _simplices.size(); ++i ) - // { - // p.push_back( uv[i] ); - // if ( _simplices[i]._nOpp ) - // p.push_back( helper.GetNodeUV( face, _simplices[i]._nOpp, _node )); - // } - // newPos = computeAngularPos( p, helper.GetNodeUV( face, _node ), refSign ); + 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 == CENTROIDAL && _simplices.size() > 3 ) + 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 ) @@ -4203,7 +4714,6 @@ bool _SmoothNode::Smooth(int& badNb, else { // Laplacian smooth - //isCentroidal = false; for ( size_t i = 0; i < _simplices.size(); ++i ) newPos += uv[i]; newPos /= _simplices.size(); @@ -4212,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; } @@ -4257,22 +4765,22 @@ gp_XY _SmoothNode::computeAngularPos(vector& uv, { uv.push_back( uv.front() ); - vector< gp_XY > edgeDir( uv.size() ); + 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]; + 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(); + edgeDir.back() = edgeDir.front(); edgeSize.back() = edgeSize.front(); - gp_XY newPos(0,0); - int nbEdges = 0; + gp_XY newPos(0,0); + int nbEdges = 0; double sumSize = 0; for ( size_t i = 1; i < edgeDir.size(); ++i ) { @@ -4292,7 +4800,7 @@ gp_XY _SmoothNode::computeAngularPos(vector& uv, } bisec /= bisecSize; - gp_XY dirToN = uvToFix - p; + gp_XY dirToN = uvToFix - p; double distToN = dirToN.Modulus(); if ( bisec * dirToN < 0 ) distToN = -distToN; @@ -4313,7 +4821,7 @@ gp_XY _SmoothNode::computeAngularPos(vector& uv, _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; @@ -4386,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() ) @@ -4424,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]; @@ -4446,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]; @@ -4465,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() ); @@ -4518,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; @@ -4589,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; } } @@ -4605,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;