X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;ds=sidebyside;f=src%2FStdMeshers%2FStdMeshers_QuadFromMedialAxis_1D2D.cxx;h=46d661191d4799816baa7972924315dbde517156;hb=7eda9ca931ed2a11cb5e4637e4ffe19f5c061115;hp=bc36be444166ea43080541212dbd8aec448bbc7e;hpb=27c5228fcfe39a6405378da49ac700ee2adca6cc;p=modules%2Fsmesh.git diff --git a/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx index bc36be444..46d661191 100644 --- a/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx +++ b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE // // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS @@ -61,6 +61,8 @@ #include #include +using namespace std; + //================================================================================ /*! * \brief 1D algo @@ -551,8 +553,8 @@ namespace allEdges, theShortEdges[ nbBranchPoints > 0 ] )) return false; - for ( size_t iS = 0; iS < theShortEdges[ nbBranchPoints ].size(); ++iS ) - shortMap.Add( theShortEdges[ nbBranchPoints ][ iS ]); + for ( size_t iS = 0; iS < theShortEdges[ nbBranchPoints > 0 ].size(); ++iS ) + shortMap.Add( theShortEdges[ nbBranchPoints > 0 ][ iS ]); ++nbBranchPoints; } @@ -949,7 +951,7 @@ namespace { const SMDS_MeshNode* _node; double _u; - int _edgeInd; // index in theSinuEdges vector + size_t _edgeInd; // index in theSinuEdges vector NodePoint(): _node(0), _u(0), _edgeInd(-1) {} NodePoint(const SMDS_MeshNode* n, double u, size_t iEdge ): _node(n), _u(u), _edgeInd(iEdge) {} @@ -1022,23 +1024,26 @@ namespace */ //================================================================================ - bool projectVertices( SMESH_MesherHelper& theHelper, - 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< bool >& theIsEdgeComputed, - TMAPar2NPoints & thePointsOnE, - SinuousFace& theSinuFace) + bool projectVertices( SMESH_MesherHelper& theHelper, + const SMESH_MAT2d::MedialAxis& theMA, + vector< SMESH_MAT2d::BranchPoint >& theDivPoints, + const vector< std::size_t > & theEdgeIDs1, + const vector< std::size_t > & theEdgeIDs2, + const vector< bool >& theIsEdgeComputed, + TMAPar2NPoints & thePointsOnE, + SinuousFace& theSinuFace) { + if ( theDivPoints.empty() ) + return true; + SMESHDS_Mesh* meshDS = theHelper.GetMeshDS(); - const vector& theSinuEdges = theSinuFace._sinuEdges; + const vector< TopoDS_Edge >& theSinuEdges = theSinuFace._sinuEdges; const vector< Handle(Geom_Curve) >& theCurves = theSinuFace._sinuCurves; double uMA; - SMESH_MAT2d::BoundaryPoint bp[2]; + SMESH_MAT2d::BoundaryPoint bp[2]; // 2 sinuous sides 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] )) return false; @@ -1048,22 +1053,39 @@ namespace findVertexAndNode( np0, theSinuEdges, meshDS ); findVertexAndNode( np1, theSinuEdges, meshDS ); thePointsOnE.insert( make_pair( -0.1, make_pair( np0, np1 ))); - + } if ( !theSinuFace.IsRing() ) { 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]; + NodePoint 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 + else + { + // project a VERTEX of outer sinuous side corresponding to branch(0.) + // which is not included into theDivPoints + if ( ! ( theDivPoints[0]._iEdge == 0 && + theDivPoints[0]._edgeParam == 0. )) // recursive call + { + SMESH_MAT2d::BranchPoint brp( &branch, 0, 0 ); + vector< SMESH_MAT2d::BranchPoint > divPoint( 1, brp ); + vector< std::size_t > edgeIDs1(2), edgeIDs2(2); + edgeIDs1[0] = theEdgeIDs1.back(); + edgeIDs1[1] = theEdgeIDs1[0]; + edgeIDs2[0] = theEdgeIDs2.back(); + edgeIDs2[1] = theEdgeIDs2[0]; + projectVertices( theHelper, theMA, divPoint, edgeIDs1, edgeIDs2, + theIsEdgeComputed, thePointsOnE, theSinuFace ); + } + } - if ( theDivPoints.empty() ) - return true; + // project theDivPoints + TMAPar2NPoints::iterator u2NP; for ( size_t i = 0; i < theDivPoints.size(); ++i ) { if ( !branch.getParameter( theDivPoints[i], uMA )) @@ -1079,15 +1101,32 @@ namespace findVertexAndNode( np[0], theSinuEdges, meshDS, theEdgeIDs1[i], theEdgeIDs1[i+1] ), findVertexAndNode( np[1], theSinuEdges, meshDS, theEdgeIDs2[i], theEdgeIDs2[i+1] ) }; + const size_t iVert = isVertex[0] ? 0 : 1; // side with a VERTEX + const size_t iNode = 1 - iVert; // opposite (meshed?) side - TMAPar2NPoints::iterator u2NP = - thePointsOnE.insert( make_pair( uMA, make_pair( np[0], np[1])));//.first; + if ( isVertex[0] != isVertex[1] ) // try to find an opposite VERTEX + { + theMA.getBoundary().moveToClosestEdgeEnd( bp[iNode] ); // EDGE -> VERTEX + SMESH_MAT2d::BranchPoint brp; + theMA.getBoundary().getBranchPoint( bp[iNode], brp ); // WIRE -> MA + SMESH_MAT2d::BoundaryPoint bp2[2]; + branch.getBoundaryPoints( brp, bp2[0], bp2[1] ); // MA -> WIRE + NodePoint np2[2] = { NodePoint( bp2[0]), NodePoint( bp2[1]) }; + findVertexAndNode( np2[0], theSinuEdges, meshDS ); + findVertexAndNode( np2[1], theSinuEdges, meshDS ); + if ( np2[ iVert ]._node == np[ iVert ]._node && + np2[ iNode ]._node) + { + np[ iNode ] = np2[ iNode ]; + isVertex[ iNode ] = true; + } + } + + u2NP = thePointsOnE.insert( make_pair( uMA, make_pair( np[0], np[1]))); if ( !isVertex[0] && !isVertex[1] ) return false; // error if ( isVertex[0] && isVertex[1] ) continue; - const size_t iVert = isVertex[0] ? 0 : 1; - const size_t iNode = 1 - iVert; bool isOppComputed = theIsEdgeComputed[ np[ iNode ]._edgeInd ]; if ( !isOppComputed ) @@ -1101,8 +1140,8 @@ namespace // projection is set to the BoundaryPoint of this projection // evaluate distance to neighbor projections - const double rShort = 0.2; - bool isShortPrev[2], isShortNext[2]; + const double rShort = 0.33; + bool isShortPrev[2], isShortNext[2], isPrevCloser[2]; TMAPar2NPoints::iterator u2NPPrev = u2NP, u2NPNext = u2NP; --u2NPPrev; ++u2NPNext; // bool hasPrev = ( u2NP != thePointsOnE.begin() ); @@ -1120,8 +1159,9 @@ namespace double distPrev = p.Distance( pPrev ); double distNext = p.Distance( pNext ); double r = distPrev / ( distPrev + distNext ); - isShortPrev[iS] = ( r < rShort ); - isShortNext[iS] = (( 1 - r ) > ( 1 - rShort )); + isShortPrev [iS] = ( r < rShort ); + isShortNext [iS] = (( 1 - r ) > ( 1 - rShort )); + isPrevCloser[iS] = (( r < 0.5 ) && ( u2NPPrev->first > 0 )); } // if ( !hasPrev ) isShortPrev[0] = isShortPrev[1] = false; // if ( !hasNext ) isShortNext[0] = isShortNext[1] = false; @@ -1131,14 +1171,14 @@ namespace if (( isShortPrev[0] && isShortPrev[1] ) || // option 2) -> remove a too close projection ( isShortNext[0] && isShortNext[1] )) { - u2NPClose = isShortPrev[0] ? u2NPPrev : u2NPNext; + u2NPClose = isPrevCloser[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, iVert ); // NP close to VERTEX if ( !npCloseV._node ) { npProj = npCloseN; - thePointsOnE.erase( isShortPrev[0] ? u2NPPrev : u2NPNext ); + thePointsOnE.erase( isPrevCloser[0] ? u2NPPrev : u2NPNext ); continue; } else @@ -1148,7 +1188,7 @@ namespace } // else: option 1) - wide enough -> "duplicate" existing node { - u2NPClose = isShortPrev[ iNode ] ? u2NPPrev : u2NPNext; + u2NPClose = isPrevCloser[ iNode ] ? u2NPPrev : u2NPNext; NodePoint& npProj = get( u2NP->second, iNode ); // NP of VERTEX projection NodePoint& npCloseN = get( u2NPClose->second, iNode ); // NP close to npProj npProj = npCloseN; @@ -1163,6 +1203,12 @@ namespace //theNodes2Merge[ npCloseN._node ].push_back( npProj._node ); } } + + // remove auxiliary NodePoint's of ends of theSinuEdges + for ( u2NP = thePointsOnE.begin(); u2NP->first < 0; ) + thePointsOnE.erase( u2NP++ ); + thePointsOnE.erase( 1.1 ); + return true; } @@ -1197,7 +1243,8 @@ namespace void separateNodes( SMESH_MesherHelper& theHelper, const SMESH_MAT2d::MedialAxis& theMA, TMAPar2NPoints & thePointsOnE, - SinuousFace& theSinuFace ) + SinuousFace& theSinuFace, + const vector< bool >& theIsComputedEdge) { if ( thePointsOnE.size() < 2 ) return; @@ -1206,8 +1253,8 @@ namespace const vector& theSinuEdges = theSinuFace._sinuEdges; const vector< Handle(Geom_Curve) >& curves = theSinuFace._sinuCurves; - SMESH_MAT2d::BoundaryPoint bp[2]; - const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0); + //SMESH_MAT2d::BoundaryPoint bp[2]; + //const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0); typedef TMAPar2NPoints::iterator TIterator; @@ -1247,7 +1294,7 @@ namespace { // find an existing node on VERTEX among sameU2NP and get underlying EDGEs const SMDS_MeshNode* existingNode = 0; - set< int > edgeInds; + set< size_t > edgeInds; NodePoint* np; for ( size_t i = 0; i < sameU2NP.size(); ++i ) { @@ -1261,10 +1308,10 @@ namespace TIterator u2NPprev = sameU2NP.front(); TIterator u2NPnext = sameU2NP.back() ; - if ( u2NPprev->first > 0. ) --u2NPprev; - if ( u2NPnext->first < 1. ) ++u2NPprev; + if ( u2NPprev->first < 0. ) ++u2NPprev; + if ( u2NPnext->first > 1. ) --u2NPnext; - set< int >::iterator edgeID = edgeInds.begin(); + set< size_t >::iterator edgeID = edgeInds.begin(); for ( ; edgeID != edgeInds.end(); ++edgeID ) { // get U range on iEdge within which the equal points will be distributed @@ -1275,8 +1322,18 @@ namespace np = &get( u2NPnext->second, iSide ); u1 = getUOnEdgeByPoint( *edgeID, np, theSinuFace ); + if ( u0 == u1 ) + { + if ( u2NPprev != thePointsOnE.begin() ) --u2NPprev; + if ( u2NPnext != --thePointsOnE.end() ) ++u2NPnext; + np = &get( u2NPprev->second, iSide ); + u0 = getUOnEdgeByPoint( *edgeID, np, theSinuFace ); + np = &get( u2NPnext->second, iSide ); + u1 = getUOnEdgeByPoint( *edgeID, np, theSinuFace ); + } + // distribute points and create nodes - double du = ( u1 - u0 ) / ( sameU2NP.size() + 1 ); + double du = ( u1 - u0 ) / ( sameU2NP.size() + 1 /*!existingNode*/ ); double u = u0 + du; for ( size_t i = 0; i < sameU2NP.size(); ++i ) { @@ -1288,7 +1345,9 @@ namespace gp_Pnt p = np->Point( curves ); np->_node = meshDS->AddNode( p.X(), p.Y(), p.Z() ); meshDS->SetNodeOnEdge( np->_node, theSinuEdges[ *edgeID ], np->_u ); - //mergeNodes.push_back( np->_node ); + + if ( theIsComputedEdge[ *edgeID ]) + mergeNodes.push_back( np->_node ); } } } @@ -1352,8 +1411,8 @@ namespace TMAPar2NPoints::const_iterator u2NPdist, u2NP = thePointsOnEdges.begin(); for ( ; u2NP != thePointsOnEdges.end(); ++u2NP ) { - SMESH_TNodeXYZ xyz( u2NP->second.first._node ); - dist = xyz.SquareDistance( u2NP->second.second._node ); + SMESH_TNodeXYZ xyz( u2NP->second.first._node ); // node out + dist = xyz.SquareDistance( u2NP->second.second._node );// node in if ( dist > maxDist ) { u2NPdist = u2NP; @@ -1399,6 +1458,32 @@ namespace theFace._quad->side[ 0 ] = StdMeshers_FaceSide::New( uvsNew ); theFace._quad->side[ 2 ] = theFace._quad->side[ 0 ]; + if ( theFace._quad->side[ 1 ].GetUVPtStruct().empty() || + theFace._quad->side[ 3 ].GetUVPtStruct().empty() ) + return false; + + // assure that the outer sinuous side starts at nOut + if ( theFace._sinuSide[0].size() > 1 ) + { + const UVPtStructVec& uvsOut = theFace._quad->side[ 3 ].GetUVPtStruct(); // _sinuSide[0] + size_t i; // find UVPtStruct holding nOut + for ( i = 0; i < uvsOut.size(); ++i ) + if ( nOut == uvsOut[i].node ) + break; + if ( i == uvsOut.size() ) + return false; + + if ( i != 0 && i != uvsOut.size()-1 ) + { + // create a new OUT quad side + uvsNew.clear(); + uvsNew.reserve( uvsOut.size() ); + uvsNew.insert( uvsNew.end(), uvsOut.begin() + i, uvsOut.end() ); + uvsNew.insert( uvsNew.end(), uvsOut.begin() + 1, uvsOut.begin() + i + 1); + theFace._quad->side[ 3 ] = StdMeshers_FaceSide::New( uvsNew ); + } + } + // rotate the IN side if opposite nodes of IN and OUT sides don't match const SMDS_MeshNode * nIn0 = theFace._quad->side[ 1 ].First().node; if ( nIn0 != nIn ) @@ -1419,7 +1504,12 @@ namespace uvsNew.insert( uvsNew.end(), uvsIn.begin() + 1, uvsIn.begin() + i + 1); theFace._quad->side[ 1 ] = StdMeshers_FaceSide::New( uvsNew ); } - } // if ( theShortEdges[0].empty() ) + + if ( theFace._quad->side[ 1 ].GetUVPtStruct().empty() || + theFace._quad->side[ 3 ].GetUVPtStruct().empty() ) + return false; + + } // if ( theFace.IsRing() ) return true; @@ -1465,6 +1555,7 @@ namespace vector< int > edgeIDs ( theSinuEdges.size() ); // IDs in the main shape vector< bool > isComputed( theSinuEdges.size() ); curves.resize( theSinuEdges.size(), 0 ); + bool allComputed = true; for ( size_t i = 0; i < theSinuEdges.size(); ++i ) { curves[i] = BRep_Tool::Curve( theSinuEdges[i], f,l ); @@ -1473,6 +1564,8 @@ namespace SMESH_subMesh* sm = mesh->GetSubMesh( theSinuEdges[i] ); edgeIDs [i] = sm->GetId(); isComputed[i] = ( !sm->IsEmpty() ); + if ( !isComputed[i] ) + allComputed = false; } const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0); @@ -1480,7 +1573,9 @@ namespace vector< std::size_t > edgeIDs1, edgeIDs2; // indices in theSinuEdges vector< SMESH_MAT2d::BranchPoint > divPoints; - branch.getOppositeGeomEdges( edgeIDs1, edgeIDs2, divPoints ); + if ( !allComputed ) + branch.getOppositeGeomEdges( edgeIDs1, edgeIDs2, divPoints ); + for ( size_t i = 0; i < edgeIDs1.size(); ++i ) if ( isComputed[ edgeIDs1[i]] && isComputed[ edgeIDs2[i]] ) @@ -1499,15 +1594,20 @@ namespace return false; } - // map param on MA to parameters of nodes on a pair of theSinuEdges + // map (param on MA) to (parameters of nodes on a pair of theSinuEdges) TMAPar2NPoints pointsOnE; vector maParams; + set projectedEdges; // treated EDGEs which 'isComputed' // compute params of nodes on EDGEs by projecting division points from MA for ( size_t iEdgePair = 0; iEdgePair < edgeIDs1.size(); ++iEdgePair ) // loop on pairs of opposite EDGEs { + if ( projectedEdges.count( edgeIDs1[ iEdgePair ]) || + projectedEdges.count( edgeIDs2[ iEdgePair ]) ) + continue; + // -------------------------------------------------------------------------------- if ( isComputed[ edgeIDs1[ iEdgePair ]] != // one EDGE is meshed isComputed[ edgeIDs2[ iEdgePair ]]) @@ -1522,6 +1622,8 @@ namespace if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, theSinuEdges[ iEdgeComputed ], /*skipMedium=*/true, nodeParams )) return false; + projectedEdges.insert( iEdgeComputed ); + SMESH_MAT2d::BoundaryPoint& bndPnt = bp[ 1-iSideComputed ]; SMESH_MAT2d::BranchPoint brp; NodePoint npN, npB; // NodePoint's initialized by node and BoundaryPoint @@ -1530,10 +1632,10 @@ namespace double maParam1st, maParamLast, maParam; if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.begin()->first, brp )) - return false; + return false; branch.getParameter( brp, maParam1st ); if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.rbegin()->first, brp )) - return false; + return false; branch.getParameter( brp, maParamLast ); map< double, const SMDS_MeshNode* >::iterator u2n = nodeParams.begin(), u2nEnd = nodeParams.end(); @@ -1553,28 +1655,6 @@ namespace npB = NodePoint( bndPnt ); pos = pointsOnE.insert( hint, make_pair( maParam, make_pair( np0, np1 ))); } - - // move iEdgePair forward; - // find divPoints most close to max MA param - if ( edgeIDs1.size() > 1 ) - { - maParamLast = pointsOnE.rbegin()->first; - int iClosest; - double minDist = 1.; - for ( ; iEdgePair < edgeIDs1.size()-1; ++iEdgePair ) - { - branch.getParameter( divPoints[iEdgePair], maParam ); - double d = Abs( maParamLast - maParam ); - if ( d < minDist ) - minDist = d, iClosest = iEdgePair; - else - break; - } - if ( Abs( maParamLast - 1. ) < minDist ) - break; // the last pair treated - else - iEdgePair = iClosest; - } } // -------------------------------------------------------------------------------- else if ( !isComputed[ edgeIDs1[ iEdgePair ]] && // none of EDGEs is meshed @@ -1649,7 +1729,7 @@ namespace isComputed, pointsOnE, theSinuFace )) return false; - separateNodes( theHelper, theMA, pointsOnE, theSinuFace ); + separateNodes( theHelper, theMA, pointsOnE, theSinuFace, isComputed ); // create nodes TMAPar2NPoints::iterator u2np = pointsOnE.begin(); @@ -1766,12 +1846,12 @@ namespace const double dksi = 0.5, deta = 0.5; const double dksi2 = dksi*dksi, deta2 = deta*deta; double err = 0., g11, g22, g12; - int nbErr = 0; + //int nbErr = 0; FaceQuadStruct& q = *quad; UVPtStruct pNew; - double refArea = area( q.UVPt(0,0), q.UVPt(1,0), q.UVPt(1,1) ); + //double refArea = area( q.UVPt(0,0), q.UVPt(1,0), q.UVPt(1,1) ); for ( int iLoop = 0; iLoop < nbLoops; ++iLoop ) { @@ -1902,6 +1982,7 @@ bool StdMeshers_QuadFromMedialAxis_1D2D::computeQuads( SMESH_MesherHelper& theHe // create quadrangles bool ok; + theHelper.SetElementsOnShape( true ); if ( nbNodesShort0 == nbNodesShort1 ) ok = StdMeshers_Quadrangle_2D::computeQuadDominant( *theHelper.GetMesh(), theQuad->face, theQuad );