X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_QuadFromMedialAxis_1D2D.cxx;h=534860754bb28827ca79fbae0d6323f897cfc155;hp=ea8fdc9361b7e15f39a74b206cb050c80e03def4;hb=38f832b912a159ca9d43b5bd031344fd7108a6d8;hpb=0e017d4c87d670c01e72c0b4858ffcb9e47aa9ba diff --git a/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx index ea8fdc936..534860754 100644 --- a/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx +++ b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx @@ -134,13 +134,14 @@ namespace */ struct SinuousFace { - FaceQuadStruct::Ptr _quad; - vector< TopoDS_Edge > _edges; - vector< TopoDS_Edge > _sinuSide[2], _shortSide[2]; - vector< TopoDS_Edge > _sinuEdges; - int _nbWires; - list< int > _nbEdgesInWire; - TMergeMap _nodesToMerge; + FaceQuadStruct::Ptr _quad; + vector< TopoDS_Edge > _edges; + vector< TopoDS_Edge > _sinuSide[2], _shortSide[2]; + vector< TopoDS_Edge > _sinuEdges; + vector< Handle(Geom_Curve) > _sinuCurves; + int _nbWires; + list< int > _nbEdgesInWire; + TMergeMap _nodesToMerge; SinuousFace( const TopoDS_Face& f ): _quad( new FaceQuadStruct ) { @@ -830,11 +831,11 @@ namespace */ //================================================================================ - bool findVertex( NodePoint& theNodePnt, - const vector& theSinuEdges, - size_t theEdgeIndPrev, - size_t theEdgeIndNext, - SMESHDS_Mesh* theMeshDS) + bool findVertexAndNode( NodePoint& theNodePnt, + const vector& theSinuEdges, + SMESHDS_Mesh* theMeshDS = 0, + size_t theEdgeIndPrev = 0, + size_t theEdgeIndNext = 0) { if ( theNodePnt._edgeInd >= theSinuEdges.size() ) return false; @@ -851,7 +852,7 @@ namespace else if ( theEdgeIndPrev != theEdgeIndNext ) TopExp::CommonVertex( theSinuEdges[theEdgeIndPrev], theSinuEdges[theEdgeIndNext], V ); - if ( !V.IsNull() ) + if ( !V.IsNull() && theMeshDS ) { theNodePnt._node = SMESH_Algo::VertexNode( V, theMeshDS ); if ( !theNodePnt._node ) @@ -860,9 +861,8 @@ namespace theNodePnt._node = theMeshDS->AddNode( p.X(), p.Y(), p.Z() ); theMeshDS->SetNodeOnVertex( theNodePnt._node, V ); } - return true; } - return false; + return !V.IsNull(); } //================================================================================ @@ -881,41 +881,44 @@ namespace //================================================================================ bool projectVertices( SMESH_MesherHelper& theHelper, - //const double theMinSegLen, const SMESH_MAT2d::MedialAxis& theMA, const vector< SMESH_MAT2d::BranchPoint >& theDivPoints, const vector< std::size_t > & theEdgeIDs1, const vector< std::size_t > & theEdgeIDs2, - const vector& theSinuEdges, - const vector< Handle(Geom_Curve) >& theCurves, const vector< bool >& theIsEdgeComputed, map< double, pair< NodePoint, NodePoint > > & thePointsOnE, - TMergeMap& theNodes2Merge) + SinuousFace& theSinuFace) { - if ( theDivPoints.empty() ) - return true; - SMESHDS_Mesh* meshDS = theHelper.GetMeshDS(); + const vector& theSinuEdges = theSinuFace._sinuEdges; + const vector< Handle(Geom_Curve) >& theCurves = theSinuFace._sinuCurves; double uMA; SMESH_MAT2d::BoundaryPoint bp[2]; const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0); - // fill a map holding NodePoint's of ends of theSinuEdges - map< double, pair< NodePoint, NodePoint > > extremaNP; - map< double, pair< NodePoint, NodePoint > >::iterator u2NP0, u2NP1; + // add to thePointsOnE NodePoint's of ends of theSinuEdges if ( !branch.getBoundaryPoints( 0., bp[0], bp[1] ) || !theMA.getBoundary().moveToClosestEdgeEnd( bp[0] ) || !theMA.getBoundary().moveToClosestEdgeEnd( bp[1] )) return false; - u2NP0 = extremaNP.insert - ( make_pair( 0., make_pair( NodePoint( bp[0]), NodePoint( bp[1])))).first; + NodePoint np0( bp[0]), np1( bp[1] ); + findVertexAndNode( np0, theSinuEdges, meshDS ); + findVertexAndNode( np1, theSinuEdges, meshDS ); + thePointsOnE.insert( make_pair( -0.1, make_pair( np0, np1 ))); + if ( !branch.getBoundaryPoints( 1., bp[0], bp[1] ) || !theMA.getBoundary().moveToClosestEdgeEnd( bp[0] ) || !theMA.getBoundary().moveToClosestEdgeEnd( bp[1] )) return false; - u2NP1 = extremaNP.insert - ( make_pair( 1., make_pair( NodePoint( bp[0]), NodePoint( bp[1])))).first; + np0 = bp[0]; np1 = bp[1]; + findVertexAndNode( np0, theSinuEdges, meshDS ); + findVertexAndNode( np1, theSinuEdges, meshDS ); + thePointsOnE.insert( make_pair( 1.1, make_pair( np0, np1))); // project theDivPoints + + if ( theDivPoints.empty() ) + return true; + for ( size_t i = 0; i < theDivPoints.size(); ++i ) { if ( !branch.getParameter( theDivPoints[i], uMA )) @@ -928,8 +931,8 @@ namespace NodePoint( bp[1] ) }; bool isVertex[2] = { - findVertex( np[0], theSinuEdges, theEdgeIDs1[i], theEdgeIDs1[i+1], meshDS ), - findVertex( np[1], theSinuEdges, theEdgeIDs2[i], theEdgeIDs2[i+1], meshDS ) + findVertexAndNode( np[0], theSinuEdges, meshDS, theEdgeIDs1[i], theEdgeIDs1[i+1] ), + findVertexAndNode( np[1], theSinuEdges, meshDS, theEdgeIDs2[i], theEdgeIDs2[i+1] ) }; map< double, pair< NodePoint, NodePoint > >::iterator u2NP = @@ -957,10 +960,10 @@ namespace bool isShortPrev[2], isShortNext[2]; map< double, pair< NodePoint, NodePoint > >::iterator u2NPPrev = u2NP, u2NPNext = u2NP; --u2NPPrev; ++u2NPNext; - bool hasPrev = ( u2NP != thePointsOnE.begin() ); - bool hasNext = ( u2NPNext != thePointsOnE.end() ); - if ( !hasPrev ) u2NPPrev = u2NP0; - if ( !hasNext ) u2NPNext = u2NP1; + // bool hasPrev = ( u2NP != thePointsOnE.begin() ); + // bool hasNext = ( u2NPNext != thePointsOnE.end() ); + // if ( !hasPrev ) u2NPPrev = u2NP0; + // if ( !hasNext ) u2NPNext = u2NP1; for ( int iS = 0; iS < 2; ++iS ) // side with Vertex and side with Nodes { NodePoint np = get( u2NP->second, iS ); @@ -1003,19 +1006,134 @@ namespace u2NPClose = isShortPrev[ iNode ] ? u2NPPrev : u2NPNext; NodePoint& npProj = get( u2NP->second, iNode ); // NP of VERTEX projection NodePoint& npCloseN = get( u2NPClose->second, iNode ); // NP close to npProj - // npProj._edgeInd = npCloseN._edgeInd; + npProj = npCloseN; + npProj._node = 0; + //npProj._edgeInd = npCloseN._edgeInd; // npProj._u = npCloseN._u + 1e-3 * Abs( get( u2NPPrev->second, iNode )._u - // get( u2NPNext->second, iNode )._u ); - gp_Pnt p = npProj.Point( theCurves ); - npProj._node = meshDS->AddNode( p.X(), p.Y(), p.Z() ); - meshDS->SetNodeOnEdge( npProj._node, theSinuEdges[ npProj._edgeInd ], npProj._u ); + // gp_Pnt p = npProj.Point( theCurves ); + // npProj._node = meshDS->AddNode( p.X(), p.Y(), p.Z() ); + // meshDS->SetNodeOnEdge( npProj._node, theSinuEdges[ npProj._edgeInd ], npProj._u ); - theNodes2Merge[ npCloseN._node ].push_back( npProj._node ); + //theNodes2Merge[ npCloseN._node ].push_back( npProj._node ); } } return true; } + //================================================================================ + /*! + * \brief Move coincident nodes to make node params on EDGE unique + * \param [in] theHelper - the helper + * \param [in] thePointsOnE - nodes on two opposite river sides + * \param [in] theSinuFace - the sinuous FACE + * \param [out] theNodes2Merge - the map of nodes to merge + */ + //================================================================================ + + void separateNodes( SMESH_MesherHelper& theHelper, + map< double, pair< NodePoint, NodePoint > > & thePointsOnE, + SinuousFace& theSinuFace ) + { + if ( thePointsOnE.size() < 2 ) + return; + + SMESHDS_Mesh* meshDS = theHelper.GetMeshDS(); + const vector& theSinuEdges = theSinuFace._sinuEdges; + + typedef map< double, pair< NodePoint, NodePoint > >::iterator TIterator; + + for ( int iSide = 0; iSide < 2; ++iSide ) + { + TIterator u2NP0, u2NP1, u2NP = thePointsOnE.begin(); + while ( u2NP != thePointsOnE.end() ) + { + while ( u2NP != thePointsOnE.end() && + get( u2NP->second, iSide )._node ) + ++u2NP; // skip NP with an existing node (VERTEXes must be meshed) + if ( u2NP == thePointsOnE.end() ) + break; + + // find a range of not meshed NP on one EDGE + u2NP0 = u2NP; + if ( !findVertexAndNode( get( u2NP0->second, iSide ), theSinuEdges )) + --u2NP0; + int iCurEdge = get( u2NP->second, iSide )._edgeInd; + int nbNP = 1; + while ( get( u2NP->second, iSide )._edgeInd == iCurEdge && + get( u2NP->second, iSide )._node == 0 ) + ++u2NP, ++nbNP; + u2NP1 = u2NP; // end of not meshed NP on iCurEdge + + // fix parameters of extremity NP of the range + NodePoint* np0 = & get( u2NP0->second, iSide ); + NodePoint* np1 = & get( u2NP1->second, iSide ); + const TopoDS_Edge& edge = TopoDS::Edge( theSinuFace._sinuEdges[ iCurEdge ]); + if ( np0->_node && np0->_edgeInd != iCurEdge ) + { + np0->_u = theHelper.GetNodeU( edge, np0->_node ); + np0->_edgeInd = iCurEdge; + } + if ( np1->_node && np1->_edgeInd != iCurEdge ) + { + np1->_u = theHelper.GetNodeU( edge, np1->_node ); + np1->_edgeInd = iCurEdge; + } + + // find coincident NPs + double f,l; + BRep_Tool::Range( edge, f,l ); + double tol = 1e-2* (l-f) / nbNP; + TIterator u2NPEq = thePointsOnE.end(); + u2NP = u2NP0; + for ( ++u2NP; u2NP0 != u2NP1; ++u2NP, ++u2NP0 ) + { + np0 = & get( u2NP0->second, iSide ); + np1 = & get( u2NP->second, iSide ); + bool coincides = ( Abs( np0->_u - np1->_u ) < tol ); + if ( coincides && u2NPEq == thePointsOnE.end() ) + u2NPEq = u2NP0; + + if (( u2NPEq != thePointsOnE.end() ) && + ( u2NP == u2NP1 || !coincides )) + { + if ( !get( u2NPEq->second, iSide )._node ) + --u2NPEq; + if ( coincides && !get( u2NP->second, iSide )._node && u2NP0 != u2NP1 ) + ++u2NP; + + // distribute nodes between u2NPEq and u2NP + size_t nbSeg = std::distance( u2NPEq, u2NP ); + double du = 1. / nbSeg * ( get( u2NP->second, iSide )._u - + get( u2NPEq->second, iSide )._u ); + double u = get( u2NPEq->second, iSide )._u + du; + + const SMDS_MeshNode* closeNode = + get(( coincides ? u2NP : u2NPEq )->second, iSide )._node; + list< const SMDS_MeshNode* >& eqNodes = theSinuFace._nodesToMerge[ closeNode ]; + + for ( ++u2NPEq; u2NPEq != u2NP; ++u2NPEq, u += du ) + { + np0 = & get( u2NPEq->second, iSide ); + np0->_u = u; + gp_Pnt p = np0->Point( theSinuFace._sinuCurves ); + np0->_node = meshDS->AddNode( p.X(), p.Y(), p.Z() ); + meshDS->SetNodeOnEdge( np0->_node, theSinuEdges[ np0->_edgeInd ], np0->_u ); + if ( !closeNode ) + eqNodes = theSinuFace._nodesToMerge[ closeNode = np0->_node ]; + else + eqNodes.push_back( np0->_node ); + } + } + } + u2NP = u2NP1; + while ( get( u2NP->second, iSide )._edgeInd != iCurEdge ) + --u2NP; + u2NP++; + } + } + } + //================================================================================ /*! * \brief Divide the sinuous EDGEs by projecting the division point of Medial @@ -1048,11 +1166,12 @@ namespace SMESHDS_Mesh* meshDS = theHelper.GetMeshDS(); double f,l; + // get data of sinuous EDGEs and remove unnecessary nodes const vector< TopoDS_Edge >& theSinuEdges = theSinuFace._sinuEdges; - vector< Handle(Geom_Curve) > curves ( theSinuEdges.size() ); + vector< Handle(Geom_Curve) >& curves = theSinuFace._sinuCurves; vector< int > edgeIDs( theSinuEdges.size() ); vector< bool > isComputed( theSinuEdges.size() ); - //bool hasComputed = false; + curves.resize( theSinuEdges.size(), 0 ); for ( size_t i = 0; i < theSinuEdges.size(); ++i ) { curves[i] = BRep_Tool::Curve( theSinuEdges[i], f,l ); @@ -1172,10 +1291,12 @@ namespace ++iEdgePair; } - if ( !projectVertices( theHelper, theMA, divPoints, edgeIDs1, edgeIDs2, theSinuEdges, - curves, isComputed, pointsOnE, theSinuFace._nodesToMerge )) + if ( !projectVertices( theHelper, theMA, divPoints, edgeIDs1, edgeIDs2, + isComputed, pointsOnE, theSinuFace )) return false; + separateNodes( theHelper, pointsOnE, theSinuFace ); + // create nodes TMAPar2NPoints::iterator u2np = pointsOnE.begin(); for ( ; u2np != pointsOnE.end(); ++u2np ) @@ -1353,6 +1474,7 @@ namespace TMergeMap::iterator n2nn = theSinuFace._nodesToMerge.begin(); for ( ; n2nn != theSinuFace._nodesToMerge.end(); ++n2nn ) { + if ( !n2nn->first ) continue; nodesGroups.push_back( list< const SMDS_MeshNode* >() ); list< const SMDS_MeshNode* > & group = nodesGroups.back();