2 //=============================================================================
3 // File : SMESH_MEFISTO_2D.cxx
4 // Created : sam mai 18 08:10:55 CEST 2002
5 // Author : Paul RASCLE, EDF
7 // Copyright : EDF 2002
9 //=============================================================================
12 #include "SMESH_MEFISTO_2D.hxx"
13 #include "SMESH_Gen.hxx"
14 #include "SMESH_Mesh.hxx"
16 #include "SMESH_MaxElementArea.hxx"
17 #include "SMESH_LengthFromEdges.hxx"
22 #include "SMESHDS_ListOfPtrHypothesis.hxx"
23 #include "SMESHDS_ListIteratorOfListOfPtrHypothesis.hxx"
24 #include "SMDS_MeshElement.hxx"
25 #include "SMDS_MeshNode.hxx"
26 #include "SMDS_EdgePosition.hxx"
27 #include "SMDS_FacePosition.hxx"
29 #include "utilities.h"
31 #include <TopoDS_Face.hxx>
32 #include <TopoDS_Edge.hxx>
33 #include <TopoDS_Shape.hxx>
34 #include <Geom_Surface.hxx>
35 #include <GeomAdaptor_Curve.hxx>
36 #include <Geom2d_Curve.hxx>
37 #include <gp_Pnt2d.hxx>
38 #include <BRep_Tool.hxx>
39 #include <BRepTools.hxx>
40 #include <BRepTools_WireExplorer.hxx>
41 #include <GCPnts_AbscissaPoint.hxx>
42 #include <GCPnts_UniformAbscissa.hxx>
43 #include <TColStd_ListIteratorOfListOfInteger.hxx>
48 //=============================================================================
52 //=============================================================================
54 SMESH_MEFISTO_2D::SMESH_MEFISTO_2D(int hypId, int studyId, SMESH_Gen* gen)
55 : SMESH_2D_Algo(hypId, studyId, gen)
57 MESSAGE("SMESH_MEFISTO_2D::SMESH_MEFISTO_2D");
59 // _shapeType = TopAbs_FACE;
60 _shapeType = (1<<TopAbs_FACE);
61 _compatibleHypothesis.push_back("MaxElementArea");
62 _compatibleHypothesis.push_back("LengthFromEdges");
66 _hypMaxElementArea = NULL;
67 _hypLengthFromEdges = NULL;
70 //=============================================================================
74 //=============================================================================
76 SMESH_MEFISTO_2D::~SMESH_MEFISTO_2D()
78 MESSAGE("SMESH_MEFISTO_2D::~SMESH_MEFISTO_2D");
81 //=============================================================================
85 //=============================================================================
87 bool SMESH_MEFISTO_2D::CheckHypothesis(SMESH_Mesh& aMesh,
88 const TopoDS_Shape& aShape)
90 //MESSAGE("SMESH_MEFISTO_2D::CheckHypothesis");
92 _hypMaxElementArea = NULL;
93 _hypLengthFromEdges = NULL;
95 list<SMESHDS_Hypothesis*>::const_iterator itl;
96 SMESHDS_Hypothesis* theHyp;
98 const list<SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
99 int nbHyp = hyps.size();
100 if (nbHyp != 1) return false; // only one compatible hypothesis allowed
105 string hypName = theHyp->GetName();
106 int hypId = theHyp->GetID();
111 if (hypName == "MaxElementArea")
113 _hypMaxElementArea = dynamic_cast<SMESH_MaxElementArea*> (theHyp);
114 ASSERT(_hypMaxElementArea);
115 _maxElementArea = _hypMaxElementArea->GetMaxArea();
120 if (hypName == "LengthFromEdges")
122 _hypLengthFromEdges = dynamic_cast<SMESH_LengthFromEdges*> (theHyp);
123 ASSERT(_hypLengthFromEdges);
133 if (_maxElementArea > 0)
135 _edgeLength = 2*sqrt(_maxElementArea); // triangles : minorant
138 else isOk = (_hypLengthFromEdges != NULL); // **** check mode
141 //SCRUTE(_edgeLength);
142 //SCRUTE(_maxElementArea);
147 //=============================================================================
151 //=============================================================================
153 bool SMESH_MEFISTO_2D::Compute(SMESH_Mesh& aMesh,
154 const TopoDS_Shape& aShape)
156 MESSAGE("SMESH_MEFISTO_2D::Compute");
158 if (_hypLengthFromEdges)
159 _edgeLength = ComputeEdgeElementLength(aMesh, aShape);
162 const Handle(SMESHDS_Mesh)& meshDS = aMesh.GetMeshDS();
163 SMESH_subMesh* theSubMesh = aMesh.GetSubMesh(aShape);
165 const TopoDS_Face& FF = TopoDS::Face(aShape);
166 bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
167 TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
169 Z nblf; //nombre de lignes fermees (enveloppe en tete)
170 Z *nudslf=NULL; //numero du dernier sommet de chaque ligne fermee
172 Z nbpti=0; //nombre points internes futurs sommets de la triangulation
181 Z nutysu=1; // 1: il existe un fonction areteideale_()
182 // Z nutysu=0; // 0: on utilise aretmx
183 R aretmx=_edgeLength; // longueur max aretes future triangulation
186 nblf = NumberOfWires(F);
189 nudslf = new Z[1+nblf];
194 const TopoDS_Wire OW1 = BRepTools::OuterWire(F);
195 nbpnt += NumberOfPoints (aMesh, OW1);
196 nudslf [iw++] = nbpnt;
199 for (TopExp_Explorer exp (F, TopAbs_WIRE); exp.More(); exp.Next())
201 const TopoDS_Wire& W = TopoDS::Wire(exp.Current());
204 nbpnt += NumberOfPoints (aMesh, W);
205 nudslf [iw++] = nbpnt;
210 uvslf = new R2[nudslf[nblf]];
211 //SCRUTE(nudslf[nblf]);
214 map<int,int> mefistoToDS; // correspondence mefisto index--> points IDNodes
215 TopoDS_Wire OW = BRepTools::OuterWire(F);
216 LoadPoints (aMesh, F, OW, uvslf, m, mefistoToDS);
219 for (TopExp_Explorer exp (F, TopAbs_WIRE); exp.More(); exp.Next())
221 const TopoDS_Wire& W = TopoDS::Wire(exp.Current());
224 LoadPoints (aMesh, F, W, uvslf, m, mefistoToDS);
228 // SCRUTE(nudslf[nblf]);
229 // for (int i=0; i<=nblf; i++)
231 // MESSAGE(" -+- " <<i<< " "<< nudslf[i]);
233 // for (int i=0; i<nudslf[nblf]; i++)
235 // MESSAGE(" -+- " <<i<< " "<< uvslf[i]);
241 MESSAGE("MEFISTO triangulation ...");
244 aptrte( nutysu, aretmx,
247 nbst, uvst, nbt, nust,
252 MESSAGE("... End Triangulation");
255 StoreResult (aMesh, nbst, uvst, nbt, nust, F,
256 faceIsForward, mefistoToDS);
261 MESSAGE("Error in Triangulation");
264 if (nudslf != NULL) delete [] nudslf;
265 if (uvslf != NULL) delete [] uvslf;
266 if (uvst != NULL) delete [] uvst;
267 if (nust != NULL) delete [] nust;
271 //=============================================================================
275 //=============================================================================
277 void SMESH_MEFISTO_2D::LoadPoints(SMESH_Mesh& aMesh,
278 const TopoDS_Face& FF,
279 const TopoDS_Wire& WW,
282 map<int,int>& mefistoToDS)
284 MESSAGE("SMESH_MEFISTO_2D::LoadPoints");
286 Handle (SMDS_Mesh) meshDS = aMesh.GetMeshDS();
290 TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
291 ComputeScaleOnFace(aMesh, F, scalex, scaley);
293 TopoDS_Wire W = TopoDS::Wire(WW.Oriented(TopAbs_FORWARD));
294 BRepTools_WireExplorer wexp(W,F);
295 for (wexp.Init(W,F);wexp.More(); wexp.Next())
297 const TopoDS_Edge& E = wexp.Current();
299 // --- IDNodes of first and last Vertex
301 TopoDS_Vertex VFirst, VLast;
302 TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
304 ASSERT(!VFirst.IsNull());
305 SMESH_subMesh* firstSubMesh = aMesh.GetSubMesh(VFirst);
306 const TColStd_ListOfInteger& lidf
307 = firstSubMesh->GetSubMeshDS()->GetIDNodes();
308 int idFirst= lidf.First();
311 ASSERT(!VLast.IsNull());
312 SMESH_subMesh* lastSubMesh = aMesh.GetSubMesh(VLast);
313 const TColStd_ListOfInteger& lidl
314 = lastSubMesh->GetSubMeshDS()->GetIDNodes();
315 int idLast= lidl.First();
318 // --- edge internal IDNodes (relies on good order storage, not checked)
320 int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
324 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E,F,f,l);
326 const TColStd_ListOfInteger& indElt
327 = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetIDNodes();
328 TColStd_ListIteratorOfListOfInteger ite(indElt);
330 //SCRUTE(indElt.Extent());
331 ASSERT(nbPoints == indElt.Extent());
332 bool isForward = (E.Orientation() == TopAbs_FORWARD);
333 map<double,int> params;
334 for (; ite.More(); ite.Next())
336 int nodeId = ite.Value();
337 Handle (SMDS_MeshElement) elt = meshDS->FindNode(nodeId);
338 Handle (SMDS_MeshNode) node = meshDS->GetNode(1, elt);
339 Handle (SMDS_EdgePosition) epos
340 = Handle (SMDS_EdgePosition)::DownCast(node->GetPosition());
341 double param = epos->GetUParameter();
342 params[param] = nodeId;
343 // MESSAGE(" " << param << " " << params[param]);
346 // --- load 2D values into MEFISTO structure,
347 // add IDNodes in mefistoToDS map
349 if (E.Orientation() == TopAbs_FORWARD)
351 gp_Pnt2d p = C2d->Value(f); // first point = Vertex Forward
352 uvslf [m].x = scalex * p.X();
353 uvslf [m].y = scaley * p.Y();
354 mefistoToDS[m+1] = idFirst;
355 //MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
356 //MESSAGE("__ f "<<f<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
358 map<double,int>::iterator itp = params.begin();
359 for (Standard_Integer i = 1; i<=nbPoints; i++) // nbPoints internal
361 double param = (*itp).first;
362 gp_Pnt2d p = C2d->Value(param);
363 uvslf [m].x = scalex * p.X();
364 uvslf [m].y = scaley * p.Y();
365 mefistoToDS[m+1] = (*itp).second;
366 // MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
367 // MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
374 gp_Pnt2d p = C2d->Value(l); // last point = Vertex Reversed
375 uvslf [m].x = scalex * p.X();
376 uvslf [m].y = scaley * p.Y();
377 mefistoToDS[m+1] = idLast;
378 // MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
379 // MESSAGE("__ l "<<l<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
381 map<double,int>::reverse_iterator itp = params.rbegin();
382 for (Standard_Integer i = nbPoints ; i >= 1; i--)
384 double param = (*itp).first;
385 gp_Pnt2d p = C2d->Value(param);
386 uvslf [m].x = scalex * p.X();
387 uvslf [m].y = scaley * p.Y();
388 mefistoToDS[m+1] = (*itp).second;
389 // MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
390 // MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
398 //=============================================================================
402 //=============================================================================
404 // **** a mettre dans SMESH_Algo ou SMESH_2D_Algo
406 void SMESH_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh& aMesh,
407 const TopoDS_Face& aFace,
411 //MESSAGE("SMESH_MEFISTO_2D::ComputeScaleOnFace");
412 TopoDS_Face F = TopoDS::Face(aFace.Oriented(TopAbs_FORWARD));
413 TopoDS_Wire W = BRepTools::OuterWire(F);
415 BRepTools_WireExplorer wexp(W,F);
417 double xmin = 1.e300; // min & max of face 2D parametric coord.
418 double xmax = -1.e300;
419 double ymin = 1.e300;
420 double ymax = -1.e300;
424 for (wexp.Init(W,F);wexp.More(); wexp.Next())
426 const TopoDS_Edge& E = wexp.Current();
428 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E,F,f,l);
429 for (int i = 0; i<= nbp; i++)
431 double param = f + (double(i)/double(nbp))*(l-f);
432 gp_Pnt2d p = C2d->Value(param);
433 if (p.X() < xmin) xmin = p.X();
434 if (p.X() > xmax) xmax = p.X();
435 if (p.Y() < ymin) ymin = p.Y();
436 if (p.Y() > ymax) ymax = p.Y();
437 // MESSAGE(" "<< f<<" "<<l<<" "<<param<<" "<<xmin<<" "<<xmax<<" "<<ymin<<" "<<ymax);
444 double xmoy = (xmax + xmin)/2.;
445 double ymoy = (ymax + ymin)/2.;
447 Handle(Geom_Surface) S = BRep_Tool::Surface(F); // 3D surface
451 gp_Pnt PX0 = S->Value(xmin, ymoy);
452 gp_Pnt PY0 = S->Value(xmoy, ymin);
453 for (Standard_Integer i = 1; i<= nbp; i++)
455 double x = xmin + (double(i)/double(nbp))*(xmax-xmin);
456 gp_Pnt PX = S->Value(x,ymoy);
457 double y = ymin + (double(i)/double(nbp))*(ymax-ymin);
458 gp_Pnt PY = S->Value(xmoy,y);
459 length_x += PX.Distance(PX0);
460 length_y += PY.Distance(PY0);
461 PX0.SetCoord(PX.X(),PX.Y(),PX.Z());
462 PY0.SetCoord(PY.X(),PY.Y(),PY.Z());
466 scalex = length_x/(xmax - xmin);
467 scaley = length_y/(ymax - ymin);
474 //=============================================================================
478 //=============================================================================
480 void SMESH_MEFISTO_2D::StoreResult (SMESH_Mesh& aMesh,
481 Z nbst, R2* uvst, Z nbt, Z* nust,
482 const TopoDS_Face& F, bool faceIsForward,
483 map<int,int>& mefistoToDS)
487 ComputeScaleOnFace(aMesh, F, scalex, scaley);
489 Handle (SMESHDS_Mesh) meshDS = aMesh.GetMeshDS();
492 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
494 for ( n=0; n<nbst; n++ )
496 double u = uvst[n][0]/scalex;
497 double v = uvst[n][1]/scaley;
498 gp_Pnt P = S->Value(u,v);
500 if (mefistoToDS.find(n+1) == mefistoToDS.end())
502 int nodeId = meshDS->AddNode(P.X(), P.Y(), P.Z());
503 Handle (SMDS_MeshElement) elt = meshDS->FindNode(nodeId);
504 Handle (SMDS_MeshNode) node = meshDS->GetNode(1, elt);
505 meshDS->SetNodeOnFace(node, F);
507 //MESSAGE(nodeId<<" "<<P.X()<<" "<<P.Y()<<" "<<P.Z());
508 mefistoToDS[n+1] = nodeId;
509 //MESSAGE(" "<<n<<" "<<mefistoToDS[n+1]);
510 Handle (SMDS_FacePosition) fpos
511 = Handle (SMDS_FacePosition)::DownCast(node->GetPosition());
512 fpos->SetUParameter(u);
513 fpos->SetVParameter(v);
520 //SCRUTE(faceIsForward);
521 for ( n=1; n<=nbt; n++ )
523 int inode1 = nust[m++];
524 int inode2 = nust[m++];
525 int inode3 = nust[m++];
527 int nodeId1 = mefistoToDS[inode1];
528 int nodeId2 = mefistoToDS[inode2];
529 int nodeId3 = mefistoToDS[inode3];
530 //MESSAGE("-- "<<inode1<<" "<<inode2<<" "<<inode3<<" ++ "<<nodeId1<<" "<<nodeId2<<" "<<nodeId3);
532 // triangle points must be in trigonometric order if face is Forward
533 // else they must be put clockwise
535 bool triangleIsWellOriented = faceIsForward;
537 if (triangleIsWellOriented)
539 faceId = meshDS->AddFace(nodeId1, nodeId2, nodeId3);
543 faceId = meshDS->AddFace(nodeId1, nodeId3, nodeId2);
545 Handle (SMDS_MeshElement) elt = meshDS->FindElement(faceId);
546 meshDS->SetMeshElementOnShape(elt, F);
551 //=============================================================================
555 //=============================================================================
557 double SMESH_MEFISTO_2D::ComputeEdgeElementLength(SMESH_Mesh& aMesh,
558 const TopoDS_Shape& aShape)
560 MESSAGE("SMESH_MEFISTO_2D::ComputeEdgeElementLength");
561 // **** a mettre dans SMESH_2D_Algo ?
563 const TopoDS_Face& FF = TopoDS::Face(aShape);
564 bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
565 TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
567 double meanElementLength = 100;
568 double wireLength =0;
569 int wireElementsNumber =0;
570 for (TopExp_Explorer exp (F, TopAbs_WIRE); exp.More(); exp.Next())
572 const TopoDS_Wire& W = TopoDS::Wire(exp.Current());
573 for (TopExp_Explorer expe(W,TopAbs_EDGE); expe.More(); expe.Next())
575 const TopoDS_Edge& E = TopoDS::Edge(expe.Current());
576 int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
577 double length = EdgeLength(E);
578 wireLength += length;
579 wireElementsNumber += nb;
582 if (wireElementsNumber)
583 meanElementLength = wireLength/wireElementsNumber;
584 //SCRUTE(meanElementLength);
585 return meanElementLength;
588 //=============================================================================
592 //=============================================================================
594 ostream & SMESH_MEFISTO_2D::SaveTo(ostream & save)
599 //=============================================================================
603 //=============================================================================
605 istream & SMESH_MEFISTO_2D::LoadFrom(istream & load)
607 return load >> (*this);
610 //=============================================================================
614 //=============================================================================
616 ostream & operator << (ostream & save, SMESH_MEFISTO_2D & hyp)
621 //=============================================================================
625 //=============================================================================
627 istream & operator >> (istream & load, SMESH_MEFISTO_2D & hyp)