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"
35 #include "StdMeshers_MaxElementArea.hxx"
36 #include "StdMeshers_LengthFromEdges.hxx"
41 #include "SMDS_MeshElement.hxx"
42 #include "SMDS_MeshNode.hxx"
43 #include "SMDS_EdgePosition.hxx"
44 #include "SMDS_FacePosition.hxx"
46 #include "utilities.h"
48 #include <TopoDS_Face.hxx>
49 #include <TopoDS_Edge.hxx>
50 #include <TopoDS_Shape.hxx>
51 #include <Geom_Surface.hxx>
52 #include <GeomAdaptor_Curve.hxx>
53 #include <Geom2d_Curve.hxx>
54 #include <gp_Pnt2d.hxx>
55 #include <BRep_Tool.hxx>
56 #include <BRepTools.hxx>
57 #include <BRepTools_WireExplorer.hxx>
58 #include <GCPnts_AbscissaPoint.hxx>
59 #include <GCPnts_UniformAbscissa.hxx>
60 #include <TColStd_ListIteratorOfListOfInteger.hxx>
65 //=============================================================================
69 //=============================================================================
71 StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D(int hypId, int studyId,
72 SMESH_Gen * gen):SMESH_2D_Algo(hypId, studyId, gen)
74 MESSAGE("StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D");
76 // _shapeType = TopAbs_FACE;
77 _shapeType = (1 << TopAbs_FACE);
78 _compatibleHypothesis.push_back("MaxElementArea");
79 _compatibleHypothesis.push_back("LengthFromEdges");
83 _hypMaxElementArea = NULL;
84 _hypLengthFromEdges = NULL;
87 //=============================================================================
91 //=============================================================================
93 StdMeshers_MEFISTO_2D::~StdMeshers_MEFISTO_2D()
95 MESSAGE("StdMeshers_MEFISTO_2D::~StdMeshers_MEFISTO_2D");
98 //=============================================================================
102 //=============================================================================
104 bool StdMeshers_MEFISTO_2D::CheckHypothesis
106 const TopoDS_Shape& aShape,
107 SMESH_Hypothesis::Hypothesis_Status& aStatus)
109 //MESSAGE("StdMeshers_MEFISTO_2D::CheckHypothesis");
111 _hypMaxElementArea = NULL;
112 _hypLengthFromEdges = NULL;
114 list <const SMESHDS_Hypothesis * >::const_iterator itl;
115 const SMESHDS_Hypothesis *theHyp;
117 const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
118 int nbHyp = hyps.size();
121 aStatus = SMESH_Hypothesis::HYP_MISSING;
122 return false; // can't work with no hypothesis
126 theHyp = (*itl); // use only the first hypothesis
128 string hypName = theHyp->GetName();
129 int hypId = theHyp->GetID();
134 if (hypName == "MaxElementArea")
136 _hypMaxElementArea = static_cast<const StdMeshers_MaxElementArea *>(theHyp);
137 ASSERT(_hypMaxElementArea);
138 _maxElementArea = _hypMaxElementArea->GetMaxArea();
141 aStatus = SMESH_Hypothesis::HYP_OK;
144 else if (hypName == "LengthFromEdges")
146 _hypLengthFromEdges = static_cast<const StdMeshers_LengthFromEdges *>(theHyp);
147 ASSERT(_hypLengthFromEdges);
151 aStatus = SMESH_Hypothesis::HYP_OK;
154 aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
159 if (_maxElementArea > 0)
161 _edgeLength = 2 * sqrt(_maxElementArea); // triangles : minorant
165 isOk = (_hypLengthFromEdges != NULL); // **** check mode
167 aStatus = SMESH_Hypothesis::HYP_BAD_PARAMETER;
170 //SCRUTE(_edgeLength);
171 //SCRUTE(_maxElementArea);
175 //=============================================================================
179 //=============================================================================
181 bool StdMeshers_MEFISTO_2D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape)
183 MESSAGE("StdMeshers_MEFISTO_2D::Compute");
185 if (_hypLengthFromEdges)
186 _edgeLength = ComputeEdgeElementLength(aMesh, aShape);
189 const SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
190 SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape);
192 const TopoDS_Face & FF = TopoDS::Face(aShape);
193 bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
194 TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
196 Z nblf; //nombre de lignes fermees (enveloppe en tete)
197 Z *nudslf = NULL; //numero du dernier sommet de chaque ligne fermee
199 Z nbpti = 0; //nombre points internes futurs sommets de la triangulation
208 Z nutysu = 1; // 1: il existe un fonction areteideale_()
209 // Z nutysu=0; // 0: on utilise aretmx
210 R aretmx = _edgeLength; // longueur max aretes future triangulation
213 nblf = NumberOfWires(F);
216 nudslf = new Z[1 + nblf];
221 const TopoDS_Wire OW1 = BRepTools::OuterWire(F);
222 nbpnt += NumberOfPoints(aMesh, OW1);
223 nudslf[iw++] = nbpnt;
226 for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next())
228 const TopoDS_Wire & W = TopoDS::Wire(exp.Current());
231 nbpnt += NumberOfPoints(aMesh, W);
232 nudslf[iw++] = nbpnt;
237 uvslf = new R2[nudslf[nblf]];
238 //SCRUTE(nudslf[nblf]);
241 map<int, const SMDS_MeshNode*> mefistoToDS; // correspondence mefisto index--> points IDNodes
242 TopoDS_Wire OW = BRepTools::OuterWire(F);
243 LoadPoints(aMesh, F, OW, uvslf, m, mefistoToDS);
246 for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next())
248 const TopoDS_Wire & W = TopoDS::Wire(exp.Current());
251 LoadPoints(aMesh, F, W, uvslf, m, mefistoToDS);
255 // SCRUTE(nudslf[nblf]);
256 // for (int i=0; i<=nblf; i++)
258 // MESSAGE(" -+- " <<i<< " "<< nudslf[i]);
260 // for (int i=0; i<nudslf[nblf]; i++)
262 // MESSAGE(" -+- " <<i<< " "<< uvslf[i]);
268 MESSAGE("MEFISTO triangulation ...");
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);
280 StoreResult(aMesh, nbst, uvst, nbt, nust, F,
281 faceIsForward, mefistoToDS);
286 MESSAGE("Error in Triangulation");
300 //=============================================================================
304 //=============================================================================
306 void StdMeshers_MEFISTO_2D::LoadPoints(SMESH_Mesh & aMesh,
307 const TopoDS_Face & FF,
308 const TopoDS_Wire & WW, R2 * uvslf, int &m,
309 map<int, const SMDS_MeshNode*>&mefistoToDS)
311 MESSAGE("StdMeshers_MEFISTO_2D::LoadPoints");
313 SMDS_Mesh * meshDS = aMesh.GetMeshDS();
317 TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
318 ComputeScaleOnFace(aMesh, F, scalex, scaley);
320 TopoDS_Wire W = TopoDS::Wire(WW.Oriented(TopAbs_FORWARD));
321 BRepTools_WireExplorer wexp(W, F);
322 for (wexp.Init(W, F); wexp.More(); wexp.Next())
324 const TopoDS_Edge & E = wexp.Current();
326 // --- IDNodes of first and last Vertex
328 TopoDS_Vertex VFirst, VLast;
329 TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
331 ASSERT(!VFirst.IsNull());
332 SMDS_NodeIteratorPtr lid=
333 aMesh.GetSubMesh(VFirst)->GetSubMeshDS()->GetNodes();
334 const SMDS_MeshNode* idFirst = lid->next();
336 ASSERT(!VLast.IsNull());
337 lid=aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes();
338 const SMDS_MeshNode* idLast = lid->next();
340 // --- edge internal IDNodes (relies on good order storage, not checked)
342 int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
346 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
348 SMDS_NodeIteratorPtr ite= aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes();
350 bool isForward = (E.Orientation() == TopAbs_FORWARD);
351 map<double, const SMDS_MeshNode*> params;
355 const SMDS_MeshNode * node = ite->next();
356 const SMDS_EdgePosition* epos =
357 static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
358 double param = epos->GetUParameter();
359 params[param] = node;
361 // --- load 2D values into MEFISTO structure,
362 // add IDNodes in mefistoToDS map
364 if (E.Orientation() == TopAbs_FORWARD)
366 gp_Pnt2d p = C2d->Value(f); // first point = Vertex Forward
367 uvslf[m].x = scalex * p.X();
368 uvslf[m].y = scaley * p.Y();
369 mefistoToDS[m + 1] = idFirst;
370 //MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
371 //MESSAGE("__ f "<<f<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
373 map<double, const SMDS_MeshNode*>::iterator itp = params.begin();
374 for (int i = 1; i <= nbPoints; i++) // nbPoints internal
376 double param = (*itp).first;
377 gp_Pnt2d p = C2d->Value(param);
378 uvslf[m].x = scalex * p.X();
379 uvslf[m].y = scaley * p.Y();
380 mefistoToDS[m + 1] = (*itp).second;
381 // MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
382 // MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
389 gp_Pnt2d p = C2d->Value(l); // last point = Vertex Reversed
390 uvslf[m].x = scalex * p.X();
391 uvslf[m].y = scaley * p.Y();
392 mefistoToDS[m + 1] = idLast;
393 // MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
394 // MESSAGE("__ l "<<l<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
396 map<double, const SMDS_MeshNode*>::reverse_iterator itp = params.rbegin();
397 for (int i = nbPoints; i >= 1; i--)
399 double param = (*itp).first;
400 gp_Pnt2d p = C2d->Value(param);
401 uvslf[m].x = scalex * p.X();
402 uvslf[m].y = scaley * p.Y();
403 mefistoToDS[m + 1] = (*itp).second;
404 // MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
405 // MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
413 //=============================================================================
417 //=============================================================================
419 // **** a mettre dans SMESH_Algo ou SMESH_2D_Algo
421 void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh & aMesh,
422 const TopoDS_Face & aFace, double &scalex, double &scaley)
424 //MESSAGE("StdMeshers_MEFISTO_2D::ComputeScaleOnFace");
425 TopoDS_Face F = TopoDS::Face(aFace.Oriented(TopAbs_FORWARD));
426 TopoDS_Wire W = BRepTools::OuterWire(F);
428 BRepTools_WireExplorer wexp(W, F);
430 double xmin = 1.e300; // min & max of face 2D parametric coord.
431 double xmax = -1.e300;
432 double ymin = 1.e300;
433 double ymax = -1.e300;
437 for (wexp.Init(W, F); wexp.More(); wexp.Next())
439 const TopoDS_Edge & E = wexp.Current();
441 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
442 for (int i = 0; i <= nbp; i++)
444 double param = f + (double (i) / double (nbp))*(l - f);
445 gp_Pnt2d p = C2d->Value(param);
454 // MESSAGE(" "<< f<<" "<<l<<" "<<param<<" "<<xmin<<" "<<xmax<<" "<<ymin<<" "<<ymax);
461 double xmoy = (xmax + xmin) / 2.;
462 double ymoy = (ymax + ymin) / 2.;
464 Handle(Geom_Surface) S = BRep_Tool::Surface(F); // 3D surface
468 gp_Pnt PX0 = S->Value(xmin, ymoy);
469 gp_Pnt PY0 = S->Value(xmoy, ymin);
470 for (int i = 1; i <= nbp; i++)
472 double x = xmin + (double (i) / double (nbp))*(xmax - xmin);
473 gp_Pnt PX = S->Value(x, ymoy);
474 double y = ymin + (double (i) / double (nbp))*(ymax - ymin);
475 gp_Pnt PY = S->Value(xmoy, y);
476 length_x += PX.Distance(PX0);
477 length_y += PY.Distance(PY0);
478 PX0.SetCoord(PX.X(), PX.Y(), PX.Z());
479 PY0.SetCoord(PY.X(), PY.Y(), PY.Z());
483 scalex = length_x / (xmax - xmin);
484 scaley = length_y / (ymax - ymin);
491 //=============================================================================
495 //=============================================================================
497 void StdMeshers_MEFISTO_2D::StoreResult(SMESH_Mesh & aMesh,
498 Z nbst, R2 * uvst, Z nbt, Z * nust,
499 const TopoDS_Face & F, bool faceIsForward,
500 map<int, const SMDS_MeshNode*>&mefistoToDS)
504 ComputeScaleOnFace(aMesh, F, scalex, scaley);
506 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
509 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
511 for (n = 0; n < nbst; n++)
513 double u = uvst[n][0] / scalex;
514 double v = uvst[n][1] / scaley;
515 gp_Pnt P = S->Value(u, v);
517 if (mefistoToDS.find(n + 1) == mefistoToDS.end())
519 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
520 meshDS->SetNodeOnFace(node, F);
522 //MESSAGE(nodeId<<" "<<P.X()<<" "<<P.Y()<<" "<<P.Z());
523 mefistoToDS[n + 1] = node;
524 //MESSAGE(" "<<n<<" "<<mefistoToDS[n+1]);
525 SMDS_FacePosition* fpos =
526 static_cast<SMDS_FacePosition*>(node->GetPosition().get());
527 fpos->SetUParameter(u);
528 fpos->SetVParameter(v);
535 //SCRUTE(faceIsForward);
536 for (n = 1; n <= nbt; n++)
538 int inode1 = nust[m++];
539 int inode2 = nust[m++];
540 int inode3 = nust[m++];
542 const SMDS_MeshNode *n1, *n2, *n3;
543 n1 = mefistoToDS[inode1];
544 n2 = mefistoToDS[inode2];
545 n3 = mefistoToDS[inode3];
546 //MESSAGE("-- "<<inode1<<" "<<inode2<<" "<<inode3<<" ++ "<<nodeId1<<" "<<nodeId2<<" "<<nodeId3);
548 // triangle points must be in trigonometric order if face is Forward
549 // else they must be put clockwise
551 bool triangleIsWellOriented = faceIsForward;
553 SMDS_MeshElement * elt;
554 if (triangleIsWellOriented)
555 elt = meshDS->AddFace(n1, n2, n3);
557 elt = meshDS->AddFace(n1, n3, n2);
559 meshDS->SetMeshElementOnShape(elt, F);
564 //=============================================================================
568 //=============================================================================
570 double StdMeshers_MEFISTO_2D::ComputeEdgeElementLength(SMESH_Mesh & aMesh,
571 const TopoDS_Shape & aShape)
573 MESSAGE("StdMeshers_MEFISTO_2D::ComputeEdgeElementLength");
574 // **** a mettre dans SMESH_2D_Algo ?
576 const TopoDS_Face & FF = TopoDS::Face(aShape);
577 bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
578 TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
580 double meanElementLength = 100;
581 double wireLength = 0;
582 int wireElementsNumber = 0;
583 for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next())
585 const TopoDS_Wire & W = TopoDS::Wire(exp.Current());
586 for (TopExp_Explorer expe(W, TopAbs_EDGE); expe.More(); expe.Next())
588 const TopoDS_Edge & E = TopoDS::Edge(expe.Current());
589 int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
590 double length = EdgeLength(E);
591 wireLength += length;
592 wireElementsNumber += nb;
595 if (wireElementsNumber)
596 meanElementLength = wireLength / wireElementsNumber;
597 //SCRUTE(meanElementLength);
598 return meanElementLength;
601 //=============================================================================
605 //=============================================================================
607 ostream & StdMeshers_MEFISTO_2D::SaveTo(ostream & save)
612 //=============================================================================
616 //=============================================================================
618 istream & StdMeshers_MEFISTO_2D::LoadFrom(istream & load)
623 //=============================================================================
627 //=============================================================================
629 ostream & operator <<(ostream & save, StdMeshers_MEFISTO_2D & hyp)
631 return hyp.SaveTo( save );
634 //=============================================================================
638 //=============================================================================
640 istream & operator >>(istream & load, StdMeshers_MEFISTO_2D & hyp)
642 return hyp.LoadFrom( load );