X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_MEFISTO_2D.cxx;h=34003c6c9df03a5557a69e43f5c0acf4499b868b;hp=ca7405ba22c0673f41e258e715648256ff5f4a89;hb=ae32dcd34f98b91cdb4f5800063a394feb0df408;hpb=c63ee099ad2b149bd70136839c973e8910137bc5 diff --git a/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx b/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx index ca7405ba2..34003c6c9 100644 --- a/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx +++ b/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx @@ -1,66 +1,76 @@ -// SMESH SMESH : implementaion of SMESH idl descriptions +// Copyright (C) 2007-2014 CEA/DEN, EDF R&D, OPEN CASCADE // -// Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, -// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS -// -// This library is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 2.1 of the License. -// -// This library is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// Lesser General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// 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.salome-platform.org/ or email : webmaster.salome@opencascade.com +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS // +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. // +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. // +// You should have received a copy of the GNU Lesser General Public +// 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.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +// SMESH SMESH : implementaion of SMESH idl descriptions // File : StdMeshers_MEFISTO_2D.cxx // Moved here from SMESH_MEFISTO_2D.cxx // Author : Paul RASCLE, EDF // Module : SMESH -// $Header$ - -using namespace std; +// #include "StdMeshers_MEFISTO_2D.hxx" + +#include "SMDS_EdgePosition.hxx" +#include "SMDS_MeshElement.hxx" +#include "SMDS_MeshNode.hxx" +#include "SMESH_Comment.hxx" #include "SMESH_Gen.hxx" #include "SMESH_Mesh.hxx" +#include "SMESH_MesherHelper.hxx" #include "SMESH_subMesh.hxx" - -#include "StdMeshers_MaxElementArea.hxx" +#include "StdMeshers_FaceSide.hxx" #include "StdMeshers_LengthFromEdges.hxx" +#include "StdMeshers_MaxElementArea.hxx" +#include "StdMeshers_ViscousLayers2D.hxx" + +#include "utilities.h" #include "Rn.h" #include "aptrte.h" -#include "SMDS_MeshElement.hxx" -#include "SMDS_MeshNode.hxx" -#include "SMDS_EdgePosition.hxx" -#include "SMDS_FacePosition.hxx" - -#include "utilities.h" - -#include -#include -#include -#include -#include -#include -#include -#include +#include #include -#include -#include -#include -#include +#include +#include +#include +#include +#include +#include +#include +#include #include #include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +#ifdef _DEBUG_ +//#define DUMP_POINTS // to print coordinates of MEFISTO input +#endif //============================================================================= /*! @@ -68,20 +78,21 @@ using namespace std; */ //============================================================================= -StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D(int hypId, int studyId, - SMESH_Gen * gen):SMESH_2D_Algo(hypId, studyId, gen) +StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D(int hypId, int studyId, SMESH_Gen * gen): + SMESH_2D_Algo(hypId, studyId, gen) { MESSAGE("StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D"); _name = "MEFISTO_2D"; _shapeType = (1 << TopAbs_FACE); _compatibleHypothesis.push_back("MaxElementArea"); _compatibleHypothesis.push_back("LengthFromEdges"); + _compatibleHypothesis.push_back("ViscousLayers2D"); _edgeLength = 0; _maxElementArea = 0; _hypMaxElementArea = NULL; _hypLengthFromEdges = NULL; - myTool = 0; + _helper = 0; } //============================================================================= @@ -102,75 +113,71 @@ StdMeshers_MEFISTO_2D::~StdMeshers_MEFISTO_2D() //============================================================================= bool StdMeshers_MEFISTO_2D::CheckHypothesis - (SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape, + (SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, SMESH_Hypothesis::Hypothesis_Status& aStatus) { - //MESSAGE("StdMeshers_MEFISTO_2D::CheckHypothesis"); + _hypMaxElementArea = NULL; + _hypLengthFromEdges = NULL; + _edgeLength = 0; + _maxElementArea = 0; - _hypMaxElementArea = NULL; - _hypLengthFromEdges = NULL; + if ( !error( StdMeshers_ViscousLayers2D::CheckHypothesis( aMesh, aShape, aStatus ))) + return false; - list ::const_iterator itl; - const SMESHDS_Hypothesis *theHyp; + list ::const_iterator itl; + const SMESHDS_Hypothesis *theHyp; - const list &hyps = GetUsedHypothesis(aMesh, aShape); - int nbHyp = hyps.size(); - if (!nbHyp) - { - aStatus = SMESH_Hypothesis::HYP_MISSING; - return false; // can't work with no hypothesis - } + const list &hyps = GetUsedHypothesis(aMesh, aShape); + int nbHyp = hyps.size(); + if (!nbHyp) + { + aStatus = SMESH_Hypothesis::HYP_OK; //SMESH_Hypothesis::HYP_MISSING; + return true; // (PAL13464) can work with no hypothesis, LengthFromEdges is default one + } - itl = hyps.begin(); - theHyp = (*itl); // use only the first hypothesis - - string hypName = theHyp->GetName(); - //int hypId = theHyp->GetID(); - //SCRUTE(hypName); - - bool isOk = false; - - if (hypName == "MaxElementArea") - { - _hypMaxElementArea = static_cast(theHyp); - ASSERT(_hypMaxElementArea); - _maxElementArea = _hypMaxElementArea->GetMaxArea(); - _edgeLength = 0; - isOk = true; - aStatus = SMESH_Hypothesis::HYP_OK; - } - - else if (hypName == "LengthFromEdges") - { - _hypLengthFromEdges = static_cast(theHyp); - ASSERT(_hypLengthFromEdges); - _edgeLength = 0; - _maxElementArea = 0; - isOk = true; - aStatus = SMESH_Hypothesis::HYP_OK; - } - else - aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE; - - if (isOk) - { - isOk = false; - if (_maxElementArea > 0) - { -// _edgeLength = 2 * sqrt(_maxElementArea); // triangles : minorant - _edgeLength = 2 * sqrt(_maxElementArea/sqrt(3.0)); - isOk = true; - } - else - isOk = (_hypLengthFromEdges != NULL); // **** check mode - if (!isOk) - aStatus = SMESH_Hypothesis::HYP_BAD_PARAMETER; - } - - //SCRUTE(_edgeLength); - //SCRUTE(_maxElementArea); - return isOk; + itl = hyps.begin(); + theHyp = (*itl); // use only the first hypothesis + + string hypName = theHyp->GetName(); + + bool isOk = false; + + if (hypName == "MaxElementArea") + { + _hypMaxElementArea = static_cast(theHyp); + ASSERT(_hypMaxElementArea); + _maxElementArea = _hypMaxElementArea->GetMaxArea(); + isOk = true; + aStatus = SMESH_Hypothesis::HYP_OK; + } + + else if (hypName == "LengthFromEdges") + { + _hypLengthFromEdges = static_cast(theHyp); + ASSERT(_hypLengthFromEdges); + isOk = true; + aStatus = SMESH_Hypothesis::HYP_OK; + } + else + aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE; + + if (isOk) + { + isOk = false; + if (_maxElementArea > 0) + { + //_edgeLength = 2 * sqrt(_maxElementArea); // triangles : minorant + _edgeLength = sqrt(2. * _maxElementArea/sqrt(3.0)); + isOk = true; + } + else + isOk = (_hypLengthFromEdges != NULL); // **** check mode + if (!isOk) + aStatus = SMESH_Hypothesis::HYP_BAD_PARAMETER; + } + + return isOk; } //============================================================================= @@ -183,21 +190,52 @@ bool StdMeshers_MEFISTO_2D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh { MESSAGE("StdMeshers_MEFISTO_2D::Compute"); - if (_hypLengthFromEdges) - _edgeLength = ComputeEdgeElementLength(aMesh, aShape); - - bool isOk = false; - //const SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); - //SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape); + TopoDS_Face F = TopoDS::Face(aShape.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)); + // helper builds quadratic mesh if necessary + SMESH_MesherHelper helper(aMesh); + _helper = &helper; + _quadraticMesh = _helper->IsQuadraticSubMesh(aShape); + const bool skipMediumNodes = _quadraticMesh; - 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 + // build viscous layers if required + SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F ); + if ( !proxyMesh ) + return false; + + // get all edges of a face + TError problem; + TWireVector wires = + StdMeshers_FaceSide::GetFaceWires( F, aMesh, skipMediumNodes, problem, proxyMesh ); + int nbWires = wires.size(); + if ( problem && !problem->IsOK() ) return error( problem ); + if ( nbWires == 0 ) return error( "Problem in StdMeshers_FaceSide::GetFaceWires()"); + if ( wires[0]->NbSegments() < 3 ) // ex: a circle with 2 segments + return error(COMPERR_BAD_INPUT_MESH, + SMESH_Comment("Too few segments: ")<NbSegments()); + + // compute average edge length + if (!_hypMaxElementArea) + { + _edgeLength = 0; + int nbSegments = 0; + for ( int iW = 0; iW < nbWires; ++iW ) + { + StdMeshers_FaceSidePtr wire = wires[ iW ]; + _edgeLength += wire->Length(); + nbSegments += wire->NbSegments(); + } + if ( nbSegments ) + _edgeLength /= nbSegments; + } + + if (/*_hypLengthFromEdges &&*/ _edgeLength < DBL_MIN ) + _edgeLength = 100; + + 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; @@ -206,101 +244,147 @@ bool StdMeshers_MEFISTO_2D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh 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 + Z nutysu = 1; // 1: il existe un fonction areteideale_() + // Z nutysu=0; // 0: on utilise aretmx + R aretmx = _edgeLength; // longueur max aretes future triangulation + if ( _hypMaxElementArea ) + aretmx *= 1.5; - nblf = NumberOfWires(F); + nblf = nbWires; 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; - - 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; - } + // count nb of input points + for ( int iW = 0; iW < nbWires; ++iW ) + { + nbpnt += wires[iW]->NbPoints() - 1; + nudslf[iw++] = nbpnt; } - // 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 ); - uvslf = new R2[nudslf[nblf]]; - int m = 0; double scalex, scaley; ComputeScaleOnFace(aMesh, F, scalex, scaley); - map mefistoToDS; // correspondence mefisto index--> points IDNodes - if ( !LoadPoints(aMesh, F, myOuterWire, uvslf, m, - mefistoToDS, scalex, scaley, VWMap) ) { - delete myTool; myTool = 0; - return false; - } + // correspondence mefisto index --> Nodes + vector< const SMDS_MeshNode*> mefistoToDS(nbpnt, (const SMDS_MeshNode*)0); + + bool isOk = false; - for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next()) + // fill input points UV + if ( LoadPoints(wires, uvslf, mefistoToDS, scalex, scaley) ) { - const TopoDS_Wire & W = TopoDS::Wire(exp.Current()); - if (!myOuterWire.IsSame(W)) + // Compute + aptrte(nutysu, aretmx, + nblf, nudslf, uvslf, nbpti, uvpti, nbst, uvst, nbt, nust, ierr); + + if (ierr == 0) { - if (! LoadPoints(aMesh, F, W, uvslf, m, - mefistoToDS, scalex, scaley, VWMap )) { - delete myTool; myTool = 0; - return false; - } + MESSAGE("... End Triangulation Generated Triangle Number " << nbt); + MESSAGE(" Node Number " << nbst); + StoreResult(nbst, uvst, nbt, nust, mefistoToDS, scalex, scaley); + isOk = true; + } + else + { + error(ierr,"Error in Triangulation (aptrte())"); } } + if (nudslf != NULL) delete[]nudslf; + if (uvslf != NULL) delete[]uvslf; + if (uvst != NULL) delete[]uvst; + if (nust != NULL) delete[]nust; - uvst = NULL; - nust = NULL; - aptrte(nutysu, aretmx, - nblf, nudslf, uvslf, nbpti, uvpti, nbst, uvst, nbt, nust, ierr); + return isOk; +} - 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; + +//============================================================================= +/*! + * + */ +//============================================================================= + +bool StdMeshers_MEFISTO_2D::Evaluate(SMESH_Mesh & aMesh, + const TopoDS_Shape & aShape, + MapShapeNbElems& aResMap) +{ + MESSAGE("StdMeshers_MEFISTO_2D::Evaluate"); + + TopoDS_Face F = TopoDS::Face(aShape.Oriented(TopAbs_FORWARD)); + + double aLen = 0.0; + int NbSeg = 0; + bool IsQuadratic = false; + bool IsFirst = true; + TopExp_Explorer exp(F,TopAbs_EDGE); + for(; exp.More(); exp.Next()) { + TopoDS_Edge E = TopoDS::Edge(exp.Current()); + MapShapeNbElemsItr anIt = aResMap.find( aMesh.GetSubMesh(E) ); + if( anIt == aResMap.end() ) continue; + std::vector aVec = (*anIt).second; + int nbe = Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]); + NbSeg += nbe; + if(IsFirst) { + IsQuadratic = ( aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge] ); + IsFirst = false; + } + double a,b; + TopLoc_Location L; + Handle(Geom_Curve) C = BRep_Tool::Curve(E,L,a,b); + gp_Pnt P1; + C->D0(a,P1); + double dp = (b-a)/nbe; + for(int i=1; i<=nbe; i++) { + gp_Pnt P2; + C->D0(a+i*dp,P2); + aLen += P1.Distance(P2); + P1 = P2; + } } - else - { - MESSAGE("Error in Triangulation"); - isOk = false; + if(NbSeg<1) { + std::vector aResVec(SMDSEntity_Last); + for(int i=SMDSEntity_Node; iGetComputeError(); + smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED, + "Submesh can not be evaluated",this)); + return 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; + aLen = aLen/NbSeg; // middle length - return isOk; + _edgeLength = Precision::Infinite(); + double tmpLength = Min( _edgeLength, aLen ); + + GProp_GProps G; + BRepGProp::SurfaceProperties(aShape,G); + double anArea = G.Mass(); + + int nbFaces = Precision::IsInfinite( tmpLength ) ? 0 : + (int)( anArea/(tmpLength*tmpLength*sqrt(3.)/4) ); + int nbNodes = (int) ( nbFaces*3 - (NbSeg-1)*2 ) / 6; + + std::vector aVec(SMDSEntity_Last); + for(int i=SMDSEntity_Node; i UList; @@ -383,7 +472,7 @@ static bool fixCommonVertexUV (gp_Pnt2d & theUV, // find edges of theW sharing theV // and find 2d normal to them at theV gp_Vec2d N(0.,0.); - TopoDS_Iterator itE( theW ); + TopoDS_Iterator itE( anInnerWire ); for ( ; itE.More(); itE.Next() ) { const TopoDS_Edge& E = TopoDS::Edge( itE.Value() ); @@ -401,7 +490,7 @@ static bool fixCommonVertexUV (gp_Pnt2d & theUV, gp_Vec2d d1; gp_Pnt2d p; C2d->D1( u, p, d1 ); - gp_Vec2d n( d1.Y(), -d1.X() ); + gp_Vec2d n( d1.Y() * theScaleX, -d1.X() * theScaleY); if ( E.Orientation() == TopAbs_REVERSED ) n.Reverse(); N += n.Normalized(); @@ -411,6 +500,7 @@ static bool fixCommonVertexUV (gp_Pnt2d & theUV, // define step size by which to move theUV gp_Pnt2d nextUV; // uv of next node on edge, most distant of the four + gp_Pnt2d thisUV( theUV.x, theUV.y ); double maxDist = -DBL_MAX; TopTools_ListIteratorOfListOfShape aEIt (EList); list< double >::iterator aUIt = UList.begin(); @@ -431,10 +521,10 @@ static bool fixCommonVertexUV (gp_Pnt2d & theUV, while ( nIt->more() ) { const SMDS_MeshNode* node = nIt->next(); // check if node is medium - if ( CreateQuadratic && SMESH_MesherHelper::IsMedium( node, SMDSAbs_Edge )) + if ( theCreateQuadratic && SMESH_MesherHelper::IsMedium( node, SMDSAbs_Edge )) continue; const SMDS_EdgePosition* epos = - static_cast(node->GetPosition().get()); + static_cast(node->GetPosition()); double u = epos->GetUParameter(); if ( u < umin ) umin = u; @@ -444,29 +534,33 @@ static bool fixCommonVertexUV (gp_Pnt2d & theUV, } bool isFirstCommon = ( *aUIt == f ); gp_Pnt2d uv = C2d->Value( isFirstCommon ? umin : umax ); - double dist = theUV.SquareDistance( uv ); + double dist = thisUV.SquareDistance( uv ); if ( dist > maxDist ) { maxDist = dist; nextUV = uv; } } R2 uv0, uv1, uv2; - uv0.x = theUV.X(); uv0.y = theUV.Y(); + uv0.x = thisUV.X(); uv0.y = thisUV.Y(); uv1.x = nextUV.X(); uv1.y = nextUV.Y(); - uv2.x = uv0.x; uv2.y = uv0.y; + uv2.x = thisUV.X(); uv2.y = thisUV.Y(); + + uv1.x *= theScaleX; uv1.y *= theScaleY; + if ( fixOverlappedLinkUV( uv0, uv1, uv2 )) { - double step = theUV.Distance( gp_Pnt2d( uv0.x, uv0.y )); + double step = thisUV.Distance( gp_Pnt2d( uv0.x, uv0.y )); // move theUV along the normal by the step N *= step; - MESSAGE("--fixCommonVertexUV move(" << theUV.X() << " " << theUV.Y() + MESSAGE("--fixCommonVertexUV move(" << theUV.x << " " << theUV.x << ") by (" << N.X() << " " << N.Y() << ")" << endl << "--- MAX DIST " << maxDist); - theUV.SetXY( theUV.XY() + N.XY() ); + theUV.x += N.X(); + theUV.y += N.Y(); return true; } @@ -479,113 +573,109 @@ static bool fixCommonVertexUV (gp_Pnt2d & theUV, */ //============================================================================= -bool StdMeshers_MEFISTO_2D::LoadPoints(SMESH_Mesh & aMesh, - const TopoDS_Face & FF, - const TopoDS_Wire & WW, - R2 * uvslf, - int & m, - map&mefistoToDS, - double scalex, - double scaley, - const TopTools_IndexedDataMapOfShapeListOfShape& VWMap) +bool StdMeshers_MEFISTO_2D::LoadPoints(TWireVector & wires, + R2 * uvslf, + vector& mefistoToDS, + double scalex, + double scaley) { - // MESSAGE("StdMeshers_MEFISTO_2D::LoadPoints"); - - TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD)); - - 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); wexp.More(); wexp.Next() ) + // to avoid passing same uv points for a vertex common to 2 wires + TopoDS_Face F; + TopTools_IndexedDataMapOfShapeListOfShape VWMap; + if ( wires.size() > 1 ) { - const TopoDS_Edge & E = wexp.Current(); - bool isForward = (E.Orientation() == TopAbs_FORWARD); - - // --- IDNodes of first and last Vertex + F = TopoDS::Face( _helper->GetSubShape() ); + TopExp::MapShapesAndAncestors( F, TopAbs_VERTEX, TopAbs_WIRE, VWMap ); + int nbVertices = 0; + for ( int iW = 0; iW < wires.size(); ++iW ) + nbVertices += wires[ iW ]->NbEdges(); + if ( nbVertices == VWMap.Extent() ) + VWMap.Clear(); // wires have no common vertices + } - TopoDS_Vertex VFirst, VLast, V; - TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l - V = isForward ? VFirst : VLast; + int m = 0; - ASSERT(!V.IsNull()); - SMDS_NodeIteratorPtr lid= aMesh.GetSubMesh(V)->GetSubMeshDS()->GetNodes(); - if ( !lid->more() ) { - MESSAGE (" NO NODE BUILT ON VERTEX "); - return false; + for ( int iW = 0; iW < wires.size(); ++iW ) + { + const vector& uvPtVec = wires[ iW ]->GetUVPtStruct(); + if ( uvPtVec.size() != wires[ iW ]->NbPoints() ) { + return error(COMPERR_BAD_INPUT_MESH,SMESH_Comment("Unexpected nb of points on wire ") + << iW << ": " << uvPtVec.size()<<" != "<NbPoints() + << ", probably because of invalid node parameters on geom edges"); + } + if ( m + uvPtVec.size()-1 > mefistoToDS.size() ) { + MESSAGE("Wrong mefistoToDS.size: "<next(); - - // --- edge internal IDNodes (relies on good order storage, not checked) - - map params; - const SMDS_MeshNode * node; - SMDS_NodeIteratorPtr nodeIt= aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes(); - while ( nodeIt->more() ) + list< int > mOnVertex; + vector::const_iterator uvPt = uvPtVec.begin(); + for ( ++uvPt; uvPt != uvPtVec.end(); ++uvPt ) { - 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(); - if ( !isForward ) param = -param; - if ( !params.insert( make_pair( param, node )).second ) + // bind mefisto ID to node + mefistoToDS[m] = uvPt->node; + // set UV + uvslf[m].x = uvPt->u * scalex; + uvslf[m].y = uvPt->v * scaley; + switch ( uvPt->node->GetPosition()->GetTypeOfPosition()) { - MESSAGE( "BAD NODE ON EDGE POSITIONS" ); - return false; + case SMDS_TOP_VERTEX: + mOnVertex.push_back( m ); + break; + case SMDS_TOP_EDGE: + // In order to detect degenerated faces easily, we replace + // nodes on a degenerated edge by node on the vertex of that edge + if ( _helper->IsDegenShape( uvPt->node->getshapeId() )) + { + int edgeID = uvPt->node->getshapeId(); + SMESH_subMesh* edgeSM = _helper->GetMesh()->GetSubMeshContaining( edgeID ); + SMESH_subMeshIteratorPtr smIt = edgeSM->getDependsOnIterator( /*includeSelf=*/0, + /*complexShapeFirst=*/0); + if ( smIt->more() ) + { + SMESH_subMesh* vertexSM = smIt->next(); + SMDS_NodeIteratorPtr nIt = vertexSM->GetSubMeshDS()->GetNodes(); + if ( nIt->more() ) + mefistoToDS[m] = nIt->next(); + } + } + break; + default:; } + m++; } - // --- load 2D values into MEFISTO structure, - // add IDNodes in mefistoToDS map - - double f, l, uFirst, u; - Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l); - uFirst = isForward ? f : l; - - // vertex node - gp_Pnt2d p = C2d->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 ) + int mFirst = mOnVertex.front(), mLast = m - 1; + list< int >::iterator mIt = mOnVertex.begin(); + for ( ; mIt != mOnVertex.end(); ++mIt) { - 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] = u_n->second; - //MESSAGE(" "<getshapeId(); + TopoDS_Vertex V = TopoDS::Vertex( _helper->GetMeshDS()->IndexToShape( vID )); + if ( fixCommonVertexUV( uvslf[m], V, F, VWMap, *_helper->GetMesh(), + scalex, scaley, _quadraticMesh )) { + myNodesOnCommonV.push_back( mefistoToDS[m] ); + continue; + } + } + // prevent failure on overlapped adjacent links, + // check only links ending in vertex nodes + int mB = m - 1, mA = m + 1; // indices Before and After + if ( mB < mFirst ) mB = mLast; + if ( mA > mLast ) mA = mFirst; + fixOverlappedLinkUV (uvslf[ mB ], uvslf[ m ], uvslf[ mA ]); } - } // for wexp - - // 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 ]); } +#ifdef DUMP_POINTS + cout << "MEFISTO INPUT************" << endl; + for ( int i =0; i < m; ++i ) + cout << i << ": \t" << uvslf[i].x << ", " << uvslf[i].y + << " Node " << mefistoToDS[i]->GetID()<< endl; +#endif + return true; } @@ -595,12 +685,12 @@ bool StdMeshers_MEFISTO_2D::LoadPoints(SMESH_Mesh & aMesh, */ //============================================================================= -void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh & aMesh, - const TopoDS_Face & aFace, double &scalex, double &scaley) +void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh & aMesh, + const TopoDS_Face & aFace, + double & scalex, + double & scaley) { - //MESSAGE("StdMeshers_MEFISTO_2D::ComputeScaleOnFace"); - TopoDS_Face F = TopoDS::Face(aFace.Oriented(TopAbs_FORWARD)); - TopoDS_Wire W = BRepTools::OuterWire(F); + TopoDS_Wire W = BRepTools::OuterWire(aFace); double xmin = 1.e300; // min & max of face 2D parametric coord. double xmax = -1.e300; @@ -615,7 +705,7 @@ void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh & aMesh, { const TopoDS_Edge & E = TopoDS::Edge( wexp.Current() ); double f, l; - Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l); + Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, aFace, f, l); if ( C2d.IsNull() ) continue; double du = (l - f) / double (nbp); for (int i = 0; i <= nbp; i++) @@ -642,7 +732,8 @@ void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh & aMesh, double xsize = xmax - xmin; double ysize = ymax - ymin; - Handle(Geom_Surface) S = BRep_Tool::Surface(F); // 3D surface + TopLoc_Location L; + Handle(Geom_Surface) S = BRep_Tool::Surface(aFace,L); // 3D surface double length_x = 0; double length_y = 0; @@ -682,64 +773,85 @@ void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh & aMesh, ASSERT(scaley); } +// namespace +// { +// bool isDegenTria( const SMDS_MeshNode * nn[3] ) +// { +// SMESH_TNodeXYZ p1( nn[0] ); +// SMESH_TNodeXYZ p2( nn[1] ); +// SMESH_TNodeXYZ p3( nn[2] ); +// gp_XYZ vec1 = p2 - p1; +// gp_XYZ vec2 = p3 - p1; +// gp_XYZ cross = vec1 ^ vec2; +// const double eps = 1e-100; +// return ( fabs( cross.X() ) < eps && +// fabs( cross.Y() ) < eps && +// fabs( cross.Z() ) < eps ); +// } +// } + //============================================================================= /*! * */ //============================================================================= -void StdMeshers_MEFISTO_2D::StoreResult(SMESH_Mesh & aMesh, - Z nbst, R2 * uvst, Z nbt, Z * nust, - const TopoDS_Face & F, bool faceIsForward, - map&mefistoToDS, +void StdMeshers_MEFISTO_2D::StoreResult(Z nbst, R2 * uvst, Z nbt, Z * nust, + vector< const SMDS_MeshNode*>&mefistoToDS, double scalex, double scaley) { - SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); - int faceID = meshDS->ShapeToIndex( F ); + _helper->SetElementsOnShape( true ); + + TopoDS_Face F = TopoDS::Face( _helper->GetSubShape() ); + Handle(Geom_Surface) S = BRep_Tool::Surface( F ); - Z n, m; - Handle(Geom_Surface) S = BRep_Tool::Surface(F); + //const size_t nbInputNodes = mefistoToDS.size(); - for (n = 0; n < nbst; n++) + Z n = mefistoToDS.size(); // nb input points + mefistoToDS.resize( nbst ); + for ( ; n < nbst; n++) { - if (mefistoToDS.find(n + 1) == mefistoToDS.end()) + if (!mefistoToDS[n]) { double u = uvst[n][0] / scalex; double v = uvst[n][1] / scaley; gp_Pnt P = S->Value(u, v); - SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(node, faceID, u, v); - - //MESSAGE(P.X()<<" "<AddNode( P.X(), P.Y(), P.Z(), 0, u, v ); } } - m = 0; + Z m = 0; // triangle points must be in trigonometric order if face is Forward // else they must be put clockwise - bool triangleIsWellOriented = faceIsForward; + int i1 = 1, i2 = 2; + if ( F.Orientation() != TopAbs_FORWARD ) + std::swap( i1, i2 ); + const SMDS_MeshNode * nn[3]; for (n = 1; n <= nbt; n++) { - const SMDS_MeshNode * n1 = mefistoToDS[ nust[m++] ]; - const SMDS_MeshNode * n2 = mefistoToDS[ nust[m++] ]; - const SMDS_MeshNode * n3 = mefistoToDS[ nust[m++] ]; - - SMDS_MeshElement * elt; - if (triangleIsWellOriented) - //elt = meshDS->AddFace(n1, n2, n3); - elt = myTool->AddFace(n1, n2, n3); - else - //elt = meshDS->AddFace(n1, n3, n2); - elt = myTool->AddFace(n1, n3, n2); - - meshDS->SetMeshElementOnShape(elt, faceID); + // const bool allNodesAreOld = ( nust[m + 0] <= nbInputNodes && + // nust[m + 1] <= nbInputNodes && + // nust[m + 2] <= nbInputNodes ); + nn[ 0 ] = mefistoToDS[ nust[m++] - 1 ]; + nn[ 1 ] = mefistoToDS[ nust[m++] - 1 ]; + nn[ 2 ] = mefistoToDS[ nust[m++] - 1 ]; m++; + + // avoid creating degenetrated faces + bool isDegen = ( _helper->HasDegeneratedEdges() && + ( nn[0] == nn[1] || nn[1] == nn[2] || nn[2] == nn[0] )); + + // It was an attemp to fix a problem of a zero area face whose all nodes + // are on one staight EDGE. But omitting this face makes a hole in the mesh :( + // if ( !isDegen && allNodesAreOld ) + // isDegen = isDegenTria( nn ); + + if ( !isDegen ) + _helper->AddFace( nn[0], nn[i1], nn[i2] ); } // remove bad elements built on vertices shared by wires @@ -759,86 +871,8 @@ void StdMeshers_MEFISTO_2D::StoreResult(SMESH_Mesh & aMesh, nbSame++; if (nbSame > 1) { MESSAGE( "RM bad element " << elem->GetID()); - meshDS->RemoveElement( elem ); + _helper->GetMeshDS()->RemoveElement( elem ); } } } - -} - -//============================================================================= -/*! - * - */ -//============================================================================= - -double StdMeshers_MEFISTO_2D::ComputeEdgeElementLength(SMESH_Mesh & aMesh, - const TopoDS_Shape & aShape) -{ - //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)); - - double meanElementLength = 100; - double wireLength = 0; - int wireElementsNumber = 0; - 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(); - double length = EdgeLength(E); - wireLength += length; - wireElementsNumber += nb; - } - if (wireElementsNumber) - meanElementLength = wireLength / wireElementsNumber; - //SCRUTE(meanElementLength); - return meanElementLength; -} - -//============================================================================= -/*! - * - */ -//============================================================================= - -ostream & StdMeshers_MEFISTO_2D::SaveTo(ostream & save) -{ - return save; -} - -//============================================================================= -/*! - * - */ -//============================================================================= - -istream & StdMeshers_MEFISTO_2D::LoadFrom(istream & load) -{ - return load; -} - -//============================================================================= -/*! - * - */ -//============================================================================= - -ostream & operator <<(ostream & save, StdMeshers_MEFISTO_2D & hyp) -{ - return hyp.SaveTo( save ); -} - -//============================================================================= -/*! - * - */ -//============================================================================= - -istream & operator >>(istream & load, StdMeshers_MEFISTO_2D & hyp) -{ - return hyp.LoadFrom( load ); }