1 // SMESH SMESH : implementaion of SMESH idl descriptions
3 // Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
24 // File : StdMeshers_MEFISTO_2D.cxx
25 // Moved here from SMESH_MEFISTO_2D.cxx
26 // Author : Paul RASCLE, EDF
31 #include "StdMeshers_MEFISTO_2D.hxx"
32 #include "SMESH_Gen.hxx"
33 #include "SMESH_Mesh.hxx"
34 #include "SMESH_subMesh.hxx"
36 #include "StdMeshers_MaxElementArea.hxx"
37 #include "StdMeshers_LengthFromEdges.hxx"
42 #include "SMDS_MeshElement.hxx"
43 #include "SMDS_MeshNode.hxx"
44 #include "SMDS_EdgePosition.hxx"
45 #include "SMDS_FacePosition.hxx"
47 #include "utilities.h"
49 #include <TopoDS_Face.hxx>
50 #include <TopoDS_Edge.hxx>
51 #include <TopoDS_Shape.hxx>
52 #include <Geom_Surface.hxx>
53 #include <GeomAdaptor_Curve.hxx>
54 #include <Geom2d_Curve.hxx>
55 #include <gp_Pnt2d.hxx>
56 #include <BRep_Tool.hxx>
57 #include <BRepTools.hxx>
58 #include <BRepTools_WireExplorer.hxx>
59 #include <GCPnts_AbscissaPoint.hxx>
60 #include <GCPnts_UniformAbscissa.hxx>
61 #include <TColStd_ListIteratorOfListOfInteger.hxx>
62 #include <TopTools_ListIteratorOfListOfShape.hxx>
63 #include <TopTools_ListOfShape.hxx>
66 //#include <algorithm>
68 //=============================================================================
72 //=============================================================================
74 StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D(int hypId, int studyId,
75 SMESH_Gen * gen):SMESH_2D_Algo(hypId, studyId, gen)
77 MESSAGE("StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D");
79 // _shapeType = TopAbs_FACE;
80 _shapeType = (1 << TopAbs_FACE);
81 _compatibleHypothesis.push_back("MaxElementArea");
82 _compatibleHypothesis.push_back("LengthFromEdges");
86 _hypMaxElementArea = NULL;
87 _hypLengthFromEdges = NULL;
90 //=============================================================================
94 //=============================================================================
96 StdMeshers_MEFISTO_2D::~StdMeshers_MEFISTO_2D()
98 MESSAGE("StdMeshers_MEFISTO_2D::~StdMeshers_MEFISTO_2D");
101 //=============================================================================
105 //=============================================================================
107 bool StdMeshers_MEFISTO_2D::CheckHypothesis
109 const TopoDS_Shape& aShape,
110 SMESH_Hypothesis::Hypothesis_Status& aStatus)
112 //MESSAGE("StdMeshers_MEFISTO_2D::CheckHypothesis");
114 _hypMaxElementArea = NULL;
115 _hypLengthFromEdges = NULL;
117 list <const SMESHDS_Hypothesis * >::const_iterator itl;
118 const SMESHDS_Hypothesis *theHyp;
120 const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
121 int nbHyp = hyps.size();
124 aStatus = SMESH_Hypothesis::HYP_MISSING;
125 return false; // can't work with no hypothesis
129 theHyp = (*itl); // use only the first hypothesis
131 string hypName = theHyp->GetName();
132 //int hypId = theHyp->GetID();
137 if (hypName == "MaxElementArea")
139 _hypMaxElementArea = static_cast<const StdMeshers_MaxElementArea *>(theHyp);
140 ASSERT(_hypMaxElementArea);
141 _maxElementArea = _hypMaxElementArea->GetMaxArea();
144 aStatus = SMESH_Hypothesis::HYP_OK;
147 else if (hypName == "LengthFromEdges")
149 _hypLengthFromEdges = static_cast<const StdMeshers_LengthFromEdges *>(theHyp);
150 ASSERT(_hypLengthFromEdges);
154 aStatus = SMESH_Hypothesis::HYP_OK;
157 aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
162 if (_maxElementArea > 0)
164 // _edgeLength = 2 * sqrt(_maxElementArea); // triangles : minorant
165 _edgeLength = 2 * sqrt(_maxElementArea/sqrt(3.0));
169 isOk = (_hypLengthFromEdges != NULL); // **** check mode
171 aStatus = SMESH_Hypothesis::HYP_BAD_PARAMETER;
174 //SCRUTE(_edgeLength);
175 //SCRUTE(_maxElementArea);
179 //=============================================================================
183 //=============================================================================
185 bool StdMeshers_MEFISTO_2D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape)
187 MESSAGE("StdMeshers_MEFISTO_2D::Compute");
189 if (_hypLengthFromEdges)
190 _edgeLength = ComputeEdgeElementLength(aMesh, aShape);
193 //const SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
194 //SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape);
196 const TopoDS_Face & FF = TopoDS::Face(aShape);
197 bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
198 TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
200 Z nblf; //nombre de lignes fermees (enveloppe en tete)
201 Z *nudslf = NULL; //numero du dernier sommet de chaque ligne fermee
203 Z nbpti = 0; //nombre points internes futurs sommets de la triangulation
212 Z nutysu = 1; // 1: il existe un fonction areteideale_()
213 // Z nutysu=0; // 0: on utilise aretmx
214 R aretmx = _edgeLength; // longueur max aretes future triangulation
216 nblf = NumberOfWires(F);
218 nudslf = new Z[1 + nblf];
223 bool QuadMode = true;
225 myTool = new StdMeshers_Helper(aMesh);
226 myCreateQuadratic = myTool->IsQuadraticSubMesh(aShape,QuadMode);
228 myOuterWire = BRepTools::OuterWire(F);
229 nbpnt += NumberOfPoints(aMesh, myOuterWire);
230 if ( nbpnt < 3 ) // ex: a circle with 2 segments
232 nudslf[iw++] = nbpnt;
234 for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next()) {
235 const TopoDS_Wire & W = TopoDS::Wire(exp.Current());
236 if (!myOuterWire.IsSame(W)) {
237 nbpnt += NumberOfPoints(aMesh, W);
238 nudslf[iw++] = nbpnt;
242 // avoid passing same uv points for a vertex common to 2 wires
243 TopTools_IndexedDataMapOfShapeListOfShape VWMap;
244 if ( iw - 1 > 1 ) // nbofWires > 1
245 TopExp::MapShapesAndAncestors( F , TopAbs_VERTEX, TopAbs_WIRE, VWMap );
247 uvslf = new R2[nudslf[nblf]];
250 double scalex, scaley;
251 ComputeScaleOnFace(aMesh, F, scalex, scaley);
253 map<int, const SMDS_MeshNode*> mefistoToDS; // correspondence mefisto index--> points IDNodes
254 if ( !LoadPoints(aMesh, F, myOuterWire, uvslf, m,
255 mefistoToDS, scalex, scaley, VWMap) )
258 for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next())
260 const TopoDS_Wire & W = TopoDS::Wire(exp.Current());
261 if (!myOuterWire.IsSame(W))
263 if (! LoadPoints(aMesh, F, W, uvslf, m,
264 mefistoToDS, scalex, scaley, VWMap ))
271 aptrte(nutysu, aretmx,
272 nblf, nudslf, uvslf, nbpti, uvpti, nbst, uvst, nbt, nust, ierr);
276 MESSAGE("... End Triangulation Generated Triangle Number " << nbt);
277 MESSAGE(" Node Number " << nbst);
278 StoreResult(aMesh, nbst, uvst, nbt, nust, F,
279 faceIsForward, mefistoToDS, scalex, scaley);
284 MESSAGE("Error in Triangulation");
298 //=======================================================================
299 //function : fixOverlappedLinkUV
300 //purpose : prevent failure due to overlapped adjacent links
301 //=======================================================================
303 static bool fixOverlappedLinkUV( R2& uv0, const R2& uv1, const R2& uv2 )
305 gp_XY v1( uv0.x - uv1.x, uv0.y - uv1.y );
306 gp_XY v2( uv2.x - uv1.x, uv2.y - uv1.y );
308 double tol2 = DBL_MIN * DBL_MIN;
309 double sqMod1 = v1.SquareModulus();
310 if ( sqMod1 <= tol2 ) return false;
311 double sqMod2 = v2.SquareModulus();
312 if ( sqMod2 <= tol2 ) return false;
316 // check sinus >= 1.e-3
317 const double minSin = 1.e-3;
318 if ( dot > 0 && 1 - dot * dot / ( sqMod1 * sqMod2 ) < minSin * minSin ) {
319 MESSAGE(" ___ FIX UV ____" << uv0.x << " " << uv0.y);
320 v1.SetCoord( -v1.Y(), v1.X() );
321 double delta = sqrt( sqMod1 ) * minSin;
330 // MESSAGE(" -> " << uv0.x << " " << uv0.y << " ");
331 // MESSAGE("v1( " << v1.X() << " " << v1.Y() << " ) " <<
332 // "v2( " << v2.X() << " " << v2.Y() << " ) ");
333 // MESSAGE("SIN: " << sqrt(1 - dot * dot / (sqMod1 * sqMod2)));
334 // v1.SetCoord( uv0.x - uv1.x, uv0.y - uv1.y );
335 // v2.SetCoord( uv2.x - uv1.x, uv2.y - uv1.y );
336 // gp_XY v3( uv2.x - uv0.x, uv2.y - uv0.y );
337 // sqMod1 = v1.SquareModulus();
338 // sqMod2 = v2.SquareModulus();
340 // double sin = sqrt(1 - dot * dot / (sqMod1 * sqMod2));
341 // MESSAGE("NEW SIN: " << sin);
347 //=======================================================================
348 //function : fixCommonVertexUV
350 //=======================================================================
352 static bool fixCommonVertexUV (gp_Pnt2d & theUV,
353 const TopoDS_Vertex& theV,
354 const TopoDS_Wire& theW,
355 const TopoDS_Wire& theOW,
356 const TopoDS_Face& theF,
357 const TopTools_IndexedDataMapOfShapeListOfShape & theVWMap,
358 SMESH_Mesh & theMesh,
359 bool CreateQuadratic)
361 if( theW.IsSame( theOW ) ||
362 !theVWMap.Contains( theV )) return false;
364 // check if there is another wire sharing theV
365 const TopTools_ListOfShape& WList = theVWMap.FindFromKey( theV );
366 TopTools_ListIteratorOfListOfShape aWIt;
367 for ( aWIt.Initialize( WList ); aWIt.More(); aWIt.Next() )
368 if ( !theW.IsSame( aWIt.Value() ))
370 if ( !aWIt.More() ) return false;
372 TopTools_ListOfShape EList;
373 list< double > UList;
375 // find edges of theW sharing theV
376 // and find 2d normal to them at theV
378 TopoDS_Iterator itE( theW );
379 for ( ; itE.More(); itE.Next() )
381 const TopoDS_Edge& E = TopoDS::Edge( itE.Value() );
382 TopoDS_Iterator itV( E );
383 for ( ; itV.More(); itV.Next() )
385 const TopoDS_Vertex & V = TopoDS::Vertex( itV.Value() );
386 if ( !V.IsSame( theV ))
389 Standard_Real u = BRep_Tool::Parameter( V, E );
390 UList.push_back( u );
392 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, theF, f, l);
396 gp_Vec2d n( d1.Y(), -d1.X() );
397 if ( E.Orientation() == TopAbs_REVERSED )
403 // define step size by which to move theUV
405 gp_Pnt2d nextUV; // uv of next node on edge, most distant of the four
406 double maxDist = -DBL_MAX;
407 TopTools_ListIteratorOfListOfShape aEIt (EList);
408 list< double >::iterator aUIt = UList.begin();
409 for ( ; aEIt.More(); aEIt.Next(), aUIt++ )
411 const TopoDS_Edge& E = TopoDS::Edge( aEIt.Value() );
413 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, theF, f, l);
415 double umin = DBL_MAX, umax = -DBL_MAX;
416 SMDS_NodeIteratorPtr nIt = theMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes();
417 if ( !nIt->more() ) // no nodes on edge, only on vertices
423 while ( nIt->more() ) {
424 const SMDS_MeshNode* node = nIt->next();
425 if(CreateQuadratic) {
426 // check if node is medium
427 bool IsMedium = false;
428 SMDS_ElemIteratorPtr itn = node->GetInverseElementIterator();
429 while (itn->more()) {
430 const SMDS_MeshElement* elem = itn->next();
431 if ( elem->GetType() != SMDSAbs_Edge )
433 if(elem->IsMediumNode(node)) {
441 const SMDS_EdgePosition* epos =
442 static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
443 double u = epos->GetUParameter();
450 bool isFirstCommon = ( *aUIt == f );
451 gp_Pnt2d uv = C2d->Value( isFirstCommon ? umin : umax );
452 double dist = theUV.SquareDistance( uv );
453 if ( dist > maxDist ) {
459 uv0.x = theUV.X(); uv0.y = theUV.Y();
460 uv1.x = nextUV.X(); uv1.y = nextUV.Y();
461 uv2.x = uv0.x; uv2.y = uv0.y;
462 if ( fixOverlappedLinkUV( uv0, uv1, uv2 ))
464 double step = theUV.Distance( gp_Pnt2d( uv0.x, uv0.y ));
466 // move theUV along the normal by the step
470 MESSAGE("--fixCommonVertexUV move(" << theUV.X() << " " << theUV.Y()
471 << ") by (" << N.X() << " " << N.Y() << ")"
472 << endl << "--- MAX DIST " << maxDist);
474 theUV.SetXY( theUV.XY() + N.XY() );
481 //=============================================================================
485 //=============================================================================
487 bool StdMeshers_MEFISTO_2D::LoadPoints(SMESH_Mesh & aMesh,
488 const TopoDS_Face & FF,
489 const TopoDS_Wire & WW,
492 map<int, const SMDS_MeshNode*>&mefistoToDS,
495 const TopTools_IndexedDataMapOfShapeListOfShape& VWMap)
497 // MESSAGE("StdMeshers_MEFISTO_2D::LoadPoints");
499 //SMDS_Mesh * meshDS = aMesh.GetMeshDS();
501 TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
503 int mInit = m, mFirst, iEdge;
504 gp_XY scale( scalex, scaley );
506 TopoDS_Wire W = TopoDS::Wire(WW.Oriented(TopAbs_FORWARD));
507 BRepTools_WireExplorer wexp(W, F);
508 for (wexp.Init(W, F), iEdge = 0; wexp.More(); wexp.Next(), iEdge++) {
509 const TopoDS_Edge & E = wexp.Current();
511 // --- IDNodes of first and last Vertex
513 TopoDS_Vertex VFirst, VLast;
514 TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
516 ASSERT(!VFirst.IsNull());
517 SMDS_NodeIteratorPtr lid =
518 aMesh.GetSubMesh(VFirst)->GetSubMeshDS()->GetNodes();
519 if ( !lid->more() ) {
520 MESSAGE (" NO NODE BUILT ON VERTEX ");
523 const SMDS_MeshNode* idFirst = lid->next();
525 ASSERT(!VLast.IsNull());
526 lid = aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes();
527 if ( !lid->more() ) {
528 MESSAGE (" NO NODE BUILT ON VERTEX ");
531 const SMDS_MeshNode* idLast = lid->next();
533 // --- edge internal IDNodes (relies on good order storage, not checked)
535 // if(myCreateQuadratic) {
536 // fill myNLinkNodeMap
537 // SMDS_ElemIteratorPtr iter = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetElements();
538 // while(iter->more()) {
539 // const SMDS_MeshElement* elem = iter->next();
540 // SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
541 // const SMDS_MeshNode* n1 = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
542 // const SMDS_MeshNode* n2 = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
543 // const SMDS_MeshNode* n3 = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
544 // NLink link(( n1 < n2 ? n1 : n2 ), ( n1 < n2 ? n2 : n1 ));
545 // myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n3));
546 // myNLinkNodeMap[link] = n3;
550 int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
553 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
555 SMDS_NodeIteratorPtr ite= aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes();
557 //bool isForward = (E.Orientation() == TopAbs_FORWARD);
558 map<double, const SMDS_MeshNode*> params;
560 if(!myCreateQuadratic) {
562 const SMDS_MeshNode * node = ite->next();
563 const SMDS_EdgePosition* epos =
564 static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
565 double param = epos->GetUParameter();
566 params[param] = node;
570 nbPoints = nbPoints/2;
572 const SMDS_MeshNode* node = ite->next();
573 // check if node is medium
574 bool IsMedium = false;
575 SMDS_ElemIteratorPtr itn = node->GetInverseElementIterator();
576 while (itn->more()) {
577 const SMDS_MeshElement* elem = itn->next();
578 if ( elem->GetType() != SMDSAbs_Edge )
580 if(elem->IsMediumNode(node)) {
587 const SMDS_EdgePosition* epos =
588 static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
589 double param = epos->GetUParameter();
590 params[param] = node;
594 if ( nbPoints != params.size()) {
595 MESSAGE( "BAD NODE ON EDGE POSITIONS" );
601 // --- load 2D values into MEFISTO structure,
602 // add IDNodes in mefistoToDS map
603 if (E.Orientation() == TopAbs_FORWARD) {
604 gp_Pnt2d p = C2d->Value(f).XY().Multiplied( scale ); // first point = Vertex Forward
605 if ( fixCommonVertexUV( p, VFirst, W, myOuterWire, F, VWMap, aMesh, myCreateQuadratic ))
606 myNodesOnCommonV.push_back( idFirst );
609 mefistoToDS[m + 1] = idFirst;
610 //MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
611 //MESSAGE("__ f "<<f<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
613 map<double, const SMDS_MeshNode*>::iterator itp = params.begin();
614 for (int i = 1; i <= nbPoints; i++) { // nbPoints internal
615 double param = (*itp).first;
616 gp_Pnt2d p = C2d->Value(param).XY().Multiplied( scale );
619 mefistoToDS[m + 1] = (*itp).second;
620 //MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
621 //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
627 gp_Pnt2d p = C2d->Value(l).XY().Multiplied( scale ); // last point = Vertex Reversed
628 if ( fixCommonVertexUV( p, VLast, W, myOuterWire, F, VWMap, aMesh, myCreateQuadratic ))
629 myNodesOnCommonV.push_back( idLast );
632 mefistoToDS[m + 1] = idLast;
633 //MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
634 //MESSAGE("__ l "<<l<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
636 map<double, const SMDS_MeshNode*>::reverse_iterator itp = params.rbegin();
637 for (int i = nbPoints; i >= 1; i--)
639 double param = (*itp).first;
640 gp_Pnt2d p = C2d->Value(param).XY().Multiplied( scale );
643 mefistoToDS[m + 1] = (*itp).second;
644 //MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
645 //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
650 // prevent failure on overlapped adjacent links
652 fixOverlappedLinkUV (uvslf[ mFirst - 1],
654 uvslf[ mFirst + 1 ]);
658 fixOverlappedLinkUV (uvslf[ m - 1],
665 //=============================================================================
669 //=============================================================================
671 void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh & aMesh,
672 const TopoDS_Face & aFace, double &scalex, double &scaley)
674 //MESSAGE("StdMeshers_MEFISTO_2D::ComputeScaleOnFace");
675 TopoDS_Face F = TopoDS::Face(aFace.Oriented(TopAbs_FORWARD));
676 TopoDS_Wire W = BRepTools::OuterWire(F);
678 double xmin = 1.e300; // min & max of face 2D parametric coord.
679 double xmax = -1.e300;
680 double ymin = 1.e300;
681 double ymax = -1.e300;
686 TopExp_Explorer wexp(W, TopAbs_EDGE);
687 for ( ; wexp.More(); wexp.Next())
689 const TopoDS_Edge & E = TopoDS::Edge( wexp.Current() );
691 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
692 if ( C2d.IsNull() ) continue;
693 double du = (l - f) / double (nbp);
694 for (int i = 0; i <= nbp; i++)
696 double param = f + double (i) * du;
697 gp_Pnt2d p = C2d->Value(param);
706 // MESSAGE(" "<< f<<" "<<l<<" "<<param<<" "<<xmin<<" "<<xmax<<" "<<ymin<<" "<<ymax);
713 double xmoy = (xmax + xmin) / 2.;
714 double ymoy = (ymax + ymin) / 2.;
715 double xsize = xmax - xmin;
716 double ysize = ymax - ymin;
718 Handle(Geom_Surface) S = BRep_Tool::Surface(F); // 3D surface
722 gp_Pnt PX0 = S->Value(xmin, ymoy);
723 gp_Pnt PY0 = S->Value(xmoy, ymin);
724 double dx = xsize / double (nbp);
725 double dy = ysize / double (nbp);
726 for (int i = 1; i <= nbp; i++)
728 double x = xmin + double (i) * dx;
729 gp_Pnt PX = S->Value(x, ymoy);
730 double y = ymin + double (i) * dy;
731 gp_Pnt PY = S->Value(xmoy, y);
732 length_x += PX.Distance(PX0);
733 length_y += PY.Distance(PY0);
737 scalex = length_x / xsize;
738 scaley = length_y / ysize;
741 double xyratio = xsize*scalex/(ysize*scaley);
742 const double maxratio = 1.e2;
744 if (xyratio > maxratio) {
746 scaley *= xyratio / maxratio;
749 else if (xyratio < 1./maxratio) {
751 scalex *= 1 / xyratio / maxratio;
758 //=============================================================================
762 //=============================================================================
764 void StdMeshers_MEFISTO_2D::StoreResult(SMESH_Mesh & aMesh,
765 Z nbst, R2 * uvst, Z nbt, Z * nust,
766 const TopoDS_Face & F, bool faceIsForward,
767 map<int, const SMDS_MeshNode*>&mefistoToDS,
768 double scalex, double scaley)
770 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
771 int faceID = meshDS->ShapeToIndex( F );
774 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
776 for (n = 0; n < nbst; n++)
778 if (mefistoToDS.find(n + 1) == mefistoToDS.end())
780 double u = uvst[n][0] / scalex;
781 double v = uvst[n][1] / scaley;
782 gp_Pnt P = S->Value(u, v);
784 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
785 meshDS->SetNodeOnFace(node, faceID, u, v);
787 //MESSAGE(P.X()<<" "<<P.Y()<<" "<<P.Z());
788 mefistoToDS[n + 1] = node;
789 //MESSAGE("NEW: "<<n<<" "<<mefistoToDS[n+1]);
796 //SCRUTE(faceIsForward);
797 for (n = 1; n <= nbt; n++)
799 int inode1 = nust[m++];
800 int inode2 = nust[m++];
801 int inode3 = nust[m++];
803 const SMDS_MeshNode *n1, *n2, *n3;
804 n1 = mefistoToDS[inode1];
805 n2 = mefistoToDS[inode2];
806 n3 = mefistoToDS[inode3];
807 //MESSAGE("-- "<<inode1<<" "<<inode2<<" "<<inode3);
809 // triangle points must be in trigonometric order if face is Forward
810 // else they must be put clockwise
812 bool triangleIsWellOriented = faceIsForward;
814 SMDS_MeshElement * elt;
815 if (triangleIsWellOriented)
816 //elt = meshDS->AddFace(n1, n2, n3);
817 elt = myTool->AddFace(n1, n2, n3);
819 //elt = meshDS->AddFace(n1, n3, n2);
820 elt = myTool->AddFace(n1, n3, n2);
822 meshDS->SetMeshElementOnShape(elt, faceID);
826 // remove bad elements build on vertices shared by wires
828 list<const SMDS_MeshNode*>::iterator itN = myNodesOnCommonV.begin();
829 for ( ; itN != myNodesOnCommonV.end(); itN++ )
831 const SMDS_MeshNode* node = *itN;
832 SMDS_ElemIteratorPtr invElemIt = node->GetInverseElementIterator();
833 while ( invElemIt->more() )
835 const SMDS_MeshElement* elem = invElemIt->next();
836 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
838 while ( itN->more() )
839 if ( itN->next() == node)
842 MESSAGE( "RM bad element " << elem->GetID());
843 meshDS->RemoveElement( elem );
850 //=============================================================================
854 //=============================================================================
856 double StdMeshers_MEFISTO_2D::ComputeEdgeElementLength(SMESH_Mesh & aMesh,
857 const TopoDS_Shape & aShape)
859 //MESSAGE("StdMeshers_MEFISTO_2D::ComputeEdgeElementLength");
860 // **** a mettre dans SMESH_2D_Algo ?
862 const TopoDS_Face & FF = TopoDS::Face(aShape);
863 //bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
864 TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
866 double meanElementLength = 100;
867 double wireLength = 0;
868 int wireElementsNumber = 0;
869 for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next())
871 const TopoDS_Wire & W = TopoDS::Wire(exp.Current());
872 for (TopExp_Explorer expe(W, TopAbs_EDGE); expe.More(); expe.Next())
874 const TopoDS_Edge & E = TopoDS::Edge(expe.Current());
875 int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
876 double length = EdgeLength(E);
877 wireLength += length;
878 wireElementsNumber += nb;
881 if (wireElementsNumber)
882 meanElementLength = wireLength / wireElementsNumber;
883 //SCRUTE(meanElementLength);
884 return meanElementLength;
887 //=============================================================================
891 //=============================================================================
893 ostream & StdMeshers_MEFISTO_2D::SaveTo(ostream & save)
898 //=============================================================================
902 //=============================================================================
904 istream & StdMeshers_MEFISTO_2D::LoadFrom(istream & load)
909 //=============================================================================
913 //=============================================================================
915 ostream & operator <<(ostream & save, StdMeshers_MEFISTO_2D & hyp)
917 return hyp.SaveTo( save );
920 //=============================================================================
924 //=============================================================================
926 istream & operator >>(istream & load, StdMeshers_MEFISTO_2D & hyp)
928 return hyp.LoadFrom( load );