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
30 #include "SMESH_MEFISTO_2D.hxx"
31 #include "SMESH_Gen.hxx"
32 #include "SMESH_Mesh.hxx"
34 #include "SMESH_MaxElementArea.hxx"
35 #include "SMESH_LengthFromEdges.hxx"
40 #include "SMDS_MeshElement.hxx"
41 #include "SMDS_MeshNode.hxx"
42 #include "SMDS_EdgePosition.hxx"
43 #include "SMDS_FacePosition.hxx"
45 #include "utilities.h"
47 #include <TopoDS_Face.hxx>
48 #include <TopoDS_Edge.hxx>
49 #include <TopoDS_Shape.hxx>
50 #include <Geom_Surface.hxx>
51 #include <GeomAdaptor_Curve.hxx>
52 #include <Geom2d_Curve.hxx>
53 #include <gp_Pnt2d.hxx>
54 #include <BRep_Tool.hxx>
55 #include <BRepTools.hxx>
56 #include <BRepTools_WireExplorer.hxx>
57 #include <GCPnts_AbscissaPoint.hxx>
58 #include <GCPnts_UniformAbscissa.hxx>
59 #include <TColStd_ListIteratorOfListOfInteger.hxx>
64 //=============================================================================
68 //=============================================================================
70 SMESH_MEFISTO_2D::SMESH_MEFISTO_2D(int hypId, int studyId,
71 SMESH_Gen * gen):SMESH_2D_Algo(hypId, studyId, gen)
73 MESSAGE("SMESH_MEFISTO_2D::SMESH_MEFISTO_2D");
75 // _shapeType = TopAbs_FACE;
76 _shapeType = (1 << TopAbs_FACE);
77 _compatibleHypothesis.push_back("MaxElementArea");
78 _compatibleHypothesis.push_back("LengthFromEdges");
82 _hypMaxElementArea = NULL;
83 _hypLengthFromEdges = NULL;
86 //=============================================================================
90 //=============================================================================
92 SMESH_MEFISTO_2D::~SMESH_MEFISTO_2D()
94 MESSAGE("SMESH_MEFISTO_2D::~SMESH_MEFISTO_2D");
97 //=============================================================================
101 //=============================================================================
103 bool SMESH_MEFISTO_2D::CheckHypothesis(SMESH_Mesh & aMesh,
104 const TopoDS_Shape & aShape)
106 //MESSAGE("SMESH_MEFISTO_2D::CheckHypothesis");
108 _hypMaxElementArea = NULL;
109 _hypLengthFromEdges = NULL;
111 list <const SMESHDS_Hypothesis * >::const_iterator itl;
112 const SMESHDS_Hypothesis *theHyp;
114 const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
115 int nbHyp = hyps.size();
116 if (nbHyp != 1) return false;// only one compatible hypothesis allowed
121 string hypName = theHyp->GetName();
122 int hypId = theHyp->GetID();
127 if (hypName == "MaxElementArea")
129 _hypMaxElementArea = static_cast<const SMESH_MaxElementArea *>(theHyp);
130 ASSERT(_hypMaxElementArea);
131 _maxElementArea = _hypMaxElementArea->GetMaxArea();
136 if (hypName == "LengthFromEdges")
138 _hypLengthFromEdges = static_cast<const SMESH_LengthFromEdges *>(theHyp);
139 ASSERT(_hypLengthFromEdges);
148 if (_maxElementArea > 0)
150 _edgeLength = 2 * sqrt(_maxElementArea); // triangles : minorant
154 isOk = (_hypLengthFromEdges != NULL); // **** check mode
157 //SCRUTE(_edgeLength);
158 //SCRUTE(_maxElementArea);
162 //=============================================================================
166 //=============================================================================
168 bool SMESH_MEFISTO_2D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape)
170 MESSAGE("SMESH_MEFISTO_2D::Compute");
172 if (_hypLengthFromEdges)
173 _edgeLength = ComputeEdgeElementLength(aMesh, aShape);
176 const SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
177 SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape);
179 const TopoDS_Face & FF = TopoDS::Face(aShape);
180 bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
181 TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
183 Z nblf; //nombre de lignes fermees (enveloppe en tete)
184 Z *nudslf = NULL; //numero du dernier sommet de chaque ligne fermee
186 Z nbpti = 0; //nombre points internes futurs sommets de la triangulation
195 Z nutysu = 1; // 1: il existe un fonction areteideale_()
196 // Z nutysu=0; // 0: on utilise aretmx
197 R aretmx = _edgeLength; // longueur max aretes future triangulation
200 nblf = NumberOfWires(F);
203 nudslf = new Z[1 + nblf];
208 const TopoDS_Wire OW1 = BRepTools::OuterWire(F);
209 nbpnt += NumberOfPoints(aMesh, OW1);
210 nudslf[iw++] = nbpnt;
213 for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next())
215 const TopoDS_Wire & W = TopoDS::Wire(exp.Current());
218 nbpnt += NumberOfPoints(aMesh, W);
219 nudslf[iw++] = nbpnt;
224 uvslf = new R2[nudslf[nblf]];
225 //SCRUTE(nudslf[nblf]);
228 map < int, int >mefistoToDS; // correspondence mefisto index--> points IDNodes
229 TopoDS_Wire OW = BRepTools::OuterWire(F);
230 LoadPoints(aMesh, F, OW, uvslf, m, mefistoToDS);
233 for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next())
235 const TopoDS_Wire & W = TopoDS::Wire(exp.Current());
238 LoadPoints(aMesh, F, W, uvslf, m, mefistoToDS);
242 // SCRUTE(nudslf[nblf]);
243 // for (int i=0; i<=nblf; i++)
245 // MESSAGE(" -+- " <<i<< " "<< nudslf[i]);
247 // for (int i=0; i<nudslf[nblf]; i++)
249 // MESSAGE(" -+- " <<i<< " "<< uvslf[i]);
255 MESSAGE("MEFISTO triangulation ...");
258 aptrte(nutysu, aretmx,
259 nblf, nudslf, uvslf, nbpti, uvpti, nbst, uvst, nbt, nust, ierr);
263 MESSAGE("... End Triangulation");
266 StoreResult(aMesh, nbst, uvst, nbt, nust, F,
267 faceIsForward, mefistoToDS);
272 MESSAGE("Error in Triangulation");
286 //=============================================================================
290 //=============================================================================
292 void SMESH_MEFISTO_2D::LoadPoints(SMESH_Mesh & aMesh,
293 const TopoDS_Face & FF,
294 const TopoDS_Wire & WW, R2 * uvslf, int &m, map < int, int >&mefistoToDS)
296 MESSAGE("SMESH_MEFISTO_2D::LoadPoints");
298 SMDS_Mesh * meshDS = aMesh.GetMeshDS();
302 TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
303 ComputeScaleOnFace(aMesh, F, scalex, scaley);
305 TopoDS_Wire W = TopoDS::Wire(WW.Oriented(TopAbs_FORWARD));
306 BRepTools_WireExplorer wexp(W, F);
307 for (wexp.Init(W, F); wexp.More(); wexp.Next())
309 const TopoDS_Edge & E = wexp.Current();
311 // --- IDNodes of first and last Vertex
313 TopoDS_Vertex VFirst, VLast;
314 TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
316 ASSERT(!VFirst.IsNull());
317 SMESH_subMesh *firstSubMesh = aMesh.GetSubMesh(VFirst);
318 const vector<int> & lidf
319 = firstSubMesh->GetSubMeshDS()->GetIDNodes();
320 int idFirst = lidf[0];
323 ASSERT(!VLast.IsNull());
324 SMESH_subMesh *lastSubMesh = aMesh.GetSubMesh(VLast);
325 const vector<int> & lidl
326 = lastSubMesh->GetSubMeshDS()->GetIDNodes();
327 int idLast = lidl[0];
330 // --- edge internal IDNodes (relies on good order storage, not checked)
332 int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
336 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
338 const vector<int> & indElt
339 = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetIDNodes();
341 ASSERT(nbPoints == indElt.size());
342 bool isForward = (E.Orientation() == TopAbs_FORWARD);
343 map < double, int >params;
344 for (int ite=0; ite<indElt.size(); ite++)
346 int nodeId = indElt[ite];
347 const SMDS_MeshNode * node = meshDS->FindNode(nodeId);
348 const SMDS_EdgePosition* epos
349 = static_cast<const SMDS_EdgePosition*>(node->GetPosition());
350 double param = epos->GetUParameter();
351 params[param] = nodeId;
354 // --- load 2D values into MEFISTO structure,
355 // add IDNodes in mefistoToDS map
357 if (E.Orientation() == TopAbs_FORWARD)
359 gp_Pnt2d p = C2d->Value(f); // first point = Vertex Forward
360 uvslf[m].x = scalex * p.X();
361 uvslf[m].y = scaley * p.Y();
362 mefistoToDS[m + 1] = idFirst;
363 //MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
364 //MESSAGE("__ f "<<f<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
366 map < double, int >::iterator itp = params.begin();
367 for (int i = 1; i <= nbPoints; i++) // nbPoints internal
369 double param = (*itp).first;
370 gp_Pnt2d p = C2d->Value(param);
371 uvslf[m].x = scalex * p.X();
372 uvslf[m].y = scaley * p.Y();
373 mefistoToDS[m + 1] = (*itp).second;
374 // MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
375 // MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
382 gp_Pnt2d p = C2d->Value(l); // last point = Vertex Reversed
383 uvslf[m].x = scalex * p.X();
384 uvslf[m].y = scaley * p.Y();
385 mefistoToDS[m + 1] = idLast;
386 // MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
387 // MESSAGE("__ l "<<l<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
389 map < double, int >::reverse_iterator itp = params.rbegin();
390 for (int i = nbPoints; i >= 1; i--)
392 double param = (*itp).first;
393 gp_Pnt2d p = C2d->Value(param);
394 uvslf[m].x = scalex * p.X();
395 uvslf[m].y = scaley * p.Y();
396 mefistoToDS[m + 1] = (*itp).second;
397 // MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
398 // MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
406 //=============================================================================
410 //=============================================================================
412 // **** a mettre dans SMESH_Algo ou SMESH_2D_Algo
414 void SMESH_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh & aMesh,
415 const TopoDS_Face & aFace, double &scalex, double &scaley)
417 //MESSAGE("SMESH_MEFISTO_2D::ComputeScaleOnFace");
418 TopoDS_Face F = TopoDS::Face(aFace.Oriented(TopAbs_FORWARD));
419 TopoDS_Wire W = BRepTools::OuterWire(F);
421 BRepTools_WireExplorer wexp(W, F);
423 double xmin = 1.e300; // min & max of face 2D parametric coord.
424 double xmax = -1.e300;
425 double ymin = 1.e300;
426 double ymax = -1.e300;
430 for (wexp.Init(W, F); wexp.More(); wexp.Next())
432 const TopoDS_Edge & E = wexp.Current();
434 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
435 for (int i = 0; i <= nbp; i++)
437 double param = f + (double (i) / double (nbp))*(l - f);
438 gp_Pnt2d p = C2d->Value(param);
447 // MESSAGE(" "<< f<<" "<<l<<" "<<param<<" "<<xmin<<" "<<xmax<<" "<<ymin<<" "<<ymax);
454 double xmoy = (xmax + xmin) / 2.;
455 double ymoy = (ymax + ymin) / 2.;
457 Handle(Geom_Surface) S = BRep_Tool::Surface(F); // 3D surface
461 gp_Pnt PX0 = S->Value(xmin, ymoy);
462 gp_Pnt PY0 = S->Value(xmoy, ymin);
463 for (int i = 1; i <= nbp; i++)
465 double x = xmin + (double (i) / double (nbp))*(xmax - xmin);
466 gp_Pnt PX = S->Value(x, ymoy);
467 double y = ymin + (double (i) / double (nbp))*(ymax - ymin);
468 gp_Pnt PY = S->Value(xmoy, y);
469 length_x += PX.Distance(PX0);
470 length_y += PY.Distance(PY0);
471 PX0.SetCoord(PX.X(), PX.Y(), PX.Z());
472 PY0.SetCoord(PY.X(), PY.Y(), PY.Z());
476 scalex = length_x / (xmax - xmin);
477 scaley = length_y / (ymax - ymin);
484 //=============================================================================
488 //=============================================================================
490 void SMESH_MEFISTO_2D::StoreResult(SMESH_Mesh & aMesh,
491 Z nbst, R2 * uvst, Z nbt, Z * nust,
492 const TopoDS_Face & F, bool faceIsForward, map < int, int >&mefistoToDS)
496 ComputeScaleOnFace(aMesh, F, scalex, scaley);
498 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
501 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
503 for (n = 0; n < nbst; n++)
505 double u = uvst[n][0] / scalex;
506 double v = uvst[n][1] / scaley;
507 gp_Pnt P = S->Value(u, v);
509 if (mefistoToDS.find(n + 1) == mefistoToDS.end())
511 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
512 meshDS->SetNodeOnFace(node, F);
514 //MESSAGE(nodeId<<" "<<P.X()<<" "<<P.Y()<<" "<<P.Z());
515 mefistoToDS[n + 1] = node->GetID();
516 //MESSAGE(" "<<n<<" "<<mefistoToDS[n+1]);
517 SMDS_FacePosition* fpos
518 = static_cast<SMDS_FacePosition*>(node->GetPosition());
519 fpos->SetUParameter(u);
520 fpos->SetVParameter(v);
527 //SCRUTE(faceIsForward);
528 for (n = 1; n <= nbt; n++)
530 int inode1 = nust[m++];
531 int inode2 = nust[m++];
532 int inode3 = nust[m++];
534 int nodeId1 = mefistoToDS[inode1];
535 int nodeId2 = mefistoToDS[inode2];
536 int nodeId3 = mefistoToDS[inode3];
537 //MESSAGE("-- "<<inode1<<" "<<inode2<<" "<<inode3<<" ++ "<<nodeId1<<" "<<nodeId2<<" "<<nodeId3);
539 // triangle points must be in trigonometric order if face is Forward
540 // else they must be put clockwise
542 bool triangleIsWellOriented = faceIsForward;
544 SMDS_MeshElement * elt;
545 if (triangleIsWellOriented)
546 elt = meshDS->AddFace(nodeId1, nodeId2, nodeId3);
548 elt = meshDS->AddFace(nodeId1, nodeId3, nodeId2);
550 meshDS->SetMeshElementOnShape(elt, F);
555 //=============================================================================
559 //=============================================================================
561 double SMESH_MEFISTO_2D::ComputeEdgeElementLength(SMESH_Mesh & aMesh,
562 const TopoDS_Shape & aShape)
564 MESSAGE("SMESH_MEFISTO_2D::ComputeEdgeElementLength");
565 // **** a mettre dans SMESH_2D_Algo ?
567 const TopoDS_Face & FF = TopoDS::Face(aShape);
568 bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
569 TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
571 double meanElementLength = 100;
572 double wireLength = 0;
573 int wireElementsNumber = 0;
574 for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next())
576 const TopoDS_Wire & W = TopoDS::Wire(exp.Current());
577 for (TopExp_Explorer expe(W, TopAbs_EDGE); expe.More(); expe.Next())
579 const TopoDS_Edge & E = TopoDS::Edge(expe.Current());
580 int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
581 double length = EdgeLength(E);
582 wireLength += length;
583 wireElementsNumber += nb;
586 if (wireElementsNumber)
587 meanElementLength = wireLength / wireElementsNumber;
588 //SCRUTE(meanElementLength);
589 return meanElementLength;
592 //=============================================================================
596 //=============================================================================
598 ostream & SMESH_MEFISTO_2D::SaveTo(ostream & save)
603 //=============================================================================
607 //=============================================================================
609 istream & SMESH_MEFISTO_2D::LoadFrom(istream & load)
611 return load >> (*this);
614 //=============================================================================
618 //=============================================================================
620 ostream & operator <<(ostream & save, SMESH_MEFISTO_2D & hyp)
625 //=============================================================================
629 //=============================================================================
631 istream & operator >>(istream & load, SMESH_MEFISTO_2D & hyp)