X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_QuadFromMedialAxis_1D2D.cxx;h=c7df8a3f7f74812aeaea025cc325bde3b3f200d3;hp=55e0d19ab892068571b2b03d540d079dd03b719a;hb=4b9725e701864a6b59195f92a1de6b4baa2b1da8;hpb=7a65c9fad427b1ccba6b9ccae612296e5092a324 diff --git a/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx index 55e0d19ab..c7df8a3f7 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 @@ -25,6 +25,7 @@ #include "StdMeshers_QuadFromMedialAxis_1D2D.hxx" +#include "SMESHDS_Mesh.hxx" #include "SMESH_Block.hxx" #include "SMESH_Gen.hxx" #include "SMESH_MAT2d.hxx" @@ -61,6 +62,8 @@ #include #include +using namespace std; + //================================================================================ /*! * \brief 1D algo @@ -119,7 +122,7 @@ public: if ( !StdMeshers_Regular_1D::computeInternalParameters( mesh, C3D, len, f, l, theParams, false)) { for ( size_t i = 1; i < 15; ++i ) - theParams.push_back( i/15 ); + theParams.push_back( i/15. ); // ???? } else { @@ -158,7 +161,7 @@ StdMeshers_QuadFromMedialAxis_1D2D::StdMeshers_QuadFromMedialAxis_1D2D(int _shapeType = (1 << TopAbs_FACE); _onlyUnaryInput = true; // FACE by FACE so far _requireDiscreteBoundary = false; // make 1D by myself - _supportSubmeshes = true; // make 1D by myself + _supportSubmeshes = true; // make 1D by myself _neededLowerHyps[ 1 ] = true; // suppress warning on hiding a global 1D algo _neededLowerHyps[ 2 ] = true; // suppress warning on hiding a global 2D algo _compatibleHypothesis.clear(); @@ -380,7 +383,7 @@ namespace tmpMesh.ShapeToMesh( theEdges[i] ); try { if ( !mesh->GetGen() ) continue; // tmp mesh - mesh->GetGen()->Compute( tmpMesh, theEdges[i], true, true ); // make nodes on VERTEXes + mesh->GetGen()->Compute( tmpMesh, theEdges[i], SMESH_Gen::SHAPE_ONLY_UPWARD ); // make nodes on VERTEXes if ( !algo->Compute( tmpMesh, theEdges[i] )) continue; } @@ -552,7 +555,7 @@ namespace return false; for ( size_t iS = 0; iS < theShortEdges[ nbBranchPoints > 0 ].size(); ++iS ) - shortMap.Add( theShortEdges[ nbBranchPoints ][ iS ]); + shortMap.Add( theShortEdges[ nbBranchPoints > 0 ][ iS ]); ++nbBranchPoints; } @@ -865,7 +868,7 @@ namespace TmpMesh tmpMesh; tmpMesh.ShapeToMesh( branchEdge ); try { - mesh->GetGen()->Compute( tmpMesh, branchEdge, true, true ); // make nodes on VERTEXes + mesh->GetGen()->Compute( tmpMesh, branchEdge, SMESH_Gen::SHAPE_ONLY_UPWARD ); // make nodes on VERTEXes if ( !algo->Compute( tmpMesh, branchEdge )) return false; } @@ -1069,7 +1072,7 @@ namespace if ( ! ( theDivPoints[0]._iEdge == 0 && theDivPoints[0]._edgeParam == 0. )) // recursive call { - SMESH_MAT2d::BranchPoint brp( &branch, 0, 0 ); + 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(); @@ -1081,9 +1084,10 @@ namespace } } - // project theDivPoints + // project theDivPoints and keep projections to merge TMAPar2NPoints::iterator u2NP; + vector< TMAPar2NPoints::iterator > projToMerge; for ( size_t i = 0; i < theDivPoints.size(); ++i ) { if ( !branch.getParameter( theDivPoints[i], uMA )) @@ -1126,9 +1130,18 @@ namespace if ( isVertex[0] && isVertex[1] ) continue; - bool isOppComputed = theIsEdgeComputed[ np[ iNode ]._edgeInd ]; - if ( !isOppComputed ) - continue; + // bool isOppComputed = theIsEdgeComputed[ np[ iNode ]._edgeInd ]; + // if ( isOppComputed ) + projToMerge.push_back( u2NP ); + } + + // merge projections + + for ( size_t i = 0; i < projToMerge.size(); ++i ) + { + u2NP = projToMerge[i]; + const size_t iVert = get( u2NP->second, 0 )._node ? 0 : 1; // side with a VERTEX + const size_t iNode = 1 - iVert; // opposite (meshed?) side // a VERTEX is projected on a meshed EDGE; there are two options: // 1) a projected point is joined with a closet node if a strip between this and neighbor @@ -1142,10 +1155,6 @@ namespace bool isShortPrev[2], isShortNext[2], isPrevCloser[2]; TMAPar2NPoints::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 ); @@ -1158,11 +1167,9 @@ namespace double distNext = p.Distance( pNext ); double r = distPrev / ( distPrev + distNext ); isShortPrev [iS] = ( r < rShort ); - isShortNext [iS] = (( 1 - r ) > ( 1 - rShort )); - isPrevCloser[iS] = (( r < 0.5 ) && ( u2NPPrev->first > 0 )); + isShortNext [iS] = (( 1 - r ) < rShort ); + isPrevCloser[iS] = (( r < 0.5 ) && ( theSinuFace.IsRing() || u2NPPrev->first > 0 )); } - // if ( !hasPrev ) isShortPrev[0] = isShortPrev[1] = false; - // if ( !hasNext ) isShortNext[0] = isShortNext[1] = false; TMAPar2NPoints::iterator u2NPClose; @@ -1171,17 +1178,18 @@ namespace { u2NPClose = isPrevCloser[0] ? u2NPPrev : u2NPNext; NodePoint& npProj = get( u2NP->second, iNode ); // NP of VERTEX projection + NodePoint& npVert = get( u2NP->second, iVert ); // NP of VERTEX NodePoint npCloseN = get( u2NPClose->second, iNode ); // NP close to npProj - NodePoint npCloseV = get( u2NPClose->second, iVert ); // NP close to VERTEX - if ( !npCloseV._node ) + NodePoint npCloseV = get( u2NPClose->second, iVert ); // NP close to npVert + if ( !npCloseV._node || npCloseV._node == npVert._node ) { npProj = npCloseN; - thePointsOnE.erase( isPrevCloser[0] ? u2NPPrev : u2NPNext ); + thePointsOnE.erase( u2NPClose ); continue; } else { - // can't remove the neighbor projection as it is also from VERTEX, -> option 1) + // can't remove the neighbor projection as it is also from VERTEX -> option 1) } } // else: option 1) - wide enough -> "duplicate" existing node @@ -1361,6 +1369,116 @@ namespace return; } // separateNodes() + + //================================================================================ + /*! + * \brief Find association of nodes existing on the sinuous sides of a ring + * + * TMAPar2NPoints filled here is used in setQuadSides() only if theSinuFace.IsRing() + * to find most distant nodes of the inner and the outer wires + */ + //================================================================================ + + void assocNodes( SMESH_MesherHelper& theHelper, + SinuousFace& theSinuFace, + const SMESH_MAT2d::MedialAxis& theMA, + TMAPar2NPoints & thePointsOnE ) + { + list< TopoDS_Edge > ee1( theSinuFace._sinuSide [0].begin(), theSinuFace._sinuSide [0].end() ); + list< TopoDS_Edge > ee2( theSinuFace._sinuSide [1].begin(), theSinuFace._sinuSide [1].end() ); + StdMeshers_FaceSide sideOut( theSinuFace.Face(), ee1, theHelper.GetMesh(), true, true ); + StdMeshers_FaceSide sideIn ( theSinuFace.Face(), ee2, theHelper.GetMesh(), true, true ); + const UVPtStructVec& uvsOut = sideOut.GetUVPtStruct(); + const UVPtStructVec& uvsIn = sideIn.GetUVPtStruct(); + // if ( uvs1.size() != uvs2.size() ) + // return; + + const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0); + SMESH_MAT2d::BoundaryPoint bp[2]; + SMESH_MAT2d::BranchPoint brp; + SMESHDS_Mesh* meshDS = theHelper.GetMeshDS(); + + map< double, const SMDS_MeshNode* > nodeParams; // params of existing nodes + map< double, const SMDS_MeshNode* >::iterator u2n; + + // find a node of sideOut most distant from sideIn + + vector< BRepAdaptor_Curve > curvesIn( theSinuFace._sinuSide[1].size() ); + for ( size_t iE = 0; iE < theSinuFace._sinuSide[1].size(); ++iE ) + curvesIn[ iE ].Initialize( theSinuFace._sinuSide[1][iE] ); + + double maxDist = 0; + SMESH_MAT2d::BoundaryPoint bpIn; // closest IN point + const SMDS_MeshNode* nOut = 0; + const size_t nbEOut = theSinuFace._sinuSide[0].size(); + for ( size_t iE = 0; iE < nbEOut; ++iE ) + { + const TopoDS_Edge& E = theSinuFace._sinuSide[0][iE]; + + if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, E, /*skipMedium=*/true, nodeParams )) + return; + for ( u2n = nodeParams.begin(); u2n != nodeParams.end(); ++u2n ) + { + // point on EDGE (u2n) --> MA point (brp) + if ( !theMA.getBoundary().getBranchPoint( iE, u2n->first, brp ) || + !branch.getBoundaryPoints( brp, bp[0], bp[1] )) + return; + gp_Pnt pOut = SMESH_TNodeXYZ( u2n->second ); + gp_Pnt pIn = curvesIn[ bp[1]._edgeIndex - nbEOut ].Value( bp[1]._param ); + double dist = pOut.SquareDistance( pIn ); + if ( dist > maxDist ) + { + maxDist = dist; + nOut = u2n->second; + bpIn = bp[1]; + } + } + } + const SMDS_MeshNode* nIn = 0; + if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, + theSinuFace._sinuEdges[ bpIn._edgeIndex ], + /*skipMedium=*/true, + nodeParams )) + return; + u2n = nodeParams.lower_bound( bpIn._param ); + if ( u2n == nodeParams.end() ) + nIn = nodeParams.rbegin()->second; + else + nIn = u2n->second; + + // find position of distant nodes in uvsOut and uvsIn + size_t iDistOut, iDistIn; + for ( iDistOut = 0; iDistOut < uvsOut.size(); ++iDistOut ) + { + if ( uvsOut[iDistOut].node == nOut ) + break; + } + for ( iDistIn = 0; iDistIn < uvsIn.size(); ++iDistIn ) + { + if ( uvsIn[iDistIn].node == nIn ) + break; + } + if ( iDistOut == uvsOut.size() || iDistIn == uvsIn.size() ) + return; + + // store opposite nodes in thePointsOnE (param and EDGE have no sense) + pair< NodePoint, NodePoint > oppNodes( NodePoint( nOut, 0, 0 ), NodePoint( nIn, 0, 0)); + thePointsOnE.insert( make_pair( uvsOut[ iDistOut ].normParam, oppNodes )); + int iOut = iDistOut, iIn = iDistIn; + int i, nbNodes = std::min( uvsOut.size(), uvsIn.size() ); + if ( nbNodes > 5 ) nbNodes = 5; + for ( i = 0, ++iOut, --iIn; i < nbNodes; ++iOut, --iIn, ++i ) + { + iOut = theHelper.WrapIndex( iOut, uvsOut.size() ); + iIn = theHelper.WrapIndex( iIn, uvsIn.size() ); + oppNodes.first._node = uvsOut[ iOut ].node; + oppNodes.second._node = uvsIn[ iIn ].node; + thePointsOnE.insert( make_pair( uvsOut[ iOut ].normParam, oppNodes )); + } + + return; + } // assocNodes() + //================================================================================ /*! * \brief Setup sides of SinuousFace::_quad @@ -1385,9 +1503,9 @@ namespace list< TopoDS_Edge > side[4]; side[0].insert( side[0].end(), theFace._shortSide[0].begin(), theFace._shortSide[0].end() ); - side[1].insert( side[1].end(), theFace._sinuSide[1].begin(), theFace._sinuSide[1].end() ); + side[1].insert( side[1].end(), theFace._sinuSide [1].begin(), theFace._sinuSide [1].end() ); side[2].insert( side[2].end(), theFace._shortSide[1].begin(), theFace._shortSide[1].end() ); - side[3].insert( side[3].end(), theFace._sinuSide[0].begin(), theFace._sinuSide[0].end() ); + side[3].insert( side[3].end(), theFace._sinuSide [0].begin(), theFace._sinuSide [0].end() ); for ( int i = 0; i < 4; ++i ) { @@ -1404,6 +1522,11 @@ namespace if ( thePointsOnEdges.size() < 4 ) return false; + int nbOut = theFace._quad->side[ 1 ].GetUVPtStruct().size(); + int nbIn = theFace._quad->side[ 3 ].GetUVPtStruct().size(); + if ( nbOut == 0 || nbIn == 0 ) + return false; + // find most distant opposite nodes double maxDist = 0, dist; TMAPar2NPoints::const_iterator u2NPdist, u2NP = thePointsOnEdges.begin(); @@ -1414,7 +1537,7 @@ namespace if ( dist > maxDist ) { u2NPdist = u2NP; - maxDist = dist; + maxDist = dist; } } // compute distribution of radial nodes @@ -1426,6 +1549,8 @@ namespace params ); // add a radial quad side + + theHelper.SetElementsOnShape( true ); u2NP = thePointsOnEdges.begin(); const SMDS_MeshNode* nOut = u2NP->second.first._node; const SMDS_MeshNode* nIn = u2NP->second.second._node; @@ -1441,7 +1566,7 @@ namespace uvsNew.push_back( uvPt ); for (list::iterator itU = params.begin(); itU != params.end(); ++itU ) { - gp_XY uv = ( 1 - *itU ) * uvOut + *itU * uvIn; + gp_XY uv = ( 1 - *itU ) * uvOut + *itU * uvIn; // applied in direction Out -> In gp_Pnt p = surface->Value( uv.X(), uv.Y() ); uvPt.node = theHelper.AddNode( p.X(), p.Y(), p.Z(), /*id=*/0, uv.X(), uv.Y() ); uvPt.u = uv.X(); @@ -1455,13 +1580,10 @@ 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; + if ( nbIn != nbOut ) + theFace._quad->side[ 2 ] = StdMeshers_FaceSide::New( uvsNew ); // 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 @@ -1483,6 +1605,7 @@ namespace } // 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 ) { @@ -1569,208 +1692,223 @@ namespace const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0); SMESH_MAT2d::BoundaryPoint bp[2]; - vector< std::size_t > edgeIDs1, edgeIDs2; // indices in theSinuEdges - vector< SMESH_MAT2d::BranchPoint > divPoints; - if ( !allComputed ) - branch.getOppositeGeomEdges( edgeIDs1, edgeIDs2, divPoints ); - - for ( size_t i = 0; i < edgeIDs1.size(); ++i ) - if ( isComputed[ edgeIDs1[i]] && - isComputed[ edgeIDs2[i]] ) - { - int nbNodes1 = meshDS->MeshElements(edgeIDs[ edgeIDs1[i]] )->NbNodes(); - int nbNodes2 = meshDS->MeshElements(edgeIDs[ edgeIDs2[i]] )->NbNodes(); - if ( nbNodes1 != nbNodes2 ) - return false; - if (( i-1 >= 0 ) && - ( edgeIDs1[i-1] == edgeIDs1[i] || - edgeIDs2[i-1] == edgeIDs2[i] )) - return false; - if (( i+1 < edgeIDs1.size() ) && - ( edgeIDs1[i+1] == edgeIDs1[i] || - edgeIDs2[i+1] == edgeIDs2[i] )) - return false; - } - - // map (param on MA) to (parameters of nodes on a pair of theSinuEdges) TMAPar2NPoints pointsOnE; - vector maParams; - set projectedEdges; // treated EDGEs which 'isComputed' + // check that computed EDGEs are opposite and equally meshed + if ( allComputed ) + { + // int nbNodes[2] = { 0, 0 }; + // for ( int iSide = 0; iSide < 2; ++iSide ) // loop on two sinuous sides + // nbNodes[ iSide ] += meshDS->MeshElements( theSinuFace._sinuSide[ iSide ])->NbNodes() - 1; - // compute params of nodes on EDGEs by projecting division points from MA + // if ( nbNodes[0] != nbNodes[1] ) + // return false; - for ( size_t iEdgePair = 0; iEdgePair < edgeIDs1.size(); ++iEdgePair ) - // loop on pairs of opposite EDGEs + if ( theSinuFace.IsRing() ) + assocNodes( theHelper, theSinuFace, theMA, pointsOnE ); + } + else { - if ( projectedEdges.count( edgeIDs1[ iEdgePair ]) || - projectedEdges.count( edgeIDs2[ iEdgePair ]) ) - continue; + vector< std::size_t > edgeIDs1, edgeIDs2; // indices in theSinuEdges + vector< SMESH_MAT2d::BranchPoint > divPoints; + branch.getOppositeGeomEdges( edgeIDs1, edgeIDs2, divPoints ); + + for ( size_t i = 0; i < edgeIDs1.size(); ++i ) + if ( isComputed[ edgeIDs1[i]] && + isComputed[ edgeIDs2[i]] ) + { + int nbNodes1 = meshDS->MeshElements(edgeIDs[ edgeIDs1[i]] )->NbNodes(); + int nbNodes2 = meshDS->MeshElements(edgeIDs[ edgeIDs2[i]] )->NbNodes(); + if ( nbNodes1 != nbNodes2 ) + return false; + if (( int(i)-1 >= 0 ) && + ( edgeIDs1[i-1] == edgeIDs1[i] || + edgeIDs2[i-1] == edgeIDs2[i] )) + return false; + if (( i+1 < edgeIDs1.size() ) && + ( edgeIDs1[i+1] == edgeIDs1[i] || + edgeIDs2[i+1] == edgeIDs2[i] )) + return false; + } + + // map (param on MA) to (parameters of nodes on a pair of theSinuEdges) + vector maParams; + set projectedEdges; // treated EDGEs which 'isComputed' + + // compute params of nodes on EDGEs by projecting division points from MA - // -------------------------------------------------------------------------------- - if ( isComputed[ edgeIDs1[ iEdgePair ]] != // one EDGE is meshed - isComputed[ edgeIDs2[ iEdgePair ]]) + for ( size_t iEdgePair = 0; iEdgePair < edgeIDs1.size(); ++iEdgePair ) + // loop on pairs of opposite EDGEs { - // "projection" from one side to the other + if ( projectedEdges.count( edgeIDs1[ iEdgePair ]) || + projectedEdges.count( edgeIDs2[ iEdgePair ]) ) + continue; - size_t iEdgeComputed = edgeIDs1[iEdgePair], iSideComputed = 0; - if ( !isComputed[ iEdgeComputed ]) - ++iSideComputed, iEdgeComputed = edgeIDs2[iEdgePair]; + // -------------------------------------------------------------------------------- + if ( isComputed[ edgeIDs1[ iEdgePair ]] != // one EDGE is meshed + isComputed[ edgeIDs2[ iEdgePair ]]) + { + // "projection" from one side to the other - map< double, const SMDS_MeshNode* > nodeParams; // params of existing nodes - if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, theSinuEdges[ iEdgeComputed ], /*skipMedium=*/true, nodeParams )) - return false; + size_t iEdgeComputed = edgeIDs1[iEdgePair], iSideComputed = 0; + if ( !isComputed[ iEdgeComputed ]) + ++iSideComputed, iEdgeComputed = edgeIDs2[iEdgePair]; - projectedEdges.insert( iEdgeComputed ); + map< double, const SMDS_MeshNode* > nodeParams; // params of existing nodes + if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, theSinuEdges[ iEdgeComputed ], /*skipMedium=*/true, nodeParams )) + return false; - SMESH_MAT2d::BoundaryPoint& bndPnt = bp[ 1-iSideComputed ]; - SMESH_MAT2d::BranchPoint brp; - NodePoint npN, npB; // NodePoint's initialized by node and BoundaryPoint - NodePoint& np0 = iSideComputed ? npB : npN; - NodePoint& np1 = iSideComputed ? npN : npB; + projectedEdges.insert( iEdgeComputed ); - double maParam1st, maParamLast, maParam; - if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.begin()->first, brp )) - return false; - branch.getParameter( brp, maParam1st ); - if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.rbegin()->first, brp )) - return false; - branch.getParameter( brp, maParamLast ); + SMESH_MAT2d::BoundaryPoint& bndPnt = bp[ 1-iSideComputed ]; + SMESH_MAT2d::BranchPoint brp; + NodePoint npN, npB; // NodePoint's initialized by node and BoundaryPoint + NodePoint& np0 = iSideComputed ? npB : npN; + NodePoint& np1 = iSideComputed ? npN : npB; - map< double, const SMDS_MeshNode* >::iterator u2n = nodeParams.begin(), u2nEnd = nodeParams.end(); - TMAPar2NPoints::iterator end = pointsOnE.end(), pos = end; - TMAPar2NPoints::iterator & hint = (maParamLast > maParam1st) ? end : pos; - for ( ++u2n, --u2nEnd; u2n != u2nEnd; ++u2n ) - { - // point on EDGE (u2n) --> MA point (brp) - if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, u2n->first, brp )) + double maParam1st, maParamLast, maParam; + if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.begin()->first, brp )) return false; - // MA point --> points on 2 EDGEs (bp) - if ( !branch.getBoundaryPoints( brp, bp[0], bp[1] ) || - !branch.getParameter( brp, maParam )) + branch.getParameter( brp, maParam1st ); + if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.rbegin()->first, brp )) return false; + branch.getParameter( brp, maParamLast ); - npN = NodePoint( u2n->second, u2n->first, iEdgeComputed ); - npB = NodePoint( bndPnt ); - pos = pointsOnE.insert( hint, make_pair( maParam, make_pair( np0, np1 ))); + map< double, const SMDS_MeshNode* >::iterator u2n = nodeParams.begin(), u2nEnd = nodeParams.end(); + TMAPar2NPoints::iterator end = pointsOnE.end(), pos = end; + TMAPar2NPoints::iterator & hint = (maParamLast > maParam1st) ? end : pos; + for ( ++u2n, --u2nEnd; u2n != u2nEnd; ++u2n ) + { + // point on EDGE (u2n) --> MA point (brp) + if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, u2n->first, brp )) + return false; + // MA point --> points on 2 EDGEs (bp) + if ( !branch.getBoundaryPoints( brp, bp[0], bp[1] ) || + !branch.getParameter( brp, maParam )) + return false; + + npN = NodePoint( u2n->second, u2n->first, iEdgeComputed ); + npB = NodePoint( bndPnt ); + pos = pointsOnE.insert( hint, make_pair( maParam, make_pair( np0, np1 ))); + } } - } - // -------------------------------------------------------------------------------- - else if ( !isComputed[ edgeIDs1[ iEdgePair ]] && // none of EDGEs is meshed - !isComputed[ edgeIDs2[ iEdgePair ]]) - { - // "projection" from MA - maParams.clear(); - if ( !getParamsForEdgePair( iEdgePair, divPoints, theMAParams, maParams )) - return false; - - for ( size_t i = 1; i < maParams.size()-1; ++i ) + // -------------------------------------------------------------------------------- + else if ( !isComputed[ edgeIDs1[ iEdgePair ]] && // none of EDGEs is meshed + !isComputed[ edgeIDs2[ iEdgePair ]]) { - if ( !branch.getBoundaryPoints( maParams[i], bp[0], bp[1] )) + // "projection" from MA + maParams.clear(); + if ( !getParamsForEdgePair( iEdgePair, divPoints, theMAParams, maParams )) return false; - pointsOnE.insert( pointsOnE.end(), make_pair( maParams[i], make_pair( NodePoint(bp[0]), - NodePoint(bp[1])))); - } - } - // -------------------------------------------------------------------------------- - else if ( isComputed[ edgeIDs1[ iEdgePair ]] && // equally meshed EDGES - isComputed[ edgeIDs2[ iEdgePair ]]) - { - // add existing nodes - - size_t iE0 = edgeIDs1[ iEdgePair ]; - size_t iE1 = edgeIDs2[ iEdgePair ]; - map< double, const SMDS_MeshNode* > nodeParams[2]; // params of existing nodes - if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, theSinuEdges[ iE0 ], - /*skipMedium=*/false, nodeParams[0] ) || - !SMESH_Algo::GetSortedNodesOnEdge( meshDS, theSinuEdges[ iE1 ], - /*skipMedium=*/false, nodeParams[1] ) || - nodeParams[0].size() != nodeParams[1].size() ) - return false; - - if ( nodeParams[0].size() <= 2 ) - continue; // nodes on VERTEXes only + for ( size_t i = 1; i < maParams.size()-1; ++i ) + { + if ( !branch.getBoundaryPoints( maParams[i], bp[0], bp[1] )) + return false; - bool reverse = ( theSinuEdges[0].Orientation() == theSinuEdges[1].Orientation() ); - double maParam; - SMESH_MAT2d::BranchPoint brp; - std::pair< NodePoint, NodePoint > npPair; - - map< double, const SMDS_MeshNode* >::iterator - u2n0F = ++nodeParams[0].begin(), - u2n1F = ++nodeParams[1].begin(); - map< double, const SMDS_MeshNode* >::reverse_iterator - u2n1R = ++nodeParams[1].rbegin(); - for ( ; u2n0F != nodeParams[0].end(); ++u2n0F ) + pointsOnE.insert( pointsOnE.end(), make_pair( maParams[i], make_pair( NodePoint(bp[0]), + NodePoint(bp[1])))); + } + } + // -------------------------------------------------------------------------------- + else if ( isComputed[ edgeIDs1[ iEdgePair ]] && // equally meshed EDGES + isComputed[ edgeIDs2[ iEdgePair ]]) { - if ( !theMA.getBoundary().getBranchPoint( iE0, u2n0F->first, brp ) || - !branch.getParameter( brp, maParam )) + // add existing nodes + + size_t iE0 = edgeIDs1[ iEdgePair ]; + size_t iE1 = edgeIDs2[ iEdgePair ]; + map< double, const SMDS_MeshNode* > nodeParams[2]; // params of existing nodes + if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, theSinuEdges[ iE0 ], + /*skipMedium=*/false, nodeParams[0] ) || + !SMESH_Algo::GetSortedNodesOnEdge( meshDS, theSinuEdges[ iE1 ], + /*skipMedium=*/false, nodeParams[1] ) || + nodeParams[0].size() != nodeParams[1].size() ) return false; - npPair.first = NodePoint( u2n0F->second, u2n0F->first, iE0 ); - if ( reverse ) - { - npPair.second = NodePoint( u2n1R->second, u2n1R->first, iE1 ); - ++u2n1R; - } - else + if ( nodeParams[0].size() <= 2 ) + continue; // nodes on VERTEXes only + + bool reverse = ( theSinuEdges[0].Orientation() == theSinuEdges[1].Orientation() ); + double maParam; + SMESH_MAT2d::BranchPoint brp; + std::pair< NodePoint, NodePoint > npPair; + + map< double, const SMDS_MeshNode* >::iterator + u2n0F = ++nodeParams[0].begin(), + u2n1F = ++nodeParams[1].begin(); + map< double, const SMDS_MeshNode* >::reverse_iterator + u2n1R = ++nodeParams[1].rbegin(); + for ( ; u2n0F != nodeParams[0].end(); ++u2n0F ) { - npPair.second = NodePoint( u2n1F->second, u2n1F->first, iE1 ); - ++u2n1F; + if ( !theMA.getBoundary().getBranchPoint( iE0, u2n0F->first, brp ) || + !branch.getParameter( brp, maParam )) + return false; + + npPair.first = NodePoint( u2n0F->second, u2n0F->first, iE0 ); + if ( reverse ) + { + npPair.second = NodePoint( u2n1R->second, u2n1R->first, iE1 ); + ++u2n1R; + } + else + { + npPair.second = NodePoint( u2n1F->second, u2n1F->first, iE1 ); + ++u2n1F; + } + pointsOnE.insert( make_pair( maParam, npPair )); } - pointsOnE.insert( make_pair( maParam, npPair )); } - } - } // loop on pairs of opposite EDGEs + } // loop on pairs of opposite EDGEs - if ( !projectVertices( theHelper, theMA, divPoints, edgeIDs1, edgeIDs2, - isComputed, pointsOnE, theSinuFace )) - return false; + if ( !projectVertices( theHelper, theMA, divPoints, edgeIDs1, edgeIDs2, + isComputed, pointsOnE, theSinuFace )) + return false; - separateNodes( theHelper, theMA, pointsOnE, theSinuFace, isComputed ); + separateNodes( theHelper, theMA, pointsOnE, theSinuFace, isComputed ); - // create nodes - TMAPar2NPoints::iterator u2np = pointsOnE.begin(); - for ( ; u2np != pointsOnE.end(); ++u2np ) - { - NodePoint* np[2] = { & u2np->second.first, & u2np->second.second }; - for ( int iSide = 0; iSide < 2; ++iSide ) + // create nodes + TMAPar2NPoints::iterator u2np = pointsOnE.begin(); + for ( ; u2np != pointsOnE.end(); ++u2np ) { - if ( np[ iSide ]->_node ) continue; - size_t iEdge = np[ iSide ]->_edgeInd; - double u = np[ iSide ]->_u; - gp_Pnt p = curves[ iEdge ]->Value( u ); - np[ iSide ]->_node = meshDS->AddNode( p.X(), p.Y(), p.Z() ); - meshDS->SetNodeOnEdge( np[ iSide ]->_node, edgeIDs[ iEdge ], u ); + NodePoint* np[2] = { & u2np->second.first, & u2np->second.second }; + for ( int iSide = 0; iSide < 2; ++iSide ) + { + if ( np[ iSide ]->_node ) continue; + size_t iEdge = np[ iSide ]->_edgeInd; + double u = np[ iSide ]->_u; + gp_Pnt p = curves[ iEdge ]->Value( u ); + np[ iSide ]->_node = meshDS->AddNode( p.X(), p.Y(), p.Z() ); + meshDS->SetNodeOnEdge( np[ iSide ]->_node, edgeIDs[ iEdge ], u ); + } } - } - // create mesh segments on EDGEs - theHelper.SetElementsOnShape( false ); - TopoDS_Face face = TopoDS::Face( theHelper.GetSubShape() ); - for ( size_t i = 0; i < theSinuEdges.size(); ++i ) - { - SMESH_subMesh* sm = mesh->GetSubMesh( theSinuEdges[i] ); - if ( sm->GetSubMeshDS() && sm->GetSubMeshDS()->NbElements() > 0 ) - continue; - - StdMeshers_FaceSide side( face, theSinuEdges[i], mesh, - /*isFwd=*/true, /*skipMediumNodes=*/true ); - vector nodes = side.GetOrderedNodes(); - for ( size_t in = 1; in < nodes.size(); ++in ) + // create mesh segments on EDGEs + theHelper.SetElementsOnShape( false ); + TopoDS_Face face = TopoDS::Face( theHelper.GetSubShape() ); + for ( size_t i = 0; i < theSinuEdges.size(); ++i ) { - const SMDS_MeshElement* seg = theHelper.AddEdge( nodes[in-1], nodes[in], 0, false ); - meshDS->SetMeshElementOnShape( seg, edgeIDs[ i ] ); + SMESH_subMesh* sm = mesh->GetSubMesh( theSinuEdges[i] ); + if ( sm->GetSubMeshDS() && sm->GetSubMeshDS()->NbElements() > 0 ) + continue; + + StdMeshers_FaceSide side( face, theSinuEdges[i], mesh, + /*isFwd=*/true, /*skipMediumNodes=*/true ); + vector nodes = side.GetOrderedNodes(); + for ( size_t in = 1; in < nodes.size(); ++in ) + { + const SMDS_MeshElement* seg = theHelper.AddEdge( nodes[in-1], nodes[in], 0, false ); + meshDS->SetMeshElementOnShape( seg, edgeIDs[ i ] ); + } } - } - // update sub-meshes on VERTEXes - for ( size_t i = 0; i < theSinuEdges.size(); ++i ) - { - mesh->GetSubMesh( theHelper.IthVertex( 0, theSinuEdges[i] )) - ->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); - mesh->GetSubMesh( theHelper.IthVertex( 1, theSinuEdges[i] )) - ->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); + // update sub-meshes on VERTEXes + for ( size_t i = 0; i < theSinuEdges.size(); ++i ) + { + mesh->GetSubMesh( theHelper.IthVertex( 0, theSinuEdges[i] )) + ->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); + mesh->GetSubMesh( theHelper.IthVertex( 1, theSinuEdges[i] )) + ->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); + } } // Setup sides of a quadrangle @@ -1797,7 +1935,8 @@ namespace { if ( !theHasRadialHyp ) // use global hyps - theHelper.GetGen()->Compute( *theHelper.GetMesh(), theShortEdges[i], true, true ); + theHelper.GetGen()->Compute( *theHelper.GetMesh(), theShortEdges[i], + SMESH_Gen::SHAPE_ONLY_UPWARD ); SMESH_subMesh* sm = theHelper.GetMesh()->GetSubMesh(theShortEdges[i] ); if ( sm->IsEmpty() ) @@ -1967,21 +2106,25 @@ bool StdMeshers_QuadFromMedialAxis_1D2D::computeQuads( SMESH_MesherHelper& theHe int nbNodesShort0 = theQuad->side[0].NbPoints(); int nbNodesShort1 = theQuad->side[2].NbPoints(); + int nbNodesSinu0 = theQuad->side[1].NbPoints(); + int nbNodesSinu1 = theQuad->side[3].NbPoints(); // compute UV of internal points myQuadList.push_back( theQuad ); - if ( !StdMeshers_Quadrangle_2D::setNormalizedGrid( theQuad )) - return false; + // if ( !StdMeshers_Quadrangle_2D::setNormalizedGrid( theQuad )) + // return false; // elliptic smooth of internal points to get boundary cell normal to the boundary bool isRing = theQuad->side[0].grid->Edge(0).IsNull(); - if ( !isRing ) + if ( !isRing ) { + if ( !StdMeshers_Quadrangle_2D::setNormalizedGrid( theQuad )) + return false; ellipticSmooth( theQuad, 1 ); - + } // create quadrangles bool ok; theHelper.SetElementsOnShape( true ); - if ( nbNodesShort0 == nbNodesShort1 ) + if ( nbNodesShort0 == nbNodesShort1 && nbNodesSinu0 == nbNodesSinu1 ) ok = StdMeshers_Quadrangle_2D::computeQuadDominant( *theHelper.GetMesh(), theQuad->face, theQuad ); else