X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ViscousLayers.cxx;h=a2a44a23829f391426f09303ed72a60a62c9baef;hp=241a21636aa89f76dc89d610accf9120ac302e8c;hb=fa9a9581d37a45a53111f7df9ae5429c9692851f;hpb=2e439615792167de7907f09cc8c897c8a3f7e211 diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index 241a21636..a2a44a238 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -95,6 +95,12 @@ namespace VISCOUS_3D enum UIndex { U_TGT = 1, U_SRC, LEN_TGT }; const double theMinSmoothCosin = 0.1; + const double theSmoothThickToElemSizeRatio = 0.3; + + bool needSmoothing( double cosin, double tgtThick, double elemSize ) + { + return cosin * tgtThick > theSmoothThickToElemSizeRatio * elemSize; + } /*! * \brief SMESH_ProxyMesh computed by _ViscousBuilder for a SOLID. @@ -295,30 +301,7 @@ namespace VISCOUS_3D double lenDelta(double len) const { return _k * ( _r + len ); } double lenDeltaByDist(double dist) const { return dist * _h2lenRatio; } }; - struct _LayerEdge; - //-------------------------------------------------------------------------------- - /*! - * Structure used to smooth a _LayerEdge (master) based on an EDGE. - */ - struct _2NearEdges - { - // target nodes of 2 neighbour _LayerEdge's based on the same EDGE - const SMDS_MeshNode* _nodes[2]; - // vectors from source nodes of 2 _LayerEdge's to the source node of master _LayerEdge - //gp_XYZ _vec[2]; - double _wgt[2]; // weights of _nodes - _LayerEdge* _edges[2]; - - // normal to plane passing through _LayerEdge._normal and tangent of EDGE - gp_XYZ* _plnNorm; - - _2NearEdges() { _nodes[0]=_nodes[1]=0; _plnNorm = 0; } - void reverse() { - std::swap( _nodes[0], _nodes[1] ); - std::swap( _wgt [0], _wgt [1] ); - std::swap( _edges[0], _edges[1] ); - } - }; + struct _2NearEdges; //-------------------------------------------------------------------------------- /*! * \brief Edge normal to surface, connecting a node on solid surface (_nodes[0]) @@ -368,6 +351,7 @@ namespace VISCOUS_3D double& dist, const double& epsilon) const; gp_Ax1 LastSegment(double& segLen) const; + gp_XY LastUV( const TopoDS_Face& F ) const; bool IsOnEdge() const { return _2neibors; } gp_XYZ Copy( _LayerEdge& other, SMESH_MesherHelper& helper ); void SetCosin( double cosin ); @@ -380,6 +364,31 @@ namespace VISCOUS_3D return cmpNodes ? ( e1->_nodes[0]->GetID() < e2->_nodes[0]->GetID()) : ( e1 < e2 ); } }; + struct _LayerEdge; + //-------------------------------------------------------------------------------- + /*! + * Structure used to smooth a _LayerEdge based on an EDGE. + */ + struct _2NearEdges + { + double _wgt [2]; // weights of _nodes + _LayerEdge* _edges[2]; + + // normal to plane passing through _LayerEdge._normal and tangent of EDGE + gp_XYZ* _plnNorm; + + _2NearEdges() { _edges[0]=_edges[1]=0; _plnNorm = 0; } + const SMDS_MeshNode* tgtNode(bool is2nd) { + return _edges[is2nd] ? _edges[is2nd]->_nodes.back() : 0; + } + const SMDS_MeshNode* srcNode(bool is2nd) { + return _edges[is2nd] ? _edges[is2nd]->_nodes[0] : 0; + } + void reverse() { + std::swap( _wgt [0], _wgt [1] ); + std::swap( _edges[0], _edges[1] ); + } + }; //-------------------------------------------------------------------------------- /*! * \brief Convex FACE whose radius of curvature is less than the thickness of @@ -422,10 +431,11 @@ namespace VISCOUS_3D set _reversedFaceIds; set _ignoreFaceIds; // WOL FACEs and FACEs of other SOLIDS - double _stepSize, _stepSizeCoeff; + double _stepSize, _stepSizeCoeff, _geomSize; const SMDS_MeshNode* _stepSizeNodes[2]; - TNode2Edge _n2eMap; + TNode2Edge _n2eMap; // nodes and _LayerEdge's based on them + // map to find _n2eMap of another _SolidData by a shrink shape shared by two _SolidData's map< TGeomID, TNode2Edge* > _s2neMap; // edges of _n2eMap. We keep same data in two containers because @@ -441,8 +451,9 @@ namespace VISCOUS_3D // Convex FACEs whose radius of curvature is less than the thickness of layers map< TGeomID, _ConvexFace > _convexFaces; - // FACE's WOL, srink on which is forbiden due to algo on the adjacent SOLID - set< TGeomID > _noShrinkFaces; + // shapes (EDGEs and VERTEXes) srink from which is forbiden due to collisions with + // the adjacent SOLID + set< TGeomID > _noShrinkShapes; // to -- for analytic smooth map< TGeomID,Handle(Geom_Curve)> _edge2curve; @@ -453,7 +464,7 @@ namespace VISCOUS_3D double _epsilon; // precision for SegTriaInter() - int _index; // for debug + TGeomID _index; // SOLID id, for debug _SolidData(const TopoDS_Shape& s=TopoDS_Shape(), const StdMeshers_ViscousLayers* h=0, @@ -584,6 +595,7 @@ namespace VISCOUS_3D const bool toSort = false); void findSimplexTestEdges( _SolidData& data, vector< vector<_LayerEdge*> >& edgesByGeom); + void computeGeomSize( _SolidData& data ); bool sortEdges( _SolidData& data, vector< vector<_LayerEdge*> >& edgesByGeom); void limitStepSizeByCurvature( _SolidData& data ); @@ -608,6 +620,7 @@ namespace VISCOUS_3D bool prepareEdgeToShrink( _LayerEdge& edge, const TopoDS_Face& F, SMESH_MesherHelper& helper, const SMESHDS_SubMesh* faceSubMesh ); + void restoreNoShrink( _LayerEdge& edge ) const; void fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& helper, const bool is2D, @@ -616,7 +629,7 @@ namespace VISCOUS_3D bool addBoundaryElements(); bool error( const string& text, int solidID=-1 ); - SMESHDS_Mesh* getMeshDS() { return _mesh->GetMeshDS(); } + SMESHDS_Mesh* getMeshDS() const { return _mesh->GetMeshDS(); } // debug void makeGroupOfLE(); @@ -659,7 +672,8 @@ namespace VISCOUS_3D virtual SMDSAbs_ElementType GetType() const { return SMDSAbs_Face; } virtual vtkIdType GetVtkType() const { return -1; } virtual SMDSAbs_EntityType GetEntityType() const { return SMDSEntity_Last; } - virtual SMDSAbs_GeometryType GetGeomType() const { return SMDSGeom_TRIANGLE; } + virtual SMDSAbs_GeometryType GetGeomType() const + { return _nn.size() == 3 ? SMDSGeom_TRIANGLE : SMDSGeom_QUADRANGLE; } virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType) const { return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));} }; @@ -717,6 +731,7 @@ namespace VISCOUS_3D return _surface->Value( uv.X(), uv.Y() ).XYZ(); } }; + } // namespace VISCOUS_3D @@ -842,22 +857,32 @@ namespace gp_Vec dir; double f,l; gp_Pnt p; Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l ); + if ( c.IsNull() ) return gp_XYZ( 1e100, 1e100, 1e100 ); double u = helper.GetNodeU( E, atNode ); c->D1( u, p, dir ); return dir.XYZ(); } //-------------------------------------------------------------------------------- + gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV, + const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok, + double* cosin=0); + //-------------------------------------------------------------------------------- gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Edge& fromE, const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok) { + double f,l; + Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l ); + if ( c.IsNull() ) + { + TopoDS_Vertex v = helper.IthVertex( 0, fromE ); + return getFaceDir( F, v, node, helper, ok ); + } gp_XY uv = helper.GetNodeUV( F, node, 0, &ok ); Handle(Geom_Surface) surface = BRep_Tool::Surface( F ); gp_Pnt p; gp_Vec du, dv, norm; surface->D1( uv.X(),uv.Y(), p, du,dv ); norm = du ^ dv; - double f,l; - Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l ); double u = helper.GetNodeU( fromE, node, 0, &ok ); c->D1( u, p, du ); TopAbs_Orientation o = helper.GetSubShapeOri( F.Oriented(TopAbs_FORWARD), fromE); @@ -881,7 +906,7 @@ namespace //-------------------------------------------------------------------------------- gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV, const SMDS_MeshNode* node, SMESH_MesherHelper& helper, - bool& ok, double* cosin=0) + bool& ok, double* cosin) { TopoDS_Face faceFrw = F; faceFrw.Orientation( TopAbs_FORWARD ); @@ -1012,6 +1037,59 @@ namespace } return false; } + //================================================================================ + /*! + * \brief Computes mimimal distance of face in-FACE nodes from an EDGE + * \param [in] face - the mesh face to treat + * \param [in] nodeOnEdge - a node on the EDGE + * \param [out] faceSize - the computed distance + * \return bool - true if faceSize computed + */ + //================================================================================ + + bool getDistFromEdge( const SMDS_MeshElement* face, + const SMDS_MeshNode* nodeOnEdge, + double & faceSize ) + { + faceSize = Precision::Infinite(); + bool done = false; + + int nbN = face->NbCornerNodes(); + int iOnE = face->GetNodeIndex( nodeOnEdge ); + int iNext[2] = { SMESH_MesherHelper::WrapIndex( iOnE+1, nbN ), + SMESH_MesherHelper::WrapIndex( iOnE-1, nbN ) }; + const SMDS_MeshNode* nNext[2] = { face->GetNode( iNext[0] ), + face->GetNode( iNext[1] ) }; + gp_XYZ segVec, segEnd = SMESH_TNodeXYZ( nodeOnEdge ); // segment on EDGE + double segLen = -1.; + // look for two neighbor not in-FACE nodes of face + for ( int i = 0; i < 2; ++i ) + { + if ( nNext[i]->GetPosition()->GetDim() != 2 && + nNext[i]->GetID() < nodeOnEdge->GetID() ) + { + // look for an in-FACE node + for ( int iN = 0; iN < nbN; ++iN ) + { + if ( iN == iOnE || iN == iNext[i] ) + continue; + SMESH_TNodeXYZ pInFace = face->GetNode( iN ); + gp_XYZ v = pInFace - segEnd; + if ( segLen < 0 ) + { + segVec = SMESH_TNodeXYZ( nNext[i] ) - segEnd; + segLen = segVec.Modulus(); + } + double distToSeg = v.Crossed( segVec ).Modulus() / segLen; + faceSize = Min( faceSize, distToSeg ); + done = true; + } + segLen = -1; + } + } + return done; + } + //-------------------------------------------------------------------------------- // DEBUG. Dump intermediate node positions into a python script // HOWTO use: run python commands written in a console to see @@ -1092,19 +1170,36 @@ _ViscousBuilder::_ViscousBuilder() bool _ViscousBuilder::error(const string& text, int solidId ) { + const string prefix = string("Viscous layers builder: "); _error->myName = COMPERR_ALGO_FAILED; - _error->myComment = string("Viscous layers builder: ") + text; + _error->myComment = prefix + text; if ( _mesh ) { SMESH_subMesh* sm = _mesh->GetSubMeshContaining( solidId ); if ( !sm && !_sdVec.empty() ) - sm = _mesh->GetSubMeshContaining( _sdVec[0]._index ); + sm = _mesh->GetSubMeshContaining( solidId = _sdVec[0]._index ); if ( sm && sm->GetSubShape().ShapeType() == TopAbs_SOLID ) { SMESH_ComputeErrorPtr& smError = sm->GetComputeError(); if ( smError && smError->myAlgo ) _error->myAlgo = smError->myAlgo; smError = _error; + sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); + } + // set KO to all solids + for ( size_t i = 0; i < _sdVec.size(); ++i ) + { + if ( _sdVec[i]._index == solidId ) + continue; + sm = _mesh->GetSubMesh( _sdVec[i]._solid ); + if ( !sm->IsEmpty() ) + continue; + SMESH_ComputeErrorPtr& smError = sm->GetComputeError(); + if ( !smError || smError->IsOK() ) + { + smError = SMESH_ComputeError::New( COMPERR_ALGO_FAILED, prefix + "failed"); + sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); + } } } makeGroupOfLE(); // debug @@ -1360,6 +1455,7 @@ bool _ViscousBuilder::findFacesWithLayers() continue; // nothing interesting TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ]; // check presence of layers on fWOL within an adjacent SOLID + bool collision = false; PShapeIteratorPtr sIt = helper.GetAncestors( fWOL, *_mesh, TopAbs_SOLID ); while ( const TopoDS_Shape* solid = sIt->next() ) if ( !solid->IsSame( _sdVec[i]._solid )) @@ -1368,8 +1464,9 @@ bool _ViscousBuilder::findFacesWithLayers() int iFace = getMeshDS()->ShapeToIndex( fWOL ); if ( iSolid > 0 && !_sdVec[ iSolid-1 ]._ignoreFaceIds.count( iFace )) { - _sdVec[i]._noShrinkFaces.insert( iFace ); - fWOL.Nullify(); + //_sdVec[i]._noShrinkShapes.insert( iFace ); + //fWOL.Nullify(); + collision = true; } } // add edge to maps @@ -1377,6 +1474,12 @@ bool _ViscousBuilder::findFacesWithLayers() { TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge ); _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL )); + if ( collision ) + { + // _shrinkShape2Shape will be used to temporary inflate _LayerEdge's based + // on the edge but shrink won't be performed + _sdVec[i]._noShrinkShapes.insert( edgeInd ); + } } } } @@ -1385,54 +1488,70 @@ bool _ViscousBuilder::findFacesWithLayers() set< string > notSupportAlgos; notSupportAlgos.insert("Hexa_3D"); for ( size_t i = 0; i < _sdVec.size(); ++i ) { - TopTools_MapOfShape noShrinkVertices; map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin(); for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f ) { const TopoDS_Shape& fWOL = e2f->second; - TGeomID edgeID = e2f->first; + const TGeomID edgeID = e2f->first; bool notShrinkFace = false; PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID); - while ( soIt->more()) + while ( soIt->more() ) { const TopoDS_Shape* solid = soIt->next(); if ( _sdVec[i]._solid.IsSame( *solid )) continue; SMESH_Algo* algo = _mesh->GetGen()->GetAlgo( *_mesh, *solid ); if ( !algo || !notSupportAlgos.count( algo->GetName() )) continue; notShrinkFace = true; - for ( size_t j = 0; j < _sdVec.size(); ++j ) + size_t iSolid = 0; + for ( ; iSolid < _sdVec.size(); ++iSolid ) { - if ( _sdVec[j]._solid.IsSame( *solid ) ) - if ( _sdVec[j]._shrinkShape2Shape.count( edgeID )) + if ( _sdVec[iSolid]._solid.IsSame( *solid ) ) { + if ( _sdVec[iSolid]._shrinkShape2Shape.count( edgeID )) notShrinkFace = false; + break; + } } - } - if ( notShrinkFace ) - { - _sdVec[i]._noShrinkFaces.insert( getMeshDS()->ShapeToIndex( fWOL )); - for ( TopExp_Explorer vExp( fWOL, TopAbs_VERTEX ); vExp.More(); vExp.Next() ) - noShrinkVertices.Add( vExp.Current() ); - } - } - // erase from _shrinkShape2Shape all srink EDGE's of a SOLID connected - // to the found not shrinked fWOL's - e2f = _sdVec[i]._shrinkShape2Shape.begin(); - for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ) - { - TGeomID edgeID = e2f->first; - TopoDS_Vertex VV[2]; - TopExp::Vertices( TopoDS::Edge( getMeshDS()->IndexToShape( edgeID )),VV[0],VV[1]); - if ( noShrinkVertices.Contains( VV[0] ) || noShrinkVertices.Contains( VV[1] )) - { - _sdVec[i]._noShrinkFaces.insert( getMeshDS()->ShapeToIndex( e2f->second )); - _sdVec[i]._shrinkShape2Shape.erase( e2f++ ); - } - else - { - e2f++; - } - } - } + if ( notShrinkFace ) + { + _sdVec[i]._noShrinkShapes.insert( edgeID ); + + // add VERTEXes of the edge in _noShrinkShapes + TopoDS_Shape edge = getMeshDS()->IndexToShape( edgeID ); + for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() ) + _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( vIt.Value() )); + + // check if there is a collision with to-shrink-from EDGEs in iSolid + if ( iSolid == _sdVec.size() ) + continue; // no VL in the solid + shapes.Clear(); + TopExp::MapShapes( fWOL, TopAbs_EDGE, shapes); + for ( int iE = 1; iE <= shapes.Extent(); ++iE ) + { + const TopoDS_Edge& E = TopoDS::Edge( shapes( iE )); + const TGeomID eID = getMeshDS()->ShapeToIndex( E ); + if ( eID == edgeID || + !_sdVec[iSolid]._shrinkShape2Shape.count( eID ) || + _sdVec[i]._noShrinkShapes.count( eID )) + continue; + for ( int is1st = 0; is1st < 2; ++is1st ) + { + TopoDS_Vertex V = helper.IthVertex( is1st, E ); + if ( _sdVec[i]._noShrinkShapes.count( getMeshDS()->ShapeToIndex( V ) )) + { + // _sdVec[i]._noShrinkShapes.insert( eID ); + // V = helper.IthVertex( !is1st, E ); + // _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( V )); + //iE = 0; // re-start the loop on EDGEs of fWOL + return error("No way to make a conformal mesh with " + "the given set of faces with layers", _sdVec[i]._index); + } + } + } + } + + } // while ( soIt->more() ) + } // loop on _sdVec[i]._shrinkShape2Shape + } // loop on _sdVec to fill in _SolidData::_noShrinkShapes // Find the SHAPE along which to inflate _LayerEdge based on VERTEX @@ -1454,8 +1573,8 @@ bool _ViscousBuilder::findFacesWithLayers() { totalNbFaces++; const int fID = getMeshDS()->ShapeToIndex( *f ); - if ( _sdVec[i]._ignoreFaceIds.count ( fID ) && - !_sdVec[i]._noShrinkFaces.count( fID )) + if ( _sdVec[i]._ignoreFaceIds.count ( fID ) /*&& + !_sdVec[i]._noShrinkShapes.count( fID )*/) facesWOL.push_back( *f ); } } @@ -1533,15 +1652,14 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) { // get all sub-shapes to make layers on set subIds, faceIds; - subIds = data._noShrinkFaces; + subIds = data._noShrinkShapes; TopExp_Explorer exp( data._solid, TopAbs_FACE ); for ( ; exp.More(); exp.Next() ) { SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() ); - if ( ! data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() ))) + if ( ! data._ignoreFaceIds.count( fSubM->GetId() )) faceIds.insert( fSubM->GetId() ); - SMESH_subMeshIteratorPtr subIt = - fSubM->getDependsOnIterator(/*includeSelf=*/true, /*complexShapeFirst=*/false); + SMESH_subMeshIteratorPtr subIt = fSubM->getDependsOnIterator(/*includeSelf=*/true); while ( subIt->more() ) subIds.insert( subIt->next()->GetId() ); } @@ -1574,7 +1692,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) SMESH_MesherHelper helper( *_mesh ); helper.SetSubShape( data._solid ); - helper.SetElementsOnShape(true); + helper.SetElementsOnShape( true ); vector< const SMDS_MeshNode*> newNodes; // of a mesh face TNode2Edge::iterator n2e2; @@ -1596,28 +1714,51 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) while ( eIt->more() ) { const SMDS_MeshElement* face = eIt->next(); + double faceMaxCosin = -1; + _LayerEdge* maxCosinEdge = 0; + int nbDegenNodes = 0; + newNodes.resize( face->NbCornerNodes() ); - double faceMaxCosin = -1; - _LayerEdge* maxCosinEdge = 0; - for ( int i = 0 ; i < face->NbCornerNodes(); ++i ) + for ( size_t i = 0 ; i < newNodes.size(); ++i ) { - const SMDS_MeshNode* n = face->GetNode(i); + const SMDS_MeshNode* n = face->GetNode( i ); + const int shapeID = n->getshapeId(); + const bool onDegenShap = helper.IsDegenShape( shapeID ); + const bool onDegenEdge = ( onDegenShap && n->GetPosition()->GetDim() == 1 ); + if ( onDegenShap ) + { + if ( onDegenEdge ) + { + // substitute n on a degenerated EDGE with a node on a corresponding VERTEX + const TopoDS_Shape& E = getMeshDS()->IndexToShape( shapeID ); + TopoDS_Vertex V = helper.IthVertex( 0, TopoDS::Edge( E )); + if ( const SMDS_MeshNode* vN = SMESH_Algo::VertexNode( V, getMeshDS() )) { + n = vN; + nbDegenNodes++; + } + } + else + { + nbDegenNodes++; + } + } TNode2Edge::iterator n2e = data._n2eMap.insert( make_pair( n, (_LayerEdge*)0 )).first; if ( !(*n2e).second ) { // add a _LayerEdge _LayerEdge* edge = new _LayerEdge(); - n2e->second = edge; edge->_nodes.push_back( n ); - const int shapeID = n->getshapeId(); + n2e->second = edge; edgesByGeom[ shapeID ].push_back( edge ); + const bool noShrink = data._noShrinkShapes.count( shapeID ); 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 = data._s2neMap.find( shapeID )) != data._s2neMap.end() && - ( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end()) + if (( !noShrink ) && + ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE ) && + (( s2ne = data._s2neMap.find( shapeID )) != data._s2neMap.end() ) && + (( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end() )) { _LayerEdge* foundEdge = (*n2e2).second; gp_XYZ lastPos = edge->Copy( *foundEdge, helper ); @@ -1628,7 +1769,10 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) } else { - edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() )); + if ( !noShrink ) + { + edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() )); + } if ( !setEdgeData( *edge, subIds, helper, data )) return false; } @@ -1641,7 +1785,13 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) } } newNodes[ i ] = n2e->second->_nodes.back(); + + if ( onDegenEdge ) + data._n2eMap.insert( make_pair( face->GetNode( i ), n2e->second )); } + if ( newNodes.size() - nbDegenNodes < 2 ) + continue; + // create a temporary face const SMDS_MeshElement* newFace = new _TmpMeshFace( newNodes, --_tmpFaceID, face->getshapeId() ); @@ -1650,6 +1800,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) // compute inflation step size by min size of element on a convex surface if ( faceMaxCosin > theMinSmoothCosin ) limitStepSize( data, face, maxCosinEdge ); + } // loop on 2D elements on a FACE } // loop on FACEs of a SOLID @@ -1664,28 +1815,61 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) // limit data._stepSize depending on surface curvature and fill data._convexFaces limitStepSizeByCurvature( data ); // !!! it must be before node substitution in _Simplex - // Set target nodes into _Simplex and _2NearEdges of _LayerEdge's + // Set target nodes into _Simplex and _LayerEdge's to _2NearEdges TNode2Edge::iterator n2e; + const SMDS_MeshNode* nn[2]; for ( size_t i = 0; i < data._edges.size(); ++i ) { - if ( data._edges[i]->IsOnEdge()) + _LayerEdge* edge = data._edges[i]; + if ( edge->IsOnEdge() ) + { + // get neighbor nodes + bool hasData = ( edge->_2neibors->_edges[0] ); + if ( hasData ) // _LayerEdge is a copy of another one + { + nn[0] = edge->_2neibors->srcNode(0); + nn[1] = edge->_2neibors->srcNode(1); + } + else if ( !findNeiborsOnEdge( edge, nn[0],nn[1], data )) + { + return false; + } + // set neighbor _LayerEdge's for ( int j = 0; j < 2; ++j ) { - if ( data._edges[i]->_nodes.back()->NbInverseElements(SMDSAbs_Volume) > 0 ) - break; // _LayerEdge is shared by two _SolidData's - const SMDS_MeshNode* & n = data._edges[i]->_2neibors->_nodes[j]; - if (( n2e = data._n2eMap.find( n )) == data._n2eMap.end() ) + if (( n2e = data._n2eMap.find( nn[j] )) == data._n2eMap.end() ) return error("_LayerEdge not found by src node", data._index); - n = (*n2e).second->_nodes.back(); - data._edges[i]->_2neibors->_edges[j] = n2e->second; + edge->_2neibors->_edges[j] = n2e->second; } - //else - for ( size_t j = 0; j < data._edges[i]->_simplices.size(); ++j ) + if ( !hasData ) + edge->SetDataByNeighbors( nn[0], nn[1], helper); + } + + for ( size_t j = 0; j < edge->_simplices.size(); ++j ) { - _Simplex& s = data._edges[i]->_simplices[j]; + _Simplex& s = edge->_simplices[j]; s._nNext = data._n2eMap[ s._nNext ]->_nodes.back(); s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back(); } + + // For an _LayerEdge on a degenerated EDGE, copy some data from + // a corresponding _LayerEdge on a VERTEX + // (issue 52453, pb on a downloaded SampleCase2-Tet-netgen-mephisto.hdf) + if ( helper.IsDegenShape( edge->_nodes[0]->getshapeId() )) + { + // Generally we should not get here + const TopoDS_Shape& E = getMeshDS()->IndexToShape( edge->_nodes[0]->getshapeId() ); + if ( E.ShapeType() != TopAbs_EDGE ) + continue; + TopoDS_Vertex V = helper.IthVertex( 0, TopoDS::Edge( E )); + const SMDS_MeshNode* vN = SMESH_Algo::VertexNode( V, getMeshDS() ); + if (( n2e = data._n2eMap.find( vN )) == data._n2eMap.end() ) + continue; + const _LayerEdge* vEdge = n2e->second; + edge->_normal = vEdge->_normal; + edge->_lenFactor = vEdge->_lenFactor; + edge->_cosin = vEdge->_cosin; + } } dumpFunctionEnd(); @@ -1767,6 +1951,10 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) for ( ; face.More(); face.Next() ) { const TopoDS_Face& F = TopoDS::Face( face.Current() ); + SMESH_subMesh * sm = _mesh->GetSubMesh( F ); + const TGeomID faceID = sm->GetId(); + if ( data._ignoreFaceIds.count( faceID )) continue; + BRepAdaptor_Surface surface( F, false ); surfProp.SetSurface( surface ); @@ -1774,9 +1962,6 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) int iBeg, iEnd; _ConvexFace cnvFace; - SMESH_subMesh * sm = _mesh->GetSubMesh( F ); - const TGeomID faceID = sm->GetId(); - if ( data._ignoreFaceIds.count( faceID )) continue; const double oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. ); SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true); while ( smIt->more() ) @@ -1791,8 +1976,8 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) continue; // check concavity and curvature and limit data._stepSize int nbLEdges = iEnd - iBeg; - int step = Max( 1, nbLEdges / nbTestPnt ); - for ( ; iBeg < iEnd; iBeg += step ) + int iStep = Max( 1, nbLEdges / nbTestPnt ); + for ( ; iBeg < iEnd; iBeg += iStep ) { gp_XY uv = helper.GetNodeUV( F, data._edges[ iBeg ]->_nodes[0] ); surfProp.SetParameters( uv.X(), uv.Y() ); @@ -1882,8 +2067,12 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) bool _ViscousBuilder::sortEdges( _SolidData& data, vector< vector<_LayerEdge*> >& edgesByGeom) { + // define allowed thickness + computeGeomSize( data ); // compute data._geomSize + const double tgtThick = Min( 0.5 * data._geomSize, data._hyp->GetTotalThickness() ); + // Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's - // boundry inclined at a sharp angle to the shape + // boundry inclined to the shape at a sharp angle list< TGeomID > shapesToSmooth; @@ -1900,31 +2089,33 @@ bool _ViscousBuilder::sortEdges( _SolidData& data, { case TopAbs_EDGE: { - bool isShrinkEdge = !eS[0]->_sWOL.IsNull(); + if ( SMESH_Algo::isDegenerated( TopoDS::Edge( S ))) + break; + //bool isShrinkEdge = !eS[0]->_sWOL.IsNull(); for ( TopoDS_Iterator vIt( S ); vIt.More() && !needSmooth; vIt.Next() ) { TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() ); vector<_LayerEdge*>& eV = edgesByGeom[ iV ]; if ( eV.empty() ) continue; - // double cosin = eV[0]->_cosin; - // bool badCosin = - // ( !eV[0]->_sWOL.IsNull() && ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE || !isShrinkEdge)); - // if ( badCosin ) - // { - // gp_Vec dir1, dir2; - // if ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE ) - // dir1 = getEdgeDir( TopoDS::Edge( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() )); - // else - // dir1 = getFaceDir( TopoDS::Face( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ), - // eV[0]->_nodes[0], helper, ok); - // dir2 = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() )); - // double angle = dir1.Angle( dir2 ); - // cosin = cos( angle ); - // } gp_Vec eDir = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() )); double angle = eDir.Angle( eV[0]->_normal ); double cosin = Cos( angle ); - needSmooth = ( cosin > theMinSmoothCosin ); + if ( cosin > theMinSmoothCosin ) + { + // compare tgtThick with the length of an end segment + SMDS_ElemIteratorPtr eIt = eV[0]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Edge); + while ( eIt->more() ) + { + const SMDS_MeshElement* endSeg = eIt->next(); + if ( endSeg->getshapeId() == iS ) + { + double segLen = + SMESH_TNodeXYZ( endSeg->GetNode(0) ).Distance( endSeg->GetNode(1 )); + needSmooth = needSmoothing( cosin, tgtThick, segLen ); + break; + } + } + } } break; } @@ -1935,25 +2126,38 @@ bool _ViscousBuilder::sortEdges( _SolidData& data, TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() ); vector<_LayerEdge*>& eE = edgesByGeom[ iE ]; if ( eE.empty() ) continue; - if ( eE[0]->_sWOL.IsNull() ) - { - for ( size_t i = 0; i < eE.size() && !needSmooth; ++i ) - needSmooth = ( eE[i]->_cosin > theMinSmoothCosin ); - } - else + // TopLoc_Location loc; + // Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face( S ), loc ); + // bool isPlane = GeomLib_IsPlanarSurface( surface ).IsPlanar(); + //if ( eE[0]->_sWOL.IsNull() ) { - const TopoDS_Face& F1 = TopoDS::Face( S ); - const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL ); - const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() ); + double faceSize; for ( size_t i = 0; i < eE.size() && !needSmooth; ++i ) - { - gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok ); - gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok ); - double angle = dir1.Angle( dir2 ); - double cosin = cos( angle ); - needSmooth = ( cosin > theMinSmoothCosin ); - } + if ( eE[i]->_cosin > theMinSmoothCosin ) + { + SMDS_ElemIteratorPtr fIt = eE[i]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face); + while ( fIt->more() && !needSmooth ) + { + const SMDS_MeshElement* face = fIt->next(); + if ( getDistFromEdge( face, eE[i]->_nodes[0], faceSize )) + needSmooth = needSmoothing( eE[i]->_cosin, tgtThick, faceSize ); + } + } } + // else + // { + // const TopoDS_Face& F1 = TopoDS::Face( S ); + // const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL ); + // const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() ); + // for ( size_t i = 0; i < eE.size() && !needSmooth; ++i ) + // { + // gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok ); + // gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok ); + // double angle = dir1.Angle( ); + // double cosin = cos( angle ); + // needSmooth = ( cosin > theMinSmoothCosin ); + // } + // } } break; } @@ -1961,6 +2165,7 @@ bool _ViscousBuilder::sortEdges( _SolidData& data, continue; default:; } + if ( needSmooth ) { if ( S.ShapeType() == TopAbs_EDGE ) shapesToSmooth.push_front( iS ); @@ -2165,14 +2370,16 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, if ( edge._sWOL.ShapeType() == TopAbs_EDGE ) { double u = helper.GetNodeU( TopoDS::Edge( edge._sWOL ), node, 0, &normOK ); - edge._pos.push_back( gp_XYZ( u, 0, 0)); - getMeshDS()->SetNodeOnEdge( tgtNode, TopoDS::Edge( edge._sWOL ), u ); + edge._pos.push_back( gp_XYZ( u, 0, 0 )); + if ( edge._nodes.size() > 1 ) + getMeshDS()->SetNodeOnEdge( tgtNode, TopoDS::Edge( edge._sWOL ), u ); } else // TopAbs_FACE { gp_XY uv = helper.GetNodeUV( TopoDS::Face( edge._sWOL ), node, 0, &normOK ); edge._pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0)); - getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( edge._sWOL ), uv.X(), uv.Y() ); + if ( edge._nodes.size() > 1 ) + getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( edge._sWOL ), uv.X(), uv.Y() ); } } else @@ -2202,14 +2409,14 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, { edge._2neibors = new _2NearEdges; // target node instead of source ones will be set later - if ( ! findNeiborsOnEdge( &edge, - edge._2neibors->_nodes[0], - edge._2neibors->_nodes[1], - data)) - return false; - edge.SetDataByNeighbors( edge._2neibors->_nodes[0], - edge._2neibors->_nodes[1], - helper); + // if ( ! findNeiborsOnEdge( &edge, + // edge._2neibors->_nodes[0], + // edge._2neibors->_nodes[1], + // data)) + // return false; + // edge.SetDataByNeighbors( edge._2neibors->_nodes[0], + // edge._2neibors->_nodes[1], + // helper); } edge.SetCosin( edge._cosin ); // to update edge._lenFactor @@ -2286,9 +2493,19 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node, Handle(Geom_Surface) surface = BRep_Tool::Surface( face ); int pointKind = GeomLib::NormEstim( surface, uv, 1e-5, normal ); enum { REGULAR = 0, QUASYSINGULAR, CONICAL, IMPOSSIBLE }; + + if ( pointKind == IMPOSSIBLE && + node->GetPosition()->GetDim() == 2 ) // node inside the FACE + { + // probably NormEstim() failed due to a too high tolerance + pointKind = GeomLib::NormEstim( surface, uv, 1e-20, normal ); + isOK = ( pointKind < IMPOSSIBLE ); + } if ( pointKind < IMPOSSIBLE ) { - if ( pointKind != REGULAR && !shiftInside ) + if ( pointKind != REGULAR && + !shiftInside && + node->GetPosition()->GetDim() < 2 ) // FACE boundary { gp_XYZ normShift = getFaceNormal( node, face, helper, isOK, /*shiftInside=*/true ); if ( normShift * normal.XYZ() < 0. ) @@ -2296,7 +2513,8 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node, } isOK = true; } - else // hard singularity, to call with shiftInside=true ? + + if ( !isOK ) // hard singularity, to call with shiftInside=true ? { const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face ); @@ -2309,7 +2527,9 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node, isOK = SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) normal.XYZ(), /*normalized=*/true ); if ( isOK ) { - if ( helper.IsReversedSubMesh( face )) + TopoDS_Face ff = face; + ff.Orientation( TopAbs_FORWARD ); + if ( helper.IsReversedSubMesh( ff )) normal.Reverse(); break; } @@ -2410,11 +2630,10 @@ bool _ViscousBuilder::findNeiborsOnEdge(const _LayerEdge* edge, _SolidData& data) { const SMDS_MeshNode* node = edge->_nodes[0]; - const int shapeInd = node->getshapeId(); - SMESHDS_SubMesh* edgeSM = 0; + const int shapeInd = node->getshapeId(); + SMESHDS_SubMesh* edgeSM = 0; if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE ) { - edgeSM = getMeshDS()->MeshElements( shapeInd ); if ( !edgeSM || edgeSM->NbElements() == 0 ) return error(SMESH_Comment("Not meshed EDGE ") << shapeInd, data._index); @@ -2462,11 +2681,11 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1, // Set _curvature - double sumLen = vec1.Modulus() + vec2.Modulus(); + double sumLen = vec1.Modulus() + vec2.Modulus(); _2neibors->_wgt[0] = 1 - vec1.Modulus() / sumLen; _2neibors->_wgt[1] = 1 - vec2.Modulus() / sumLen; double avgNormProj = 0.5 * ( _normal * vec1 + _normal * vec2 ); - double avgLen = 0.5 * ( vec1.Modulus() + vec2.Modulus() ); + double avgLen = 0.5 * ( vec1.Modulus() + vec2.Modulus() ); if ( _curvature ) delete _curvature; _curvature = _Curvature::New( avgNormProj, avgLen ); // if ( _curvature ) @@ -2480,10 +2699,13 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1, if ( _sWOL.IsNull() ) { TopoDS_Shape S = helper.GetSubShapeByNode( _nodes[0], helper.GetMeshDS() ); - gp_XYZ dirE = getEdgeDir( TopoDS::Edge( S ), _nodes[0], helper ); + TopoDS_Edge E = TopoDS::Edge( S ); + // if ( SMESH_Algo::isDegenerated( E )) + // return; + gp_XYZ dirE = getEdgeDir( E, _nodes[0], helper ); gp_XYZ plnNorm = dirE ^ _normal; - double proj0 = plnNorm * vec1; - double proj1 = plnNorm * vec2; + double proj0 = plnNorm * vec1; + double proj1 = plnNorm * vec2; if ( fabs( proj0 ) > 1e-10 || fabs( proj1 ) > 1e-10 ) { if ( _2neibors->_plnNorm ) delete _2neibors->_plnNorm; @@ -2650,6 +2872,31 @@ void _ViscousBuilder::makeGroupOfLE() #endif } +//================================================================================ +/*! + * \brief Find maximal _LayerEdge length (layer thickness) limited by geometry + */ +//================================================================================ + +void _ViscousBuilder::computeGeomSize( _SolidData& data ) +{ + data._geomSize = Precision::Infinite(); + double intersecDist; + auto_ptr searcher + ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), + data._proxyMesh->GetFaces( data._solid )) ); + + TNode2Edge::iterator n2e = data._n2eMap.begin(), n2eEnd = data._n2eMap.end(); + for ( ; n2e != n2eEnd; ++n2e ) + { + _LayerEdge* edge = n2e->second; + if ( edge->IsOnEdge() ) continue; + edge->FindIntersection( *searcher, intersecDist, data._epsilon ); + if ( data._geomSize > intersecDist && intersecDist > 0 ) + data._geomSize = intersecDist; + } +} + //================================================================================ /*! * \brief Increase length of _LayerEdge's to reach the required thickness of layers @@ -2662,19 +2909,8 @@ 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; - auto_ptr searcher - ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), - data._proxyMesh->GetFaces( data._solid )) ); - for ( size_t i = 0; i < data._edges.size(); ++i ) - { - if ( data._edges[i]->IsOnEdge() ) continue; - data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon ); - if ( geomSize > intersecDist && intersecDist > 0 ) - geomSize = intersecDist; - } - if ( data._stepSize > 0.3 * geomSize ) - limitStepSize( data, 0.3 * geomSize ); + if ( data._stepSize > 0.3 * data._geomSize ) + limitStepSize( data, 0.3 * data._geomSize ); const double tgtThick = data._hyp->GetTotalThickness(); if ( data._stepSize > tgtThick ) @@ -2683,7 +2919,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) if ( data._stepSize < 1. ) data._epsilon = data._stepSize * 1e-7; - debugMsg( "-- geomSize = " << geomSize << ", stepSize = " << data._stepSize ); + debugMsg( "-- geomSize = " << data._geomSize << ", stepSize = " << data._stepSize ); double avgThick = 0, curThick = 0, distToIntersection = Precision::Infinite(); int nbSteps = 0, nbRepeats = 0; @@ -2761,6 +2997,22 @@ bool _ViscousBuilder::inflate(_SolidData& data) " average reached thickness is " << avgThick )); } + + // Restore position of src nodes moved by infaltion on _noShrinkShapes + dumpFunction(SMESH_Comment("restoNoShrink_So")<_nodes.size() == 1 ) + for ( ; iBeg < iEnd; ++iBeg ) + { + restoreNoShrink( *data._edges[ iBeg ] ); + } + } + dumpFunctionEnd(); + return true; } @@ -2980,8 +3232,8 @@ Handle(Geom_Curve) _SolidData::CurveForSmooth( const TopoDS_Edge& E, bndBox.Add( SMESH_TNodeXYZ( nIt->next() )); gp_XYZ size = bndBox.CornerMax() - bndBox.CornerMin(); - SMESH_TNodeXYZ p0( _edges[iFrom]->_2neibors->_nodes[0] ); - SMESH_TNodeXYZ p1( _edges[iFrom]->_2neibors->_nodes[1] ); + SMESH_TNodeXYZ p0( _edges[iFrom]->_2neibors->tgtNode(0) ); + SMESH_TNodeXYZ p1( _edges[iFrom]->_2neibors->tgtNode(1) ); const double lineTol = 1e-2 * ( p0 - p1 ).Modulus(); for ( int i = 0; i < 3 && !isLine; ++i ) isLine = ( size.Coord( i+1 ) <= lineTol ); @@ -3064,10 +3316,10 @@ void _SolidData::SortOnEdge( const TopoDS_Edge& E, // set _2neibors according to the new order for ( int i = iFrom; i < iTo-1; ++i ) - if ( _edges[i]->_2neibors->_nodes[1] != _edges[i+1]->_nodes.back() ) + if ( _edges[i]->_2neibors->tgtNode(1) != _edges[i+1]->_nodes.back() ) _edges[i]->_2neibors->reverse(); if ( u2edge.size() > 1 && - _edges[iTo-1]->_2neibors->_nodes[0] != _edges[iTo-2]->_nodes.back() ) + _edges[iTo-1]->_2neibors->tgtNode(0) != _edges[iTo-2]->_nodes.back() ) _edges[iTo-1]->_2neibors->reverse(); } @@ -3193,8 +3445,8 @@ bool _ViscousBuilder::smoothAnalyticEdge( _SolidData& data, { if ( F.IsNull() ) // 3D { - SMESH_TNodeXYZ p0( data._edges[iFrom]->_2neibors->_nodes[0]); - SMESH_TNodeXYZ p1( data._edges[iTo-1]->_2neibors->_nodes[1]); + SMESH_TNodeXYZ p0( data._edges[iFrom]->_2neibors->tgtNode(0)); + SMESH_TNodeXYZ p1( data._edges[iTo-1]->_2neibors->tgtNode(1)); for ( int i = iFrom; i < iTo; ++i ) { double r = len[i-iFrom] / len.back(); @@ -3207,10 +3459,12 @@ bool _ViscousBuilder::smoothAnalyticEdge( _SolidData& data, } else { - gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]); - gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]); - if ( data._edges[iFrom]->_2neibors->_nodes[0] == - data._edges[iTo-1]->_2neibors->_nodes[1] ) // closed edge + // gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->tgtNode(0)); + // gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->tgtNode(1)); + gp_XY uv0 = data._edges[iFrom]->_2neibors->_edges[0]->LastUV( F ); + gp_XY uv1 = data._edges[iTo-1]->_2neibors->_edges[1]->LastUV( F ); + if ( data._edges[iFrom]->_2neibors->tgtNode(0) == + data._edges[iTo-1]->_2neibors->tgtNode(1) ) // closed edge { int iPeriodic = helper.GetPeriodicIndex(); if ( iPeriodic == 1 || iPeriodic == 2 ) @@ -3247,8 +3501,8 @@ bool _ViscousBuilder::smoothAnalyticEdge( _SolidData& data, if ( F.IsNull() ) // 3D { - if ( data._edges[iFrom]->_2neibors->_nodes[0] == - data._edges[iTo-1]->_2neibors->_nodes[1] ) + if ( data._edges[iFrom]->_2neibors->tgtNode(0) == + data._edges[iTo-1]->_2neibors->tgtNode(1) ) return true; // closed EDGE - nothing to do return false; // TODO ??? @@ -3257,9 +3511,12 @@ bool _ViscousBuilder::smoothAnalyticEdge( _SolidData& data, { const gp_XY center( center3D.X(), center3D.Y() ); - gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]); - gp_XY uvM = helper.GetNodeUV( F, data._edges[iFrom]->_nodes.back()); - gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]); + gp_XY uv0 = data._edges[iFrom]->_2neibors->_edges[0]->LastUV( F ); + gp_XY uvM = data._edges[iFrom]->LastUV( F ); + gp_XY uv1 = data._edges[iTo-1]->_2neibors->_edges[1]->LastUV( F ); + // gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->tgtNode(0)); + // gp_XY uvM = helper.GetNodeUV( F, data._edges[iFrom]->_nodes.back()); + // gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->tgtNode(1)); gp_Vec2d vec0( center, uv0 ); gp_Vec2d vecM( center, uvM ); gp_Vec2d vec1( center, uv1 ); @@ -3323,7 +3580,7 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, const SMDS_MeshNode* tgt1 = edge->_nodes.back(); for ( int j = 0; j < 2; ++j ) // loop on _2NearEdges { - const SMDS_MeshNode* tgt2 = edge->_2neibors->_nodes[j]; + const SMDS_MeshNode* tgt2 = edge->_2neibors->tgtNode(j); pair< set< SMESH_TLink >::iterator, bool > link_isnew = extrudedLinks.insert( SMESH_TLink( tgt1, tgt2 )); if ( !link_isnew.second ) @@ -3535,8 +3792,8 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, edge1->SetNewLength( data._stepSize, helper ); if ( edge1->IsOnEdge() ) { - const SMDS_MeshNode * n1 = edge1->_2neibors->_edges[0]->_nodes[0]; - const SMDS_MeshNode * n2 = edge1->_2neibors->_edges[1]->_nodes[0]; + const SMDS_MeshNode * n1 = edge1->_2neibors->srcNode(0); + const SMDS_MeshNode * n2 = edge1->_2neibors->srcNode(1); edge1->SetDataByNeighbors( n1, n2, helper ); } @@ -4211,6 +4468,31 @@ gp_Ax1 _LayerEdge::LastSegment(double& segLen) const return segDir; } +//================================================================================ +/*! + * \brief Return the last position of the target node on a FACE. + * \param [in] F - the FACE this _LayerEdge is inflated along + * \return gp_XY - result UV + */ +//================================================================================ + +gp_XY _LayerEdge::LastUV( const TopoDS_Face& F ) const +{ + if ( F.IsSame( _sWOL )) // F is my FACE + return gp_XY( _pos.back().X(), _pos.back().Y() ); + + if ( _sWOL.IsNull() || _sWOL.ShapeType() != TopAbs_EDGE ) // wrong call + return gp_XY( 1e100, 1e100 ); + + // _sWOL is EDGE of F; _pos.back().X() is the last U on the EDGE + double f, l, u = _pos.back().X(); + Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( TopoDS::Edge(_sWOL), F, f,l); + if ( !C2d.IsNull() && f <= u && u <= l ) + return C2d->Value( u ).XY(); + + return gp_XY( 1e100, 1e100 ); +} + //================================================================================ /*! * \brief Test intersection of the last segment with a given triangle @@ -4294,9 +4576,9 @@ bool _LayerEdge::SmoothOnEdge(Handle(Geom_Surface)& surface, SMESH_TNodeXYZ oldPos( tgtNode ); double dist01, distNewOld; - SMESH_TNodeXYZ p0( _2neibors->_nodes[0]); - SMESH_TNodeXYZ p1( _2neibors->_nodes[1]); - dist01 = p0.Distance( _2neibors->_nodes[1] ); + SMESH_TNodeXYZ p0( _2neibors->tgtNode(0)); + SMESH_TNodeXYZ p1( _2neibors->tgtNode(1)); + dist01 = p0.Distance( _2neibors->tgtNode(1) ); gp_Pnt newPos = p0 * _2neibors->_wgt[0] + p1 * _2neibors->_wgt[1]; double lenDelta = 0; @@ -4436,17 +4718,23 @@ void _LayerEdge::SetNewLength( double len, SMESH_MesherHelper& helper ) double u = Precision::Infinite(); // to force projection w/o distance check helper.CheckNodeU( TopoDS::Edge( _sWOL ), n, u, 1e-10, /*force=*/true, distXYZ ); _pos.back().SetCoord( u, 0, 0 ); - SMDS_EdgePosition* pos = static_cast( n->GetPosition() ); - pos->SetUParameter( u ); + if ( _nodes.size() > 1 ) + { + SMDS_EdgePosition* pos = static_cast( n->GetPosition() ); + pos->SetUParameter( u ); + } } else // TopAbs_FACE { gp_XY uv( Precision::Infinite(), 0 ); helper.CheckNodeUV( TopoDS::Face( _sWOL ), n, uv, 1e-10, /*force=*/true, distXYZ ); _pos.back().SetCoord( uv.X(), uv.Y(), 0 ); - SMDS_FacePosition* pos = static_cast( n->GetPosition() ); - pos->SetUParameter( uv.X() ); - pos->SetVParameter( uv.Y() ); + if ( _nodes.size() > 1 ) + { + SMDS_FacePosition* pos = static_cast( n->GetPosition() ); + pos->SetUParameter( uv.X() ); + pos->SetVParameter( uv.Y() ); + } } n->setXYZ( distXYZ[1], distXYZ[2], distXYZ[3]); } @@ -4519,10 +4807,15 @@ bool _ViscousBuilder::refine(_SolidData& data) TNode2Edge* n2eMap = 0; TNode2Edge::iterator n2e; + // Create intermediate nodes on each _LayerEdge + for ( size_t i = 0; i < data._edges.size(); ++i ) { _LayerEdge& edge = *data._edges[i]; + if ( edge._nodes.size() < 2 ) + continue; // on _noShrinkShapes + // get accumulated length of segments vector< double > segLen( edge._pos.size() ); segLen[0] = 0.0; @@ -4561,21 +4854,22 @@ bool _ViscousBuilder::refine(_SolidData& data) n2eMap = ( s2ne == data._s2neMap.end() ) ? 0 : n2eMap = s2ne->second; prevBaseId = baseShapeId; } + _LayerEdge* edgeOnSameNode = 0; 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(); + edgeOnSameNode = n2e->second; + const gp_XYZ& otherTgtPos = edgeOnSameNode->_pos.back(); + SMDS_PositionPtr lastPos = tgtNode->GetPosition(); if ( isOnEdge ) { SMDS_EdgePosition* epos = static_cast( lastPos ); - epos->SetUParameter( foundPos.X() ); + epos->SetUParameter( otherTgtPos.X() ); } else { SMDS_FacePosition* fpos = static_cast( lastPos ); - fpos->SetUParameter( foundPos.X() ); - fpos->SetVParameter( foundPos.Y() ); + fpos->SetUParameter( otherTgtPos.X() ); + fpos->SetVParameter( otherTgtPos.Y() ); } } // calculate height of the first layer @@ -4607,7 +4901,7 @@ bool _ViscousBuilder::refine(_SolidData& data) double r = ( segLen[iSeg] - hSum ) / ( segLen[iSeg] - segLen[iPrevSeg] ); gp_Pnt pos = r * edge._pos[iPrevSeg] + (1-r) * edge._pos[iSeg]; - SMDS_MeshNode*& node = const_cast< SMDS_MeshNode*& >(edge._nodes[ iStep ]); + SMDS_MeshNode*& node = const_cast< SMDS_MeshNode*& >( edge._nodes[ iStep ]); if ( !edge._sWOL.IsNull() ) { // compute XYZ by parameters @@ -4665,8 +4959,20 @@ bool _ViscousBuilder::refine(_SolidData& data) } node->setXYZ( pos.X(), pos.Y(), pos.Z() ); } + } // loop on edge._nodes + + if ( !edge._sWOL.IsNull() ) // prepare for shrink() + { + if ( isOnEdge ) + edge._pos.back().SetCoord( u, 0,0); + else + edge._pos.back().SetCoord( uv.X(), uv.Y() ,0); + + if ( edgeOnSameNode ) + edgeOnSameNode->_pos.back() = edge._pos.back(); } - } + + } // loop on data._edges to create nodes if ( !getMeshDS()->IsEmbeddedMode() ) // Log node movement @@ -4677,50 +4983,151 @@ bool _ViscousBuilder::refine(_SolidData& data) getMeshDS()->MoveNode( p._node, p.X(), p.Y(), p.Z() ); } - // TODO: make quadratic prisms and polyhedrons(?) + // Create volumes helper.SetElementsOnShape(true); + vector< vector* > nnVec; + set< vector* > nnSet; + set< int > degenEdgeInd; + vector degenVols; + TopExp_Explorer exp( data._solid, TopAbs_FACE ); for ( ; exp.More(); exp.Next() ) { if ( data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() ))) continue; - SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() ); + SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() ); SMDS_ElemIteratorPtr fIt = fSubM->GetElements(); - vector< vector* > nnVec; while ( fIt->more() ) { const SMDS_MeshElement* face = fIt->next(); - int nbNodes = face->NbCornerNodes(); + const int nbNodes = face->NbCornerNodes(); nnVec.resize( nbNodes ); - SMDS_ElemIteratorPtr nIt = face->nodesIterator(); + nnSet.clear(); + degenEdgeInd.clear(); + int nbZ = 0; + SMDS_NodeIteratorPtr nIt = face->nodeIterator(); for ( int iN = 0; iN < nbNodes; ++iN ) { - const SMDS_MeshNode* n = static_cast( nIt->next() ); + const SMDS_MeshNode* n = nIt->next(); nnVec[ iN ] = & data._n2eMap[ n ]->_nodes; + if ( nnVec[ iN ]->size() < 2 ) + degenEdgeInd.insert( iN ); + else + nbZ = nnVec[ iN ]->size(); + + if ( helper.HasDegeneratedEdges() ) + nnSet.insert( nnVec[ iN ]); } + if ( nbZ == 0 ) + continue; + if ( 0 < nnSet.size() && nnSet.size() < 3 ) + continue; - int nbZ = nnVec[0]->size(); switch ( nbNodes ) { case 3: - for ( int iZ = 1; iZ < nbZ; ++iZ ) - helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1], - (*nnVec[0])[iZ], (*nnVec[1])[iZ], (*nnVec[2])[iZ]); + switch ( degenEdgeInd.size() ) + { + case 0: // PENTA + { + for ( int iZ = 1; iZ < nbZ; ++iZ ) + helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1], + (*nnVec[0])[iZ], (*nnVec[1])[iZ], (*nnVec[2])[iZ]); + break; + } + case 1: // PYRAM + { + int i2 = *degenEdgeInd.begin(); + int i0 = helper.WrapIndex( i2 - 1, nbNodes ); + int i1 = helper.WrapIndex( i2 + 1, nbNodes ); + for ( int iZ = 1; iZ < nbZ; ++iZ ) + helper.AddVolume( (*nnVec[i0])[iZ-1], (*nnVec[i1])[iZ-1], + (*nnVec[i1])[iZ], (*nnVec[i0])[iZ], (*nnVec[i2])[0]); + break; + } + case 2: // TETRA + { + int i3 = !degenEdgeInd.count(0) ? 0 : !degenEdgeInd.count(1) ? 1 : 2; + for ( int iZ = 1; iZ < nbZ; ++iZ ) + helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1], + (*nnVec[i3])[iZ]); + break; + } + } break; + case 4: - for ( int iZ = 1; iZ < nbZ; ++iZ ) - helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], - (*nnVec[2])[iZ-1], (*nnVec[3])[iZ-1], - (*nnVec[0])[iZ], (*nnVec[1])[iZ], - (*nnVec[2])[iZ], (*nnVec[3])[iZ]); + switch ( degenEdgeInd.size() ) + { + case 0: // HEX + { + for ( int iZ = 1; iZ < nbZ; ++iZ ) + helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], + (*nnVec[2])[iZ-1], (*nnVec[3])[iZ-1], + (*nnVec[0])[iZ], (*nnVec[1])[iZ], + (*nnVec[2])[iZ], (*nnVec[3])[iZ]); + break; + } + case 2: // PENTA? + { + int i2 = *degenEdgeInd.begin(); + int i3 = *degenEdgeInd.rbegin(); + bool ok = ( i3 - i2 == 1 ); + if ( i2 == 0 && i3 == 3 ) { i2 = 3; i3 = 0; ok = true; } + int i0 = helper.WrapIndex( i3 + 1, nbNodes ); + int i1 = helper.WrapIndex( i0 + 1, nbNodes ); + for ( int iZ = 1; iZ < nbZ; ++iZ ) + { + const SMDS_MeshElement* vol = + helper.AddVolume( (*nnVec[i3])[0], (*nnVec[i0])[iZ], (*nnVec[i0])[iZ-1], + (*nnVec[i2])[0], (*nnVec[i1])[iZ], (*nnVec[i1])[iZ-1]); + if ( !ok && vol ) + degenVols.push_back( vol ); + } + break; + } + case 3: // degen HEX + { + const SMDS_MeshNode* nn[8]; + for ( int iZ = 1; iZ < nbZ; ++iZ ) + { + const SMDS_MeshElement* vol = + helper.AddVolume( nnVec[0]->size() > 1 ? (*nnVec[0])[iZ-1] : (*nnVec[0])[0], + nnVec[1]->size() > 1 ? (*nnVec[1])[iZ-1] : (*nnVec[1])[0], + nnVec[2]->size() > 1 ? (*nnVec[2])[iZ-1] : (*nnVec[2])[0], + nnVec[3]->size() > 1 ? (*nnVec[3])[iZ-1] : (*nnVec[3])[0], + nnVec[0]->size() > 1 ? (*nnVec[0])[iZ] : (*nnVec[0])[0], + nnVec[1]->size() > 1 ? (*nnVec[1])[iZ] : (*nnVec[1])[0], + nnVec[2]->size() > 1 ? (*nnVec[2])[iZ] : (*nnVec[2])[0], + nnVec[3]->size() > 1 ? (*nnVec[3])[iZ] : (*nnVec[3])[0]); + degenVols.push_back( vol ); + } + } break; + } + break; + default: return error("Not supported type of element", data._index); - } + + } // switch ( nbNodes ) + } // while ( fIt->more() ) + } // loop on FACEs + + if ( !degenVols.empty() ) + { + SMESH_ComputeErrorPtr& err = _mesh->GetSubMesh( data._solid )->GetComputeError(); + if ( !err || err->IsOK() ) + { + err.reset( new SMESH_ComputeError( COMPERR_WARNING, + "Degenerated volumes created" )); + err->myBadElements.insert( err->myBadElements.end(), + degenVols.begin(),degenVols.end() ); } } + return true; } @@ -4766,20 +5173,19 @@ bool _ViscousBuilder::shrink() // EDGE's to shrink map< TGeomID, _Shrinker1D > e2shrMap; + vector< _LayerEdge* > lEdges; // loop on FACES to srink mesh on map< TGeomID, _SolidData* >::iterator f2sd = f2sdMap.begin(); for ( ; f2sd != f2sdMap.end(); ++f2sd ) { - _SolidData& data = *f2sd->second; - TNode2Edge& n2eMap = data._n2eMap; - const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first )); - - Handle(Geom_Surface) surface = BRep_Tool::Surface(F); - + _SolidData& data = *f2sd->second; + const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first )); SMESH_subMesh* sm = _mesh->GetSubMesh( F ); SMESHDS_SubMesh* smDS = sm->GetSubMeshDS(); + Handle(Geom_Surface) surface = BRep_Tool::Surface(F); + helper.SetSubShape(F); // =========================== @@ -4814,23 +5220,30 @@ bool _ViscousBuilder::shrink() } // Find _LayerEdge's inflated along F - vector< _LayerEdge* > lEdges; + lEdges.clear(); { - SMESH_subMeshIteratorPtr subIt = - sm->getDependsOnIterator(/*includeSelf=*/false, /*complexShapeFirst=*/false); + set< TGeomID > subIDs; + SMESH_subMeshIteratorPtr subIt = sm->getDependsOnIterator(/*includeSelf=*/false); while ( subIt->more() ) + subIDs.insert( subIt->next()->GetId() ); + + int iBeg, iEnd = 0; + for ( int iS = 0; iS < data._endEdgeOnShape.size() && !subIDs.empty(); ++iS ) { - SMESH_subMesh* sub = subIt->next(); - SMESHDS_SubMesh* subDS = sub->GetSubMeshDS(); - if ( subDS->NbNodes() == 0 || !n2eMap.count( subDS->GetNodes()->next() )) - continue; - SMDS_NodeIteratorPtr nIt = subDS->GetNodes(); - while ( nIt->more() ) - { - _LayerEdge* edge = n2eMap[ nIt->next() ]; - lEdges.push_back( edge ); - prepareEdgeToShrink( *edge, F, helper, smDS ); - } + iBeg = iEnd; + iEnd = data._endEdgeOnShape[ iS ]; + TGeomID shapeID = data._edges[ iBeg ]->_nodes[0]->getshapeId(); + set< TGeomID >::iterator idIt = subIDs.find( shapeID ); + if ( idIt == subIDs.end() || + data._edges[ iBeg ]->_sWOL.IsNull() ) continue; + subIDs.erase( idIt ); + + if ( !data._noShrinkShapes.count( shapeID )) + for ( ; iBeg < iEnd; ++iBeg ) + { + lEdges.push_back( data._edges[ iBeg ] ); + prepareEdgeToShrink( *data._edges[ iBeg ], F, helper, smDS ); + } } } @@ -5087,12 +5500,10 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, const SMDS_MeshNode* srcNode = edge._nodes[0]; const SMDS_MeshNode* tgtNode = edge._nodes.back(); - edge._pos.clear(); - if ( edge._sWOL.ShapeType() == TopAbs_FACE ) { - gp_XY srcUV = helper.GetNodeUV( F, srcNode ); - gp_XY tgtUV = helper.GetNodeUV( F, tgtNode ); + gp_XY srcUV( edge._pos[0].X(), edge._pos[0].Y() );//helper.GetNodeUV( F, srcNode ); + gp_XY tgtUV = edge.LastUV( F ); //helper.GetNodeUV( F, tgtNode ); gp_Vec2d uvDir( srcUV, tgtUV ); double uvLen = uvDir.Magnitude(); uvDir /= uvLen; @@ -5109,7 +5520,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, } else // _sWOL is TopAbs_EDGE { - TopoDS_Edge E = TopoDS::Edge( edge._sWOL); + const TopoDS_Edge& E = TopoDS::Edge( edge._sWOL ); SMESHDS_SubMesh* edgeSM = getMeshDS()->MeshElements( E ); if ( !edgeSM || edgeSM->NbElements() == 0 ) return error(SMESH_Comment("Not meshed EDGE ") << getMeshDS()->ShapeToIndex( E )); @@ -5128,7 +5539,9 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, double uSrc = helper.GetNodeU( E, srcNode, n2 ); double uTgt = helper.GetNodeU( E, tgtNode, srcNode ); - double u2 = helper.GetNodeU( E, n2, srcNode ); + double u2 = helper.GetNodeU( E, n2, srcNode ); + + edge._pos.clear(); if ( fabs( uSrc-uTgt ) < 0.99 * fabs( uSrc-u2 )) { @@ -5143,13 +5556,56 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, edge._simplices.resize( 1 ); edge._simplices[0]._nPrev = n2; - // set UV of source node to target node + // set U of source node to the target node SMDS_EdgePosition* pos = static_cast( tgtNode->GetPosition() ); pos->SetUParameter( uSrc ); } return true; } +//================================================================================ +/*! + * \brief Restore position of a sole node of a _LayerEdge based on _noShrinkShapes + */ +//================================================================================ + +void _ViscousBuilder::restoreNoShrink( _LayerEdge& edge ) const +{ + if ( edge._nodes.size() == 1 ) + { + edge._pos.clear(); + edge._len = 0; + + const SMDS_MeshNode* srcNode = edge._nodes[0]; + TopoDS_Shape S = SMESH_MesherHelper::GetSubShapeByNode( srcNode, getMeshDS() ); + if ( S.IsNull() ) return; + + gp_Pnt p; + + switch ( S.ShapeType() ) + { + case TopAbs_EDGE: + { + double f,l; + TopLoc_Location loc; + Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( S ), loc, f, l ); + if ( curve.IsNull() ) return; + SMDS_EdgePosition* ePos = static_cast( srcNode->GetPosition() ); + p = curve->Value( ePos->GetUParameter() ); + break; + } + case TopAbs_VERTEX: + { + p = BRep_Tool::Pnt( TopoDS::Vertex( S )); + break; + } + default: return; + } + getMeshDS()->MoveNode( srcNode, p.X(), p.Y(), p.Z() ); + dumpMove( srcNode ); + } +} + //================================================================================ /*! * \brief Try to fix triangles with high aspect ratio by swaping diagonals @@ -5366,11 +5822,11 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface, } else // _sWOL is TopAbs_EDGE { - TopoDS_Edge E = TopoDS::Edge( _sWOL ); - const SMDS_MeshNode* n2 = _simplices[0]._nPrev; + const 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 u2 = helper.GetNodeU( E, n2, tgtNode ); const double uSrc = _pos[0].Coord( U_SRC ); const double lenTgt = _pos[0].Coord( LEN_TGT ); @@ -5777,6 +6233,8 @@ bool _ViscousBuilder::addBoundaryElements() { SMESH_MesherHelper helper( *_mesh ); + vector< const SMDS_MeshNode* > faceNodes; + for ( size_t i = 0; i < _sdVec.size(); ++i ) { _SolidData& data = _sdVec[i]; @@ -5785,6 +6243,8 @@ bool _ViscousBuilder::addBoundaryElements() for ( int iE = 1; iE <= geomEdges.Extent(); ++iE ) { const TopoDS_Edge& E = TopoDS::Edge( geomEdges(iE)); + if ( data._noShrinkShapes.count( getMeshDS()->ShapeToIndex( E ))) + continue; // Get _LayerEdge's based on E @@ -5817,8 +6277,13 @@ bool _ViscousBuilder::addBoundaryElements() if ( nbSharedPyram > 1 ) continue; // not free border of the pyramid - if ( getMeshDS()->FindFace( ledges[0]->_nodes[0], ledges[0]->_nodes[1], - ledges[1]->_nodes[0], ledges[1]->_nodes[1])) + faceNodes.clear(); + faceNodes.push_back( ledges[0]->_nodes[0] ); + faceNodes.push_back( ledges[1]->_nodes[0] ); + if ( ledges[0]->_nodes.size() > 1 ) faceNodes.push_back( ledges[0]->_nodes[1] ); + if ( ledges[1]->_nodes.size() > 1 ) faceNodes.push_back( ledges[1]->_nodes[1] ); + + if ( getMeshDS()->FindElement( faceNodes, SMDSAbs_Face, /*noMedium=*/true)) continue; // faces already created } for ( ++u2n; u2n != u2nodes.end(); ++u2n ) @@ -5868,12 +6333,33 @@ bool _ViscousBuilder::addBoundaryElements() { vector< const SMDS_MeshNode*>& nn1 = ledges[j-dj1]->_nodes; vector< const SMDS_MeshNode*>& nn2 = ledges[j-dj2]->_nodes; - if ( isOnFace ) - for ( size_t z = 1; z < nn1.size(); ++z ) - sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] )); + if ( nn1.size() == nn2.size() ) + { + if ( isOnFace ) + for ( size_t z = 1; z < nn1.size(); ++z ) + sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] )); + else + for ( size_t z = 1; z < nn1.size(); ++z ) + sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z] )); + } + else if ( nn1.size() == 1 ) + { + if ( isOnFace ) + for ( size_t z = 1; z < nn2.size(); ++z ) + sm->AddElement( getMeshDS()->AddFace( nn1[0], nn2[z-1], nn2[z] )); + else + for ( size_t z = 1; z < nn2.size(); ++z ) + sm->AddElement( new SMDS_FaceOfNodes( nn1[0], nn2[z-1], nn2[z] )); + } else - for ( size_t z = 1; z < nn1.size(); ++z ) - sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z])); + { + if ( isOnFace ) + for ( size_t z = 1; z < nn1.size(); ++z ) + sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[0], nn1[z] )); + else + for ( size_t z = 1; z < nn1.size(); ++z ) + sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[0], nn2[z] )); + } } // Make edges @@ -5883,7 +6369,7 @@ bool _ViscousBuilder::addBoundaryElements() if ( !edge->_sWOL.IsNull() && edge->_sWOL.ShapeType() == TopAbs_EDGE ) { vector< const SMDS_MeshNode*>& nn = edge->_nodes; - if ( nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() ) + if ( nn.size() < 2 || nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() ) continue; helper.SetSubShape( edge->_sWOL ); helper.SetElementsOnShape( true );