X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_MEFISTO_2D.cxx;h=ca7405ba22c0673f41e258e715648256ff5f4a89;hp=b7da7cbe615ec8ffe2c8069932703e36e0cc6a67;hb=104ff7b2818ce4d0f8a38d840abd3e5c70190668;hpb=e4737e85f0da6d3f90fd08f6be1c2825195fe16f diff --git a/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx b/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx index b7da7cbe6..ca7405ba2 100644 --- a/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx +++ b/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx @@ -17,7 +17,7 @@ // License along with this library; if not, write to the Free Software // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // -// See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // // // @@ -62,9 +62,6 @@ using namespace std; #include #include -#include -//#include - //============================================================================= /*! * @@ -72,19 +69,19 @@ using namespace std; //============================================================================= StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D(int hypId, int studyId, - SMESH_Gen * gen):SMESH_2D_Algo(hypId, studyId, gen) + SMESH_Gen * gen):SMESH_2D_Algo(hypId, studyId, gen) { - MESSAGE("StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D"); - _name = "MEFISTO_2D"; -// _shapeType = TopAbs_FACE; - _shapeType = (1 << TopAbs_FACE); - _compatibleHypothesis.push_back("MaxElementArea"); - _compatibleHypothesis.push_back("LengthFromEdges"); - - _edgeLength = 0; - _maxElementArea = 0; - _hypMaxElementArea = NULL; - _hypLengthFromEdges = NULL; + MESSAGE("StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D"); + _name = "MEFISTO_2D"; + _shapeType = (1 << TopAbs_FACE); + _compatibleHypothesis.push_back("MaxElementArea"); + _compatibleHypothesis.push_back("LengthFromEdges"); + + _edgeLength = 0; + _maxElementArea = 0; + _hypMaxElementArea = NULL; + _hypLengthFromEdges = NULL; + myTool = 0; } //============================================================================= @@ -95,7 +92,7 @@ StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D(int hypId, int studyId, StdMeshers_MEFISTO_2D::~StdMeshers_MEFISTO_2D() { - MESSAGE("StdMeshers_MEFISTO_2D::~StdMeshers_MEFISTO_2D"); + MESSAGE("StdMeshers_MEFISTO_2D::~StdMeshers_MEFISTO_2D"); } //============================================================================= @@ -129,7 +126,7 @@ bool StdMeshers_MEFISTO_2D::CheckHypothesis theHyp = (*itl); // use only the first hypothesis string hypName = theHyp->GetName(); - int hypId = theHyp->GetID(); + //int hypId = theHyp->GetID(); //SCRUTE(hypName); bool isOk = false; @@ -184,112 +181,124 @@ bool StdMeshers_MEFISTO_2D::CheckHypothesis bool StdMeshers_MEFISTO_2D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape) { - MESSAGE("StdMeshers_MEFISTO_2D::Compute"); + MESSAGE("StdMeshers_MEFISTO_2D::Compute"); - if (_hypLengthFromEdges) - _edgeLength = ComputeEdgeElementLength(aMesh, aShape); + if (_hypLengthFromEdges) + _edgeLength = ComputeEdgeElementLength(aMesh, aShape); + + bool isOk = false; + //const SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); + //SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape); - bool isOk = false; - const SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); - SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape); - - const TopoDS_Face & FF = TopoDS::Face(aShape); - bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD); - TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD)); - - Z nblf; //nombre de lignes fermees (enveloppe en tete) - Z *nudslf = NULL; //numero du dernier sommet de chaque ligne fermee - R2 *uvslf = NULL; - Z nbpti = 0; //nombre points internes futurs sommets de la triangulation - R2 *uvpti = NULL; - - Z nbst; - R2 *uvst = NULL; - Z nbt; - Z *nust = NULL; - Z ierr = 0; - - Z nutysu = 1; // 1: il existe un fonction areteideale_() - // Z nutysu=0; // 0: on utilise aretmx - R aretmx = _edgeLength; // longueur max aretes future triangulation - - nblf = NumberOfWires(F); - - nudslf = new Z[1 + nblf]; - nudslf[0] = 0; - int iw = 1; - int nbpnt = 0; - - myOuterWire = BRepTools::OuterWire(F); - nbpnt += NumberOfPoints(aMesh, myOuterWire); - if ( nbpnt < 3 ) // ex: a circle with 2 segments - return false; - nudslf[iw++] = nbpnt; - - for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next()) - { - const TopoDS_Wire & W = TopoDS::Wire(exp.Current()); - if (!myOuterWire.IsSame(W)) - { - nbpnt += NumberOfPoints(aMesh, W); - nudslf[iw++] = nbpnt; - } - } + const TopoDS_Face & FF = TopoDS::Face(aShape); + bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD); + TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD)); - // avoid passing same uv points for a vertex common to 2 wires - TopTools_IndexedDataMapOfShapeListOfShape VWMap; - if ( iw - 1 > 1 ) // nbofWires > 1 - TopExp::MapShapesAndAncestors( F , TopAbs_VERTEX, TopAbs_WIRE, VWMap ); + Z nblf; //nombre de lignes fermees (enveloppe en tete) + Z *nudslf = NULL; //numero du dernier sommet de chaque ligne fermee + R2 *uvslf = NULL; + Z nbpti = 0; //nombre points internes futurs sommets de la triangulation + R2 *uvpti = NULL; + + Z nbst; + R2 *uvst = NULL; + Z nbt; + Z *nust = NULL; + Z ierr = 0; + + Z nutysu = 1; // 1: il existe un fonction areteideale_() + // Z nutysu=0; // 0: on utilise aretmx + R aretmx = _edgeLength; // longueur max aretes future triangulation + + nblf = NumberOfWires(F); + + nudslf = new Z[1 + nblf]; + nudslf[0] = 0; + int iw = 1; + int nbpnt = 0; + + myTool = new SMESH_MesherHelper(aMesh); + _quadraticMesh = myTool->IsQuadraticSubMesh(aShape); + + if ( _quadraticMesh && _hypLengthFromEdges ) + aretmx *= 2.; + + myOuterWire = BRepTools::OuterWire(F); + nbpnt += NumberOfPoints(aMesh, myOuterWire); + if ( nbpnt < 3 ) { // ex: a circle with 2 segments + delete myTool; myTool = 0; + return false; + } + nudslf[iw++] = nbpnt; - uvslf = new R2[nudslf[nblf]]; - int m = 0; + for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next()) { + const TopoDS_Wire & W = TopoDS::Wire(exp.Current()); + if (!myOuterWire.IsSame(W)) { + nbpnt += NumberOfPoints(aMesh, W); + nudslf[iw++] = nbpnt; + } + } - double scalex, scaley; - ComputeScaleOnFace(aMesh, F, scalex, scaley); + // avoid passing same uv points for a vertex common to 2 wires + TopTools_IndexedDataMapOfShapeListOfShape VWMap; + if ( iw - 1 > 1 ) // nbofWires > 1 + TopExp::MapShapesAndAncestors( F , TopAbs_VERTEX, TopAbs_WIRE, VWMap ); - map mefistoToDS; // correspondence mefisto index--> points IDNodes - if ( !LoadPoints(aMesh, F, myOuterWire, uvslf, m, - mefistoToDS, scalex, scaley, VWMap)) - return false; + uvslf = new R2[nudslf[nblf]]; + int m = 0; - for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next()) - { - const TopoDS_Wire & W = TopoDS::Wire(exp.Current()); - if (!myOuterWire.IsSame(W)) - { - if (! LoadPoints(aMesh, F, W, uvslf, m, - mefistoToDS, scalex, scaley, VWMap )) - return false; - } - } + double scalex, scaley; + ComputeScaleOnFace(aMesh, F, scalex, scaley); - uvst = NULL; - nust = NULL; - aptrte(nutysu, aretmx, - nblf, nudslf, uvslf, nbpti, uvpti, nbst, uvst, nbt, nust, ierr); - - if (ierr == 0) - { - MESSAGE("... End Triangulation Generated Triangle Number " << nbt); - MESSAGE(" Node Number " << nbst); - StoreResult(aMesh, nbst, uvst, nbt, nust, F, - faceIsForward, mefistoToDS, scalex, scaley); - isOk = true; - } - else - { - MESSAGE("Error in Triangulation"); - isOk = false; - } - if (nudslf != NULL) - delete[]nudslf; - if (uvslf != NULL) - delete[]uvslf; - if (uvst != NULL) - delete[]uvst; - if (nust != NULL) - delete[]nust; - return isOk; + map mefistoToDS; // correspondence mefisto index--> points IDNodes + if ( !LoadPoints(aMesh, F, myOuterWire, uvslf, m, + mefistoToDS, scalex, scaley, VWMap) ) { + delete myTool; myTool = 0; + return false; + } + + for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next()) + { + const TopoDS_Wire & W = TopoDS::Wire(exp.Current()); + if (!myOuterWire.IsSame(W)) + { + if (! LoadPoints(aMesh, F, W, uvslf, m, + mefistoToDS, scalex, scaley, VWMap )) { + delete myTool; myTool = 0; + return false; + } + } + } + + uvst = NULL; + nust = NULL; + aptrte(nutysu, aretmx, + nblf, nudslf, uvslf, nbpti, uvpti, nbst, uvst, nbt, nust, ierr); + + if (ierr == 0) + { + MESSAGE("... End Triangulation Generated Triangle Number " << nbt); + MESSAGE(" Node Number " << nbst); + StoreResult(aMesh, nbst, uvst, nbt, nust, F, + faceIsForward, mefistoToDS, scalex, scaley); + isOk = true; + } + else + { + MESSAGE("Error in Triangulation"); + isOk = false; + } + if (nudslf != NULL) + delete[]nudslf; + if (uvslf != NULL) + delete[]uvslf; + if (uvst != NULL) + delete[]uvst; + if (nust != NULL) + delete[]nust; + delete myTool; myTool = 0; + + return isOk; } //======================================================================= @@ -324,7 +333,8 @@ static bool fixOverlappedLinkUV( R2& uv0, const R2& uv1, const R2& uv2 ) uv0.y -= delta; else uv0.y += delta; -// MESSAGE(" -> " << uv0.x << " " << uv0.y << " "); +// #ifdef _DEBUG_ +// MESSAGE(" -> " << uv0.x << " " << uv0.y << " "); // MESSAGE("v1( " << v1.X() << " " << v1.Y() << " ) " << // "v2( " << v2.X() << " " << v2.Y() << " ) "); // MESSAGE("SIN: " << sqrt(1 - dot * dot / (sqMod1 * sqMod2))); @@ -336,6 +346,7 @@ static bool fixOverlappedLinkUV( R2& uv0, const R2& uv1, const R2& uv2 ) // dot = v1*v2; // double sin = sqrt(1 - dot * dot / (sqMod1 * sqMod2)); // MESSAGE("NEW SIN: " << sin); +// #endif return true; } return false; @@ -352,7 +363,8 @@ static bool fixCommonVertexUV (gp_Pnt2d & theUV, const TopoDS_Wire& theOW, const TopoDS_Face& theF, const TopTools_IndexedDataMapOfShapeListOfShape & theVWMap, - SMESH_Mesh & theMesh) + SMESH_Mesh & theMesh, + bool CreateQuadratic) { if( theW.IsSame( theOW ) || !theVWMap.Contains( theV )) return false; @@ -415,10 +427,12 @@ static bool fixCommonVertexUV (gp_Pnt2d & theUV, umin = l; umax = f; } - else - { + else { while ( nIt->more() ) { const SMDS_MeshNode* node = nIt->next(); + // check if node is medium + if ( CreateQuadratic && SMESH_MesherHelper::IsMedium( node, SMDSAbs_Edge )) + continue; const SMDS_EdgePosition* epos = static_cast(node->GetPosition().get()); double u = epos->GetUParameter(); @@ -475,134 +489,102 @@ bool StdMeshers_MEFISTO_2D::LoadPoints(SMESH_Mesh & aMesh, double scaley, const TopTools_IndexedDataMapOfShapeListOfShape& VWMap) { -// MESSAGE("StdMeshers_MEFISTO_2D::LoadPoints"); - - SMDS_Mesh * meshDS = aMesh.GetMeshDS(); + // MESSAGE("StdMeshers_MEFISTO_2D::LoadPoints"); TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD)); - int mInit = m, mFirst, iEdge; + list< int > mOnVertex; + gp_XY scale( scalex, scaley ); TopoDS_Wire W = TopoDS::Wire(WW.Oriented(TopAbs_FORWARD)); BRepTools_WireExplorer wexp(W, F); - for (wexp.Init(W, F), iEdge = 0; wexp.More(); wexp.Next(), iEdge++) + for ( wexp.Init(W, F); wexp.More(); wexp.Next() ) { const TopoDS_Edge & E = wexp.Current(); + bool isForward = (E.Orientation() == TopAbs_FORWARD); // --- IDNodes of first and last Vertex - TopoDS_Vertex VFirst, VLast; + TopoDS_Vertex VFirst, VLast, V; TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l + V = isForward ? VFirst : VLast; - ASSERT(!VFirst.IsNull()); - SMDS_NodeIteratorPtr lid= - aMesh.GetSubMesh(VFirst)->GetSubMeshDS()->GetNodes(); + ASSERT(!V.IsNull()); + SMDS_NodeIteratorPtr lid= aMesh.GetSubMesh(V)->GetSubMeshDS()->GetNodes(); if ( !lid->more() ) { MESSAGE (" NO NODE BUILT ON VERTEX "); return false; } const SMDS_MeshNode* idFirst = lid->next(); - ASSERT(!VLast.IsNull()); - lid=aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes(); - if ( !lid->more() ) { - MESSAGE (" NO NODE BUILT ON VERTEX "); - return false; - } - const SMDS_MeshNode* idLast = lid->next(); - // --- edge internal IDNodes (relies on good order storage, not checked) - int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes(); - - double f, l; - Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l); - - SMDS_NodeIteratorPtr ite= aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes(); - - bool isForward = (E.Orientation() == TopAbs_FORWARD); map params; + const SMDS_MeshNode * node; - while(ite->more()) + SMDS_NodeIteratorPtr nodeIt= aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes(); + while ( nodeIt->more() ) { - const SMDS_MeshNode * node = ite->next(); + node = nodeIt->next(); + if ( _quadraticMesh && SMESH_MesherHelper::IsMedium( node, SMDSAbs_Edge )) + continue; const SMDS_EdgePosition* epos = static_cast(node->GetPosition().get()); double param = epos->GetUParameter(); - params[param] = node; - } - if ( nbPoints != params.size()) - { - MESSAGE( "BAD NODE ON EDGE POSITIONS" ); - return false; + if ( !isForward ) param = -param; + if ( !params.insert( make_pair( param, node )).second ) + { + MESSAGE( "BAD NODE ON EDGE POSITIONS" ); + return false; + } } - mFirst = m; - // --- load 2D values into MEFISTO structure, // add IDNodes in mefistoToDS map - if (E.Orientation() == TopAbs_FORWARD) - { - gp_Pnt2d p = C2d->Value(f).XY().Multiplied( scale ); // first point = Vertex Forward - if ( fixCommonVertexUV( p, VFirst, W, myOuterWire, F, VWMap, aMesh )) - myNodesOnCommonV.push_back( idFirst ); - uvslf[m].x = p.X(); - uvslf[m].y = p.Y(); - mefistoToDS[m + 1] = idFirst; - //MESSAGE(" "<::iterator itp = params.begin(); - for (int i = 1; i <= nbPoints; i++) // nbPoints internal - { - double param = (*itp).first; - gp_Pnt2d p = C2d->Value(param).XY().Multiplied( scale ); - uvslf[m].x = p.X(); - uvslf[m].y = p.Y(); - mefistoToDS[m + 1] = (*itp).second; - //MESSAGE(" "<Value( uFirst ).XY().Multiplied( scale ); + if ( fixCommonVertexUV( p, V, W, myOuterWire, F, VWMap, aMesh, _quadraticMesh )) + myNodesOnCommonV.push_back( idFirst ); + mOnVertex.push_back( m ); + uvslf[m].x = p.X(); + uvslf[m].y = p.Y(); + mefistoToDS[m + 1] = idFirst; + //MESSAGE(" "<::iterator u_n = params.begin(); + for ( int i = 0; u_n != params.end(); ++u_n, ++i ) { - gp_Pnt2d p = C2d->Value(l).XY().Multiplied( scale ); // last point = Vertex Reversed - if ( fixCommonVertexUV( p, VLast, W, myOuterWire, F, VWMap, aMesh )) - myNodesOnCommonV.push_back( idLast ); + u = isForward ? u_n->first : - u_n->first; + gp_Pnt2d p = C2d->Value( u ).XY().Multiplied( scale ); uvslf[m].x = p.X(); uvslf[m].y = p.Y(); - mefistoToDS[m + 1] = idLast; + mefistoToDS[m + 1] = u_n->second; //MESSAGE(" "<::reverse_iterator itp = params.rbegin(); - for (int i = nbPoints; i >= 1; i--) - { - double param = (*itp).first; - gp_Pnt2d p = C2d->Value(param).XY().Multiplied( scale ); - uvslf[m].x = p.X(); - uvslf[m].y = p.Y(); - mefistoToDS[m + 1] = (*itp).second; - //MESSAGE(" "< 0 ) - fixOverlappedLinkUV (uvslf[ mFirst - 1], - uvslf[ mFirst ], - uvslf[ mFirst + 1 ]); - } // for wexp - fixOverlappedLinkUV (uvslf[ m - 1], - uvslf[ mInit ], - uvslf[ mInit + 1 ]); + // prevent failure on overlapped adjacent links, + // check only links ending in vertex nodes + int mFirst = mOnVertex.front(), mLast = m - 1; + list< int >::iterator mIt = mOnVertex.begin(); + for ( ; mIt != mOnVertex.end(); ++mIt ) { + int i = *mIt; + int iB = i - 1, iA = i + 1; // indices Before and After + if ( iB < mFirst ) iB = mLast; + if ( iA > mLast ) iA = mFirst; + fixOverlappedLinkUV (uvslf[ iB ], uvslf[ i ], uvslf[ iA ]); + } return true; } @@ -713,6 +695,7 @@ void StdMeshers_MEFISTO_2D::StoreResult(SMESH_Mesh & aMesh, double scalex, double scaley) { SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); + int faceID = meshDS->ShapeToIndex( F ); Z n, m; Handle(Geom_Surface) S = BRep_Tool::Surface(F); @@ -726,50 +709,40 @@ void StdMeshers_MEFISTO_2D::StoreResult(SMESH_Mesh & aMesh, gp_Pnt P = S->Value(u, v); SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(node, F); + meshDS->SetNodeOnFace(node, faceID, u, v); //MESSAGE(P.X()<<" "<(node->GetPosition().get()); - fpos->SetUParameter(u); - fpos->SetVParameter(v); } } m = 0; - int mt = 0; - - //SCRUTE(faceIsForward); - for (n = 1; n <= nbt; n++) - { - int inode1 = nust[m++]; - int inode2 = nust[m++]; - int inode3 = nust[m++]; - const SMDS_MeshNode *n1, *n2, *n3; - n1 = mefistoToDS[inode1]; - n2 = mefistoToDS[inode2]; - n3 = mefistoToDS[inode3]; - //MESSAGE("-- "<AddFace(n1, n2, n3); + //elt = meshDS->AddFace(n1, n2, n3); + elt = myTool->AddFace(n1, n2, n3); else - elt = meshDS->AddFace(n1, n3, n2); + //elt = meshDS->AddFace(n1, n3, n2); + elt = myTool->AddFace(n1, n3, n2); - meshDS->SetMeshElementOnShape(elt, F); + meshDS->SetMeshElementOnShape(elt, faceID); m++; } - // remove bad elements build on vertices shared by wires + // remove bad elements built on vertices shared by wires list::iterator itN = myNodesOnCommonV.begin(); for ( ; itN != myNodesOnCommonV.end(); itN++ ) @@ -805,17 +778,14 @@ double StdMeshers_MEFISTO_2D::ComputeEdgeElementLength(SMESH_Mesh & aMesh, //MESSAGE("StdMeshers_MEFISTO_2D::ComputeEdgeElementLength"); // **** a mettre dans SMESH_2D_Algo ? - const TopoDS_Face & FF = TopoDS::Face(aShape); - bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD); - TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD)); + //const TopoDS_Face & FF = TopoDS::Face(aShape); + //bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD); + //TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD)); double meanElementLength = 100; double wireLength = 0; int wireElementsNumber = 0; - for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next()) - { - const TopoDS_Wire & W = TopoDS::Wire(exp.Current()); - for (TopExp_Explorer expe(W, TopAbs_EDGE); expe.More(); expe.Next()) + for (TopExp_Explorer expe(aShape, TopAbs_EDGE); expe.More(); expe.Next()) { const TopoDS_Edge & E = TopoDS::Edge(expe.Current()); int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes(); @@ -823,7 +793,6 @@ double StdMeshers_MEFISTO_2D::ComputeEdgeElementLength(SMESH_Mesh & aMesh, wireLength += length; wireElementsNumber += nb; } - } if (wireElementsNumber) meanElementLength = wireLength / wireElementsNumber; //SCRUTE(meanElementLength);