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=4ac71c0b818a8b09e31b0d0d8d7a618978c24738;hb=38f832b912a159ca9d43b5bd031344fd7108a6d8;hpb=23d90107acec5e54ded86d9f309fe5cb42720b78 diff --git a/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx index 4ac71c0b8..534860754 100644 --- a/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx +++ b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx @@ -70,7 +70,7 @@ public: void SetSegmentLength( double len ) { _value[ BEG_LENGTH_IND ] = len; - _value[ PRECISION_IND ] = 1e-7; + _value[ PRECISION_IND ] = 1e-7; _hypType = LOCAL_LENGTH; } }; @@ -95,7 +95,7 @@ StdMeshers_QuadFromMedialAxis_1D2D::StdMeshers_QuadFromMedialAxis_1D2D(int _neededLowerHyps[ 1 ] = true; // suppress warning on hiding a global 1D algo _neededLowerHyps[ 2 ] = true; // suppress warning on hiding a global 2D algo _compatibleHypothesis.clear(); - //_compatibleHypothesis.push_back("ViscousLayers2D"); + _compatibleHypothesis.push_back("ViscousLayers2D"); } //================================================================================ @@ -120,6 +120,7 @@ bool StdMeshers_QuadFromMedialAxis_1D2D::CheckHypothesis(SMESH_Mesh& aMe const TopoDS_Shape& aShape, Hypothesis_Status& aStatus) { + aStatus = HYP_OK; return true; // does not require hypothesis } @@ -133,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 ) { @@ -591,27 +593,45 @@ namespace //================================================================================ TopoDS_Edge makeEdgeFromMA( SMESH_MesherHelper& theHelper, - const SMESH_MAT2d::MedialAxis& theMA ) + const SMESH_MAT2d::MedialAxis& theMA, + const double theMinSegLen) { - if ( theMA.getBranches().size() != 1 ) + if ( theMA.nbBranches() != 1 ) return TopoDS_Edge(); vector< gp_XY > uv; - theMA.getPoints( theMA.getBranches()[0], uv ); + theMA.getPoints( theMA.getBranch(0), uv ); if ( uv.size() < 2 ) return TopoDS_Edge(); TopoDS_Face face = TopoDS::Face( theHelper.GetSubShape() ); Handle(Geom_Surface) surface = BRep_Tool::Surface( face ); + vector< gp_Pnt > pnt; + pnt.reserve( uv.size() * 2 ); + pnt.push_back( surface->Value( uv[0].X(), uv[0].Y() )); + for ( size_t i = 1; i < uv.size(); ++i ) + { + gp_Pnt p = surface->Value( uv[i].X(), uv[i].Y() ); + int nbDiv = int( p.Distance( pnt.back() ) / theMinSegLen ); + for ( int iD = 1; iD < nbDiv; ++iD ) + { + double R = iD / double( nbDiv ); + gp_XY uvR = uv[i-1] * (1 - R) + uv[i] * R; + pnt.push_back( surface->Value( uvR.X(), uvR.Y() )); + } + pnt.push_back( p ); + } + // cout << "from salome.geom import geomBuilder" << endl; // cout << "geompy = geomBuilder.New(salome.myStudy)" << endl; - Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1, uv.size()); - for ( size_t i = 0; i < uv.size(); ++i ) + Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1, pnt.size()); + for ( size_t i = 0; i < pnt.size(); ++i ) { - gp_Pnt p = surface->Value( uv[i].X(), uv[i].Y() ); + gp_Pnt& p = pnt[i]; points->SetValue( i+1, p ); - //cout << "geompy.MakeVertex( "<< p.X()<<", " << p.Y()<<", " << p.Z()<<" )" << endl; + // cout << "geompy.MakeVertex( "<< p.X()<<", " << p.Y()<<", " << p.Z() + // <<" theName = 'p_" << i << "')" << endl; } GeomAPI_Interpolate interpol( points, /*isClosed=*/false, gp::Resolution()); @@ -657,6 +677,7 @@ namespace const SMESH_MAT2d::MedialAxis& theMA, const SinuousFace& theSinuFace, SMESH_Algo* the1dAlgo, + const double theMinSegLen, vector& theMAParams ) { // check if all EDGEs of one size are meshed, then MA discretization is not needed @@ -673,7 +694,7 @@ namespace return true; // discretization is not needed - TopoDS_Edge branchEdge = makeEdgeFromMA( theHelper, theMA ); + TopoDS_Edge branchEdge = makeEdgeFromMA( theHelper, theMA, theMinSegLen ); if ( branchEdge.IsNull() ) return false; @@ -810,9 +831,11 @@ namespace */ //================================================================================ - bool findVertex( NodePoint& theNodePnt, - const vector& theSinuEdges, - 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; @@ -826,8 +849,10 @@ namespace V = SMESH_MesherHelper::IthVertex( 0, theSinuEdges[ theNodePnt._edgeInd ], /*CumOri=*/false); else if ( Abs( l - theNodePnt._u ) < tol ) V = SMESH_MesherHelper::IthVertex( 1, theSinuEdges[ theNodePnt._edgeInd ], /*CumOri=*/false); + 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 ) @@ -836,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(); } //================================================================================ @@ -857,23 +881,43 @@ namespace //================================================================================ bool projectVertices( SMESH_MesherHelper& theHelper, - //const double theMinSegLen, const SMESH_MAT2d::MedialAxis& theMA, const vector< SMESH_MAT2d::BranchPoint >& theDivPoints, - const vector& theSinuEdges, - const vector< Handle(Geom_Curve) >& theCurves, + const vector< std::size_t > & theEdgeIDs1, + const vector< std::size_t > & theEdgeIDs2, 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.getBranches()[0]; + const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0); + + // 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; + 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; + 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 ) { @@ -882,10 +926,14 @@ namespace if ( !branch.getBoundaryPoints( theDivPoints[i], bp[0], bp[1] )) return false; - NodePoint np[2] = { NodePoint( bp[0] ), - NodePoint( bp[1] )}; - bool isVertex[2] = { findVertex( np[0], theSinuEdges, meshDS ), - findVertex( np[1], theSinuEdges, meshDS )}; + NodePoint np[2] = { + NodePoint( bp[0] ), + NodePoint( bp[1] ) + }; + bool isVertex[2] = { + 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 = thePointsOnE.insert( make_pair( uMA, make_pair( np[0], np[1]))).first; @@ -893,8 +941,8 @@ namespace if ( !isVertex[0] && !isVertex[1] ) return false; // error if ( isVertex[0] && isVertex[1] ) continue; - const size_t iVertex = isVertex[0] ? 0 : 1; - const size_t iNode = 1 - iVertex; + const size_t iVert = isVertex[0] ? 0 : 1; + const size_t iNode = 1 - iVert; bool isOppComputed = theIsEdgeComputed[ np[ iNode ]._edgeInd ]; if ( !isOppComputed ) @@ -912,6 +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; for ( int iS = 0; iS < 2; ++iS ) // side with Vertex and side with Nodes { NodePoint np = get( u2NP->second, iS ); @@ -926,6 +978,8 @@ namespace isShortPrev[iS] = ( r < rShort ); isShortNext[iS] = (( 1 - r ) > ( 1 - rShort )); } + // if ( !hasPrev ) isShortPrev[0] = isShortPrev[1] = false; + // if ( !hasNext ) isShortNext[0] = isShortNext[1] = false; map< double, pair< NodePoint, NodePoint > >::iterator u2NPClose; @@ -933,9 +987,9 @@ namespace ( isShortNext[0] && isShortNext[1] )) { u2NPClose = isShortPrev[0] ? u2NPPrev : u2NPNext; - NodePoint& npProj = get( u2NP->second, iNode ); // NP of VERTEX projection - NodePoint npCloseN = get( u2NPClose->second, iNode); // NP close to npProj - NodePoint npCloseV = get( u2NPClose->second, iVertex); // NP close to VERTEX + NodePoint& npProj = get( u2NP->second, iNode ); // NP of VERTEX projection + NodePoint npCloseN = get( u2NPClose->second, iNode ); // NP close to npProj + NodePoint npCloseV = get( u2NPClose->second, iVert ); // NP close to VERTEX if ( !npCloseV._node ) { npProj = npCloseN; @@ -947,24 +1001,139 @@ namespace // can't remove the neighbor projection as it is also from VERTEX, -> option 1) } } - // else option 1) - wide enough -> "duplicate" existing node + // else: option 1) - wide enough -> "duplicate" existing node { 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; + NodePoint& npProj = get( u2NP->second, iNode ); // NP of VERTEX projection + NodePoint& npCloseN = get( u2NPClose->second, iNode ); // NP close to npProj + 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 @@ -985,7 +1154,7 @@ namespace vector& theMAParams, SinuousFace& theSinuFace) { - if ( theMA.getBranches().size() != 1 ) + if ( theMA.nbBranches() != 1 ) return false; // normalize theMAParams @@ -997,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 ); @@ -1029,7 +1199,7 @@ namespace } } - const SMESH_MAT2d::Branch& branch = theMA.getBranches()[0]; + const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0); SMESH_MAT2d::BoundaryPoint bp[2]; vector< std::size_t > edgeIDs1, edgeIDs2; @@ -1121,10 +1291,12 @@ namespace ++iEdgePair; } - if ( !projectVertices( theHelper, theMA, divPoints, 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 ) @@ -1302,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(); @@ -1414,7 +1587,7 @@ bool StdMeshers_QuadFromMedialAxis_1D2D::Compute(SMESH_Mesh& theMesh, _regular1D->SetSegmentLength( minSegLen ); vector maParams; - if ( ! divideMA( helper, ma, sinuFace, _regular1D, maParams )) + if ( ! divideMA( helper, ma, sinuFace, _regular1D, minSegLen, maParams )) return error(COMPERR_BAD_SHAPE); _progress = 0.4;