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 : SMESH_MEFISTO_2D.cxx
25 // Author : Paul RASCLE, EDF
31 #include "SMESH_MEFISTO_2D.hxx"
32 #include "SMESH_Gen.hxx"
33 #include "SMESH_Mesh.hxx"
35 #include "SMESH_MaxElementArea.hxx"
36 #include "SMESH_LengthFromEdges.hxx"
41 #include "SMESHDS_ListOfPtrHypothesis.hxx"
42 #include "SMESHDS_ListIteratorOfListOfPtrHypothesis.hxx"
43 #include "SMDS_MeshElement.hxx"
44 #include "SMDS_MeshNode.hxx"
45 #include "SMDS_EdgePosition.hxx"
46 #include "SMDS_FacePosition.hxx"
48 #include "utilities.h"
50 #include <TopoDS_Face.hxx>
51 #include <TopoDS_Edge.hxx>
52 #include <TopoDS_Shape.hxx>
53 #include <Geom_Surface.hxx>
54 #include <GeomAdaptor_Curve.hxx>
55 #include <Geom2d_Curve.hxx>
56 #include <gp_Pnt2d.hxx>
57 #include <BRep_Tool.hxx>
58 #include <BRepTools.hxx>
59 #include <BRepTools_WireExplorer.hxx>
60 #include <GCPnts_AbscissaPoint.hxx>
61 #include <GCPnts_UniformAbscissa.hxx>
62 #include <TColStd_ListIteratorOfListOfInteger.hxx>
67 //=============================================================================
71 //=============================================================================
73 SMESH_MEFISTO_2D::SMESH_MEFISTO_2D(int hypId, int studyId, SMESH_Gen* gen)
74 : SMESH_2D_Algo(hypId, studyId, gen)
76 MESSAGE("SMESH_MEFISTO_2D::SMESH_MEFISTO_2D");
78 // _shapeType = TopAbs_FACE;
79 _shapeType = (1<<TopAbs_FACE);
80 _compatibleHypothesis.push_back("MaxElementArea");
81 _compatibleHypothesis.push_back("LengthFromEdges");
85 _hypMaxElementArea = NULL;
86 _hypLengthFromEdges = NULL;
89 //=============================================================================
93 //=============================================================================
95 SMESH_MEFISTO_2D::~SMESH_MEFISTO_2D()
97 MESSAGE("SMESH_MEFISTO_2D::~SMESH_MEFISTO_2D");
100 //=============================================================================
104 //=============================================================================
106 bool SMESH_MEFISTO_2D::CheckHypothesis(SMESH_Mesh& aMesh,
107 const TopoDS_Shape& aShape)
109 //MESSAGE("SMESH_MEFISTO_2D::CheckHypothesis");
111 _hypMaxElementArea = NULL;
112 _hypLengthFromEdges = NULL;
114 list<SMESHDS_Hypothesis*>::const_iterator itl;
115 SMESHDS_Hypothesis* theHyp;
117 const list<SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
118 int nbHyp = hyps.size();
119 if (nbHyp != 1) return false; // only one compatible hypothesis allowed
124 string hypName = theHyp->GetName();
125 int hypId = theHyp->GetID();
130 if (hypName == "MaxElementArea")
132 _hypMaxElementArea = dynamic_cast<SMESH_MaxElementArea*> (theHyp);
133 ASSERT(_hypMaxElementArea);
134 _maxElementArea = _hypMaxElementArea->GetMaxArea();
139 if (hypName == "LengthFromEdges")
141 _hypLengthFromEdges = dynamic_cast<SMESH_LengthFromEdges*> (theHyp);
142 ASSERT(_hypLengthFromEdges);
152 if (_maxElementArea > 0)
154 _edgeLength = 2*sqrt(_maxElementArea); // triangles : minorant
157 else isOk = (_hypLengthFromEdges != NULL); // **** check mode
160 //SCRUTE(_edgeLength);
161 //SCRUTE(_maxElementArea);
166 //=============================================================================
170 //=============================================================================
172 bool SMESH_MEFISTO_2D::Compute(SMESH_Mesh& aMesh,
173 const TopoDS_Shape& aShape)
175 MESSAGE("SMESH_MEFISTO_2D::Compute");
177 if (_hypLengthFromEdges)
178 _edgeLength = ComputeEdgeElementLength(aMesh, aShape);
181 const Handle(SMESHDS_Mesh)& meshDS = aMesh.GetMeshDS();
182 SMESH_subMesh* theSubMesh = aMesh.GetSubMesh(aShape);
184 const TopoDS_Face& FF = TopoDS::Face(aShape);
185 bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
186 TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
188 Z nblf; //nombre de lignes fermees (enveloppe en tete)
189 Z *nudslf=NULL; //numero du dernier sommet de chaque ligne fermee
191 Z nbpti=0; //nombre points internes futurs sommets de la triangulation
200 Z nutysu=1; // 1: il existe un fonction areteideale_()
201 // Z nutysu=0; // 0: on utilise aretmx
202 R aretmx=_edgeLength; // longueur max aretes future triangulation
205 nblf = NumberOfWires(F);
208 nudslf = new Z[1+nblf];
213 const TopoDS_Wire OW1 = BRepTools::OuterWire(F);
214 nbpnt += NumberOfPoints (aMesh, OW1);
215 nudslf [iw++] = nbpnt;
218 for (TopExp_Explorer exp (F, TopAbs_WIRE); exp.More(); exp.Next())
220 const TopoDS_Wire& W = TopoDS::Wire(exp.Current());
223 nbpnt += NumberOfPoints (aMesh, W);
224 nudslf [iw++] = nbpnt;
229 uvslf = new R2[nudslf[nblf]];
230 //SCRUTE(nudslf[nblf]);
233 map<int,int> mefistoToDS; // correspondence mefisto index--> points IDNodes
234 TopoDS_Wire OW = BRepTools::OuterWire(F);
235 LoadPoints (aMesh, F, OW, uvslf, m, mefistoToDS);
238 for (TopExp_Explorer exp (F, TopAbs_WIRE); exp.More(); exp.Next())
240 const TopoDS_Wire& W = TopoDS::Wire(exp.Current());
243 LoadPoints (aMesh, F, W, uvslf, m, mefistoToDS);
247 // SCRUTE(nudslf[nblf]);
248 // for (int i=0; i<=nblf; i++)
250 // MESSAGE(" -+- " <<i<< " "<< nudslf[i]);
252 // for (int i=0; i<nudslf[nblf]; i++)
254 // MESSAGE(" -+- " <<i<< " "<< uvslf[i]);
260 MESSAGE("MEFISTO triangulation ...");
263 aptrte( nutysu, aretmx,
266 nbst, uvst, nbt, nust,
271 MESSAGE("... End Triangulation");
274 StoreResult (aMesh, nbst, uvst, nbt, nust, F,
275 faceIsForward, mefistoToDS);
280 MESSAGE("Error in Triangulation");
283 if (nudslf != NULL) delete [] nudslf;
284 if (uvslf != NULL) delete [] uvslf;
285 if (uvst != NULL) delete [] uvst;
286 if (nust != NULL) delete [] nust;
290 //=============================================================================
294 //=============================================================================
296 void SMESH_MEFISTO_2D::LoadPoints(SMESH_Mesh& aMesh,
297 const TopoDS_Face& FF,
298 const TopoDS_Wire& WW,
301 map<int,int>& mefistoToDS)
303 MESSAGE("SMESH_MEFISTO_2D::LoadPoints");
305 Handle (SMDS_Mesh) meshDS = aMesh.GetMeshDS();
309 TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
310 ComputeScaleOnFace(aMesh, F, scalex, scaley);
312 TopoDS_Wire W = TopoDS::Wire(WW.Oriented(TopAbs_FORWARD));
313 BRepTools_WireExplorer wexp(W,F);
314 for (wexp.Init(W,F);wexp.More(); wexp.Next())
316 const TopoDS_Edge& E = wexp.Current();
318 // --- IDNodes of first and last Vertex
320 TopoDS_Vertex VFirst, VLast;
321 TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
323 ASSERT(!VFirst.IsNull());
324 SMESH_subMesh* firstSubMesh = aMesh.GetSubMesh(VFirst);
325 const TColStd_ListOfInteger& lidf
326 = firstSubMesh->GetSubMeshDS()->GetIDNodes();
327 int idFirst= lidf.First();
330 ASSERT(!VLast.IsNull());
331 SMESH_subMesh* lastSubMesh = aMesh.GetSubMesh(VLast);
332 const TColStd_ListOfInteger& lidl
333 = lastSubMesh->GetSubMeshDS()->GetIDNodes();
334 int idLast= lidl.First();
337 // --- edge internal IDNodes (relies on good order storage, not checked)
339 int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
343 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E,F,f,l);
345 const TColStd_ListOfInteger& indElt
346 = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetIDNodes();
347 TColStd_ListIteratorOfListOfInteger ite(indElt);
349 //SCRUTE(indElt.Extent());
350 ASSERT(nbPoints == indElt.Extent());
351 bool isForward = (E.Orientation() == TopAbs_FORWARD);
352 map<double,int> params;
353 for (; ite.More(); ite.Next())
355 int nodeId = ite.Value();
356 Handle (SMDS_MeshElement) elt = meshDS->FindNode(nodeId);
357 Handle (SMDS_MeshNode) node = meshDS->GetNode(1, elt);
358 Handle (SMDS_EdgePosition) epos
359 = Handle (SMDS_EdgePosition)::DownCast(node->GetPosition());
360 double param = epos->GetUParameter();
361 params[param] = nodeId;
362 // MESSAGE(" " << param << " " << params[param]);
365 // --- load 2D values into MEFISTO structure,
366 // add IDNodes in mefistoToDS map
368 if (E.Orientation() == TopAbs_FORWARD)
370 gp_Pnt2d p = C2d->Value(f); // first point = Vertex Forward
371 uvslf [m].x = scalex * p.X();
372 uvslf [m].y = scaley * p.Y();
373 mefistoToDS[m+1] = idFirst;
374 //MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
375 //MESSAGE("__ f "<<f<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
377 map<double,int>::iterator itp = params.begin();
378 for (Standard_Integer i = 1; i<=nbPoints; i++) // nbPoints internal
380 double param = (*itp).first;
381 gp_Pnt2d p = C2d->Value(param);
382 uvslf [m].x = scalex * p.X();
383 uvslf [m].y = scaley * p.Y();
384 mefistoToDS[m+1] = (*itp).second;
385 // MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
386 // MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
393 gp_Pnt2d p = C2d->Value(l); // last point = Vertex Reversed
394 uvslf [m].x = scalex * p.X();
395 uvslf [m].y = scaley * p.Y();
396 mefistoToDS[m+1] = idLast;
397 // MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
398 // MESSAGE("__ l "<<l<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
400 map<double,int>::reverse_iterator itp = params.rbegin();
401 for (Standard_Integer i = nbPoints ; i >= 1; i--)
403 double param = (*itp).first;
404 gp_Pnt2d p = C2d->Value(param);
405 uvslf [m].x = scalex * p.X();
406 uvslf [m].y = scaley * p.Y();
407 mefistoToDS[m+1] = (*itp).second;
408 // MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
409 // MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
417 //=============================================================================
421 //=============================================================================
423 // **** a mettre dans SMESH_Algo ou SMESH_2D_Algo
425 void SMESH_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh& aMesh,
426 const TopoDS_Face& aFace,
430 //MESSAGE("SMESH_MEFISTO_2D::ComputeScaleOnFace");
431 TopoDS_Face F = TopoDS::Face(aFace.Oriented(TopAbs_FORWARD));
432 TopoDS_Wire W = BRepTools::OuterWire(F);
434 BRepTools_WireExplorer wexp(W,F);
436 double xmin = 1.e300; // min & max of face 2D parametric coord.
437 double xmax = -1.e300;
438 double ymin = 1.e300;
439 double ymax = -1.e300;
443 for (wexp.Init(W,F);wexp.More(); wexp.Next())
445 const TopoDS_Edge& E = wexp.Current();
447 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E,F,f,l);
448 for (int i = 0; i<= nbp; i++)
450 double param = f + (double(i)/double(nbp))*(l-f);
451 gp_Pnt2d p = C2d->Value(param);
452 if (p.X() < xmin) xmin = p.X();
453 if (p.X() > xmax) xmax = p.X();
454 if (p.Y() < ymin) ymin = p.Y();
455 if (p.Y() > ymax) ymax = p.Y();
456 // MESSAGE(" "<< f<<" "<<l<<" "<<param<<" "<<xmin<<" "<<xmax<<" "<<ymin<<" "<<ymax);
463 double xmoy = (xmax + xmin)/2.;
464 double ymoy = (ymax + ymin)/2.;
466 Handle(Geom_Surface) S = BRep_Tool::Surface(F); // 3D surface
470 gp_Pnt PX0 = S->Value(xmin, ymoy);
471 gp_Pnt PY0 = S->Value(xmoy, ymin);
472 for (Standard_Integer i = 1; i<= nbp; i++)
474 double x = xmin + (double(i)/double(nbp))*(xmax-xmin);
475 gp_Pnt PX = S->Value(x,ymoy);
476 double y = ymin + (double(i)/double(nbp))*(ymax-ymin);
477 gp_Pnt PY = S->Value(xmoy,y);
478 length_x += PX.Distance(PX0);
479 length_y += PY.Distance(PY0);
480 PX0.SetCoord(PX.X(),PX.Y(),PX.Z());
481 PY0.SetCoord(PY.X(),PY.Y(),PY.Z());
485 scalex = length_x/(xmax - xmin);
486 scaley = length_y/(ymax - ymin);
493 //=============================================================================
497 //=============================================================================
499 void SMESH_MEFISTO_2D::StoreResult (SMESH_Mesh& aMesh,
500 Z nbst, R2* uvst, Z nbt, Z* nust,
501 const TopoDS_Face& F, bool faceIsForward,
502 map<int,int>& mefistoToDS)
506 ComputeScaleOnFace(aMesh, F, scalex, scaley);
508 Handle (SMESHDS_Mesh) meshDS = aMesh.GetMeshDS();
511 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
513 for ( n=0; n<nbst; n++ )
515 double u = uvst[n][0]/scalex;
516 double v = uvst[n][1]/scaley;
517 gp_Pnt P = S->Value(u,v);
519 if (mefistoToDS.find(n+1) == mefistoToDS.end())
521 int nodeId = meshDS->AddNode(P.X(), P.Y(), P.Z());
522 Handle (SMDS_MeshElement) elt = meshDS->FindNode(nodeId);
523 Handle (SMDS_MeshNode) node = meshDS->GetNode(1, elt);
524 meshDS->SetNodeOnFace(node, F);
526 //MESSAGE(nodeId<<" "<<P.X()<<" "<<P.Y()<<" "<<P.Z());
527 mefistoToDS[n+1] = nodeId;
528 //MESSAGE(" "<<n<<" "<<mefistoToDS[n+1]);
529 Handle (SMDS_FacePosition) fpos
530 = Handle (SMDS_FacePosition)::DownCast(node->GetPosition());
531 fpos->SetUParameter(u);
532 fpos->SetVParameter(v);
539 //SCRUTE(faceIsForward);
540 for ( n=1; n<=nbt; n++ )
542 int inode1 = nust[m++];
543 int inode2 = nust[m++];
544 int inode3 = nust[m++];
546 int nodeId1 = mefistoToDS[inode1];
547 int nodeId2 = mefistoToDS[inode2];
548 int nodeId3 = mefistoToDS[inode3];
549 //MESSAGE("-- "<<inode1<<" "<<inode2<<" "<<inode3<<" ++ "<<nodeId1<<" "<<nodeId2<<" "<<nodeId3);
551 // triangle points must be in trigonometric order if face is Forward
552 // else they must be put clockwise
554 bool triangleIsWellOriented = faceIsForward;
556 if (triangleIsWellOriented)
558 faceId = meshDS->AddFace(nodeId1, nodeId2, nodeId3);
562 faceId = meshDS->AddFace(nodeId1, nodeId3, nodeId2);
564 Handle (SMDS_MeshElement) elt = meshDS->FindElement(faceId);
565 meshDS->SetMeshElementOnShape(elt, F);
570 //=============================================================================
574 //=============================================================================
576 double SMESH_MEFISTO_2D::ComputeEdgeElementLength(SMESH_Mesh& aMesh,
577 const TopoDS_Shape& aShape)
579 MESSAGE("SMESH_MEFISTO_2D::ComputeEdgeElementLength");
580 // **** a mettre dans SMESH_2D_Algo ?
582 const TopoDS_Face& FF = TopoDS::Face(aShape);
583 bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
584 TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
586 double meanElementLength = 100;
587 double wireLength =0;
588 int wireElementsNumber =0;
589 for (TopExp_Explorer exp (F, TopAbs_WIRE); exp.More(); exp.Next())
591 const TopoDS_Wire& W = TopoDS::Wire(exp.Current());
592 for (TopExp_Explorer expe(W,TopAbs_EDGE); expe.More(); expe.Next())
594 const TopoDS_Edge& E = TopoDS::Edge(expe.Current());
595 int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
596 double length = EdgeLength(E);
597 wireLength += length;
598 wireElementsNumber += nb;
601 if (wireElementsNumber)
602 meanElementLength = wireLength/wireElementsNumber;
603 //SCRUTE(meanElementLength);
604 return meanElementLength;
607 //=============================================================================
611 //=============================================================================
613 ostream & SMESH_MEFISTO_2D::SaveTo(ostream & save)
618 //=============================================================================
622 //=============================================================================
624 istream & SMESH_MEFISTO_2D::LoadFrom(istream & load)
626 return load >> (*this);
629 //=============================================================================
633 //=============================================================================
635 ostream & operator << (ostream & save, SMESH_MEFISTO_2D & hyp)
640 //=============================================================================
644 //=============================================================================
646 istream & operator >> (istream & load, SMESH_MEFISTO_2D & hyp)