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=b7da7cbe615ec8ffe2c8069932703e36e0cc6a67;hb=ae32dcd34f98b91cdb4f5800063a394feb0df408;hpb=e4737e85f0da6d3f90fd08f6be1c2825195fe16f diff --git a/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx b/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx index b7da7cbe6..34003c6c9 100644 --- a/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx +++ b/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx @@ -1,69 +1,76 @@ -// SMESH SMESH : implementaion of SMESH idl descriptions +// Copyright (C) 2007-2014 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 +// +// 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. // -// 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.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org +// 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; -#include -//#include +#ifdef _DEBUG_ +//#define DUMP_POINTS // to print coordinates of MEFISTO input +#endif //============================================================================= /*! @@ -71,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 = 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"); + _compatibleHypothesis.push_back("ViscousLayers2D"); + + _edgeLength = 0; + _maxElementArea = 0; + _hypMaxElementArea = NULL; + _hypLengthFromEdges = NULL; + _helper = 0; } //============================================================================= @@ -95,7 +103,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"); } //============================================================================= @@ -105,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(); + + 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; + } - 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; + return isOk; } //============================================================================= @@ -184,114 +188,203 @@ bool StdMeshers_MEFISTO_2D::CheckHypothesis bool StdMeshers_MEFISTO_2D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape) { - 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); - - 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; - } - } - - // 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)) - 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 )) - 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; - return isOk; + MESSAGE("StdMeshers_MEFISTO_2D::Compute"); + + TopoDS_Face F = TopoDS::Face(aShape.Oriented(TopAbs_FORWARD)); + + // helper builds quadratic mesh if necessary + SMESH_MesherHelper helper(aMesh); + _helper = &helper; + _quadraticMesh = _helper->IsQuadraticSubMesh(aShape); + const bool skipMediumNodes = _quadraticMesh; + + // 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; + 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 + if ( _hypMaxElementArea ) + aretmx *= 1.5; + + nblf = nbWires; + + nudslf = new Z[1 + nblf]; + nudslf[0] = 0; + int iw = 1; + int nbpnt = 0; + + // count nb of input points + for ( int iW = 0; iW < nbWires; ++iW ) + { + nbpnt += wires[iW]->NbPoints() - 1; + nudslf[iw++] = nbpnt; + } + + uvslf = new R2[nudslf[nblf]]; + + double scalex, scaley; + ComputeScaleOnFace(aMesh, F, scalex, scaley); + + // correspondence mefisto index --> Nodes + vector< const SMDS_MeshNode*> mefistoToDS(nbpnt, (const SMDS_MeshNode*)0); + + bool isOk = false; + + // fill input points UV + if ( LoadPoints(wires, uvslf, mefistoToDS, scalex, scaley) ) + { + // Compute + 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(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; + + return isOk; +} + + +//============================================================================= +/*! + * + */ +//============================================================================= + +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; + } + } + 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; + } + aLen = aLen/NbSeg; // middle length + + _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 " << 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 +430,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; @@ -346,24 +441,30 @@ static bool fixOverlappedLinkUV( R2& uv0, const R2& uv1, const R2& uv2 ) //purpose : //======================================================================= -static bool fixCommonVertexUV (gp_Pnt2d & theUV, +static bool fixCommonVertexUV (R2 & theUV, const TopoDS_Vertex& theV, - const TopoDS_Wire& theW, - const TopoDS_Wire& theOW, const TopoDS_Face& theF, const TopTools_IndexedDataMapOfShapeListOfShape & theVWMap, - SMESH_Mesh & theMesh) + SMESH_Mesh & theMesh, + const double theScaleX, + const double theScaleY, + const bool theCreateQuadratic) { - if( theW.IsSame( theOW ) || - !theVWMap.Contains( theV )) return false; + if( !theVWMap.Contains( theV )) return false; // check if there is another wire sharing theV const TopTools_ListOfShape& WList = theVWMap.FindFromKey( theV ); TopTools_ListIteratorOfListOfShape aWIt; + TopTools_MapOfShape aWires; for ( aWIt.Initialize( WList ); aWIt.More(); aWIt.Next() ) - if ( !theW.IsSame( aWIt.Value() )) - break; - if ( !aWIt.More() ) return false; + aWires.Add( aWIt.Value() ); + if ( aWires.Extent() < 2 ) return false; + + TopoDS_Shape anOuterWire = BRepTools::OuterWire(theF); + TopoDS_Shape anInnerWire; + for ( aWIt.Initialize( WList ); aWIt.More() && anInnerWire.IsNull(); aWIt.Next() ) + if ( !anOuterWire.IsSame( aWIt.Value() )) + anInnerWire = aWIt.Value(); TopTools_ListOfShape EList; list< double > UList; @@ -371,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() ); @@ -389,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(); @@ -399,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(); @@ -415,12 +517,14 @@ 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 ( 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; @@ -430,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; } @@ -465,144 +573,108 @@ 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"); - - SMDS_Mesh * meshDS = aMesh.GetMeshDS(); - - TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD)); - - int mInit = m, mFirst, iEdge; - 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++) + // 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(); - - // --- IDNodes of first and last Vertex - - TopoDS_Vertex VFirst, VLast; - TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l - - ASSERT(!VFirst.IsNull()); - SMDS_NodeIteratorPtr lid= - aMesh.GetSubMesh(VFirst)->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(); + 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 + } - bool isForward = (E.Orientation() == TopAbs_FORWARD); - map params; + int m = 0; - while(ite->more()) - { - const SMDS_MeshNode * node = ite->next(); - const SMDS_EdgePosition* epos = - static_cast(node->GetPosition().get()); - double param = epos->GetUParameter(); - params[param] = node; + 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 ( nbPoints != params.size()) - { - MESSAGE( "BAD NODE ON EDGE POSITIONS" ); - return false; + if ( m + uvPtVec.size()-1 > mefistoToDS.size() ) { + MESSAGE("Wrong mefistoToDS.size: "< mOnVertex; + vector::const_iterator uvPt = uvPtVec.begin(); + for ( ++uvPt; uvPt != uvPtVec.end(); ++uvPt ) { - 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 + // 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()) { - 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(" "<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++; } - else + + int mFirst = mOnVertex.front(), mLast = m - 1; + list< int >::iterator mIt = mOnVertex.begin(); + for ( ; mIt != mOnVertex.end(); ++mIt) { - 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 ); - uvslf[m].x = p.X(); - uvslf[m].y = p.Y(); - mefistoToDS[m + 1] = idLast; - //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(" "<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 ]); } - // prevent failure on overlapped adjacent links - if ( iEdge > 0 ) - fixOverlappedLinkUV (uvslf[ mFirst - 1], - uvslf[ mFirst ], - uvslf[ mFirst + 1 ]); - - } // for wexp - - fixOverlappedLinkUV (uvslf[ m - 1], - uvslf[ mInit ], - uvslf[ mInit + 1 ]); + } + +#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; } @@ -613,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; @@ -633,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++) @@ -660,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; @@ -700,76 +773,88 @@ 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(); + _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, F); - - //MESSAGE(P.X()<<" "<(node->GetPosition().get()); - fpos->SetUParameter(u); - fpos->SetVParameter(v); + mefistoToDS[n] = _helper->AddNode( P.X(), P.Y(), P.Z(), 0, u, v ); } } - m = 0; - int mt = 0; + Z m = 0; - //SCRUTE(faceIsForward); - for (n = 1; n <= nbt; n++) - { - int inode1 = nust[m++]; - int inode2 = nust[m++]; - int inode3 = nust[m++]; + // triangle points must be in trigonometric order if face is Forward + // else they must be put clockwise - const SMDS_MeshNode *n1, *n2, *n3; - n1 = mefistoToDS[inode1]; - n2 = mefistoToDS[inode2]; - n3 = mefistoToDS[inode3]; - //MESSAGE("-- "<HasDegeneratedEdges() && + ( nn[0] == nn[1] || nn[1] == nn[2] || nn[2] == nn[0] )); - SMDS_MeshElement * elt; - if (triangleIsWellOriented) - elt = meshDS->AddFace(n1, n2, n3); - else - elt = meshDS->AddFace(n1, n3, n2); + // 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 ); - meshDS->SetMeshElementOnShape(elt, F); - m++; + if ( !isDegen ) + _helper->AddFace( nn[0], nn[i1], nn[i2] ); } - // 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++ ) @@ -786,90 +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 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()) - { - 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 ); }