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_Quadrangle_2D.cxx
25 // Moved here from SMESH_Quadrangle_2D.cxx
26 // Author : Paul RASCLE, EDF
31 #include "StdMeshers_Quadrangle_2D.hxx"
32 #include "SMESH_Gen.hxx"
33 #include "SMESH_Mesh.hxx"
35 #include "SMDS_MeshElement.hxx"
36 #include "SMDS_MeshNode.hxx"
37 #include "SMDS_EdgePosition.hxx"
38 #include "SMDS_FacePosition.hxx"
40 #include <BRep_Tool.hxx>
41 #include <BRepTools.hxx>
42 #include <BRepTools_WireExplorer.hxx>
43 #include <Geom_Surface.hxx>
44 #include <Geom_Curve.hxx>
45 #include <Geom2d_Curve.hxx>
46 #include <Handle_Geom2d_Curve.hxx>
47 #include <Handle_Geom_Curve.hxx>
48 #include <gp_Pnt2d.hxx>
49 #include <TColStd_ListIteratorOfListOfInteger.hxx>
51 #include "utilities.h"
52 #include "Utils_ExceptHandlers.hxx"
55 //=============================================================================
59 //=============================================================================
61 StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D(int hypId,
62 int studyId, SMESH_Gen * gen):SMESH_2D_Algo(hypId, studyId, gen)
64 MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D");
65 _name = "Quadrangle_2D";
66 // _shapeType = TopAbs_FACE;
67 _shapeType = (1 << TopAbs_FACE);
70 //=============================================================================
74 //=============================================================================
76 StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D()
78 MESSAGE("StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D");
81 //=============================================================================
85 //=============================================================================
87 bool StdMeshers_Quadrangle_2D::CheckHypothesis
89 const TopoDS_Shape& aShape,
90 SMESH_Hypothesis::Hypothesis_Status& aStatus)
92 //MESSAGE("StdMeshers_Quadrangle_2D::CheckHypothesis");
95 aStatus = SMESH_Hypothesis::HYP_OK;
102 //=============================================================================
106 //=============================================================================
108 bool StdMeshers_Quadrangle_2D::Compute(SMESH_Mesh & aMesh,
109 const TopoDS_Shape & aShape)throw(SALOME_Exception)
111 Unexpect aCatch(SalomeException);
112 //MESSAGE("StdMeshers_Quadrangle_2D::Compute");
113 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
114 SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape);
116 FaceQuadStruct *quad = CheckAnd2Dcompute(aMesh, aShape);
120 // --- compute 3D values on points, store points & quadrangles
122 int nbdown = quad->nbPts[0];
123 int nbright = quad->nbPts[1];
124 int nbVertices = nbdown * nbright;
125 int nbQuad = (nbdown - 1) * (nbright - 1);
126 //SCRUTE(nbVertices);
129 // const TopoDS_Face& FF = TopoDS::Face(aShape);
130 // bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
131 // TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
132 const TopoDS_Face & F = TopoDS::Face(aShape);
133 bool faceIsForward = (F.Orientation() == TopAbs_FORWARD);
134 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
136 for (int i = 1; i < nbdown - 1; i++)
137 for (int j = 1; j < nbright - 1; j++) // internal points
139 int ij = j * nbdown + i;
140 double u = quad->uv_grid[ij].u;
141 double v = quad->uv_grid[ij].v;
142 gp_Pnt P = S->Value(u, v);
143 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
144 meshDS->SetNodeOnFace(node, F);
145 quad->uv_grid[ij].node = node;
146 // Handle (SMDS_FacePosition) fpos
147 // = new SMDS_FacePosition(theSubMesh->GetId(),i,j); // easier than u,v
148 // node->SetPosition(fpos);
149 SMDS_FacePosition* fpos =
150 dynamic_cast<SMDS_FacePosition*>(node->GetPosition().get());
151 fpos->SetUParameter(i);
152 fpos->SetVParameter(j);
155 // bool isQuadForward = ( faceIsForward == quad->isEdgeForward[0]);
156 for (int i = 0; i < nbdown - 1; i++)
157 for (int j = 0; j < nbright - 1; j++) // faces
159 const SMDS_MeshNode *a, *b, *c, *d;
160 a = quad->uv_grid[j * nbdown + i].node;
161 b = quad->uv_grid[j * nbdown + i + 1].node;
162 c = quad->uv_grid[(j + 1) * nbdown + i + 1].node;
163 d = quad->uv_grid[(j + 1) * nbdown + i].node;
164 // if (isQuadForward) faceId = meshDS->AddFace(a,b,c,d);
165 // else faceId = meshDS->AddFace(a,d,c,b);
166 SMDS_MeshFace * face = meshDS->AddFace(a, b, c, d);
167 meshDS->SetMeshElementOnShape(face, F);
175 //=============================================================================
179 //=============================================================================
181 FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute(SMESH_Mesh & aMesh,
182 const TopoDS_Shape & aShape)throw(SALOME_Exception)
184 Unexpect aCatch(SalomeException);
185 //MESSAGE("StdMeshers_Quadrangle_2D::ComputeWithoutStore");
187 SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape);
189 // const TopoDS_Face& FF = TopoDS::Face(aShape);
190 // bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
191 // TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
192 const TopoDS_Face & F = TopoDS::Face(aShape);
193 bool faceIsForward = (F.Orientation() == TopAbs_FORWARD);
195 // verify 1 wire only, with 4 edges, same number of points on opposite edges
197 if (NumberOfWires(F) != 1)
199 MESSAGE("only 1 wire by face (quadrangles)");
201 //throw SALOME_Exception(LOCALIZED("only 1 wire by face (quadrangles)"));
203 // const TopoDS_Wire WW = BRepTools::OuterWire(F);
204 // TopoDS_Wire W = TopoDS::Wire(WW.Oriented(TopAbs_FORWARD));
205 const TopoDS_Wire & W = BRepTools::OuterWire(F);
206 BRepTools_WireExplorer wexp(W, F);
208 FaceQuadStruct *quad = new FaceQuadStruct;
209 for (int i = 0; i < 4; i++)
210 quad->uv_edges[i] = 0;
214 for (wexp.Init(W, F); wexp.More(); wexp.Next())
216 // const TopoDS_Edge& EE = wexp.Current();
217 // TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
218 const TopoDS_Edge & E = wexp.Current();
219 int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
222 quad->edge[nbEdges] = E;
223 quad->nbPts[nbEdges] = nb + 2; // internal points + 2 extrema
230 MESSAGE("face must have 4 edges /quadrangles");
233 //throw SALOME_Exception(LOCALIZED("face must have 4 edges /quadrangles"));
236 if (quad->nbPts[0] != quad->nbPts[2])
238 MESSAGE("different point number-opposed edge");
241 //throw SALOME_Exception(LOCALIZED("different point number-opposed edge"));
244 if (quad->nbPts[1] != quad->nbPts[3])
246 MESSAGE("different point number-opposed edge");
249 //throw SALOME_Exception(LOCALIZED("different point number-opposed edge"));
252 // set normalized grid on unit square in parametric domain
254 SetNormalizedGrid(aMesh, F, quad);
259 //=============================================================================
263 //=============================================================================
265 void StdMeshers_Quadrangle_2D::QuadDelete(FaceQuadStruct * quad)
267 //MESSAGE("StdMeshers_Quadrangle_2D::QuadDelete");
270 for (int i = 0; i < 4; i++)
272 if (quad->uv_edges[i])
273 delete[]quad->uv_edges[i];
274 quad->edge[i].Nullify();
277 delete[]quad->uv_grid;
282 //=============================================================================
286 //=============================================================================
288 void StdMeshers_Quadrangle_2D::SetNormalizedGrid(SMESH_Mesh & aMesh,
289 const TopoDS_Shape & aShape, FaceQuadStruct * quad) throw(SALOME_Exception)
291 Unexpect aCatch(SalomeException);
292 // Algorithme décrit dans "Génération automatique de maillages"
293 // P.L. GEORGE, MASSON, § 6.4.1 p. 84-85
294 // traitement dans le domaine paramétrique 2d u,v
295 // transport - projection sur le carré unité
297 const TopoDS_Face & F = TopoDS::Face(aShape);
299 // 1 --- find orientation of the 4 edges, by test on extrema
302 // |<----north-2-------^ a3 -------------> a2
304 // west-3 east-1 =right | |
308 // v----south-0--------> a0 -------------> a1
313 Handle(Geom2d_Curve) c2d[4];
316 for (int i = 0; i < 4; i++)
318 c2d[i] = BRep_Tool::CurveOnSurface(quad->edge[i],
319 F, quad->first[i], quad->last[i]);
320 pf[i] = c2d[i]->Value(quad->first[i]);
321 pl[i] = c2d[i]->Value(quad->last[i]);
322 quad->isEdgeForward[i] = false;
325 double eps2d = 1.e-3; // *** utiliser plutot TopExp::CommonVertex, puis
326 // distances si piece fausse
328 if ((pf[1].Distance(pl[0]) < eps2d) || (pl[1].Distance(pl[0]) < eps2d))
330 quad->isEdgeForward[0] = true;
334 double tmp = quad->first[0];
335 quad->first[0] = quad->last[0];
337 pf[0] = c2d[0]->Value(quad->first[0]);
338 pl[0] = c2d[0]->Value(quad->last[0]);
340 for (int i = 1; i < 4; i++)
342 quad->isEdgeForward[i] = (pf[i].Distance(pl[i - 1]) < eps2d);
343 if (!quad->isEdgeForward[i])
345 double tmp = quad->first[i];
346 quad->first[i] = quad->last[i];
348 pf[i] = c2d[i]->Value(quad->first[i]);
349 pl[i] = c2d[i]->Value(quad->last[i]);
350 //SCRUTE(pf[i].Distance(pl[i-1]));
351 ASSERT(pf[i].Distance(pl[i - 1]) < eps2d);
354 //SCRUTE(pf[0].Distance(pl[3]));
355 ASSERT(pf[0].Distance(pl[3]) < eps2d);
357 // for (int i=0; i<4; i++)
359 // SCRUTE(quad->isEdgeForward[i]);
360 // MESSAGE(" -first "<<i<<" "<<pf[i].X()<<" "<<pf[i].Y());
361 // MESSAGE(" -last "<<i<<" "<<pl[i].X()<<" "<<pl[i].Y());
364 // 2 --- load 2d edge points (u,v) with orientation and value on unit square
366 for (int i = 0; i < 2; i++)
368 quad->uv_edges[i] = LoadEdgePoints(aMesh, F,
369 quad->edge[i], quad->first[i], quad->last[i]);
371 // quad->isEdgeForward[i]);
373 for (int i = 2; i < 4; i++)
375 quad->uv_edges[i] = LoadEdgePoints(aMesh, F,
376 quad->edge[i], quad->last[i], quad->first[i]);
378 // !quad->isEdgeForward[i]);
381 // 3 --- 2D normalized values on unit square [0..1][0..1]
383 int nbdown = quad->nbPts[0];
384 int nbright = quad->nbPts[1];
385 quad->uv_grid = new UVPtStruct[nbright * nbdown];
387 UVPtStruct *uv_grid = quad->uv_grid;
388 UVPtStruct *uv_e0 = quad->uv_edges[0];
389 UVPtStruct *uv_e1 = quad->uv_edges[1];
390 UVPtStruct *uv_e2 = quad->uv_edges[2];
391 UVPtStruct *uv_e3 = quad->uv_edges[3];
400 for (int i = 0; i < nbdown; i++)
402 int ij = j * nbdown + i;
403 uv_grid[ij].node = uv_e0[i].node;
406 for (int j = 0; j < nbright; j++)
408 int ij = j * nbdown + i;
409 uv_grid[ij].node = uv_e1[j].node;
412 for (int i = 0; i < nbdown; i++)
414 int ij = j * nbdown + i;
415 uv_grid[ij].node = uv_e2[i].node;
418 for (int j = 0; j < nbright; j++)
420 int ij = j * nbdown + i;
421 uv_grid[ij].node = uv_e3[j].node;
424 // normalized 2d values on grid
426 for (int i = 0; i < nbdown; i++)
427 for (int j = 0; j < nbright; j++)
429 int ij = j * nbdown + i;
430 // --- droite i cste : x = x0 + y(x1-x0)
431 double x0 = uv_e0[i].normParam; // bas - sud
432 double x1 = uv_e2[i].normParam; // haut - nord
433 // --- droite j cste : y = y0 + x(y1-y0)
434 double y0 = uv_e3[j].normParam; // gauche-ouest
435 double y1 = uv_e1[j].normParam; // droite - est
436 // --- intersection : x=x0+(y0+x(y1-y0))(x1-x0)
437 double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
438 double y = y0 + x * (y1 - y0);
441 //MESSAGE("-xy-01 "<<x0<<" "<<x1<<" "<<y0<<" "<<y1);
442 //MESSAGE("-xy-norm "<<i<<" "<<j<<" "<<x<<" "<<y);
445 // 4 --- projection on 2d domain (u,v)
447 for (int i = 0; i < nbdown; i++)
448 for (int j = 0; j < nbright; j++)
450 int ij = j * nbdown + i;
451 double x = uv_grid[ij].x;
452 double y = uv_grid[ij].y;
453 double param_0 = uv_e0[0].param + x * (uv_e0[nbdown - 1].param - uv_e0[0].param); // sud
454 double param_2 = uv_e2[0].param + x * (uv_e2[nbdown - 1].param - uv_e2[0].param); // nord
455 double param_1 = uv_e1[0].param + y * (uv_e1[nbright - 1].param - uv_e1[0].param); // est
456 double param_3 = uv_e3[0].param + y * (uv_e3[nbright - 1].param - uv_e3[0].param); // ouest
458 //MESSAGE("params "<<param_0<<" "<<param_1<<" "<<param_2<<" "<<param_3);
459 gp_Pnt2d p0 = c2d[0]->Value(param_0);
460 gp_Pnt2d p1 = c2d[1]->Value(param_1);
461 gp_Pnt2d p2 = c2d[2]->Value(param_2);
462 gp_Pnt2d p3 = c2d[3]->Value(param_3);
465 (1 - y) * p0.X() + x * p1.X() + y * p2.X() + (1 - x) * p3.X();
467 (1 - y) * p0.Y() + x * p1.Y() + y * p2.Y() + (1 - x) * p3.Y();
469 u -= (1 - x) * (1 - y) * a0.X() + x * (1 - y) * a1.X() +
470 x * y * a2.X() + (1 - x) * y * a3.X();
471 v -= (1 - x) * (1 - y) * a0.Y() + x * (1 - y) * a1.Y() +
472 x * y * a2.Y() + (1 - x) * y * a3.Y();
477 //MESSAGE("-uv- "<<i<<" "<<j<<" "<<uv_grid[ij].u<<" "<<uv_grid[ij].v);
481 //=============================================================================
485 //=============================================================================
487 UVPtStruct *StdMeshers_Quadrangle_2D::LoadEdgePoints(SMESH_Mesh & aMesh,
488 const TopoDS_Face & F, const TopoDS_Edge & E, double first, double last)
491 //MESSAGE("StdMeshers_Quadrangle_2D::LoadEdgePoints");
493 SMDS_Mesh * meshDS = aMesh.GetMeshDS();
495 // --- IDNodes of first and last Vertex
497 TopoDS_Vertex VFirst, VLast;
498 TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
500 ASSERT(!VFirst.IsNull());
501 SMDS_NodeIteratorPtr lid= aMesh.GetSubMesh(VFirst)->GetSubMeshDS()->GetNodes();
502 const SMDS_MeshNode * idFirst = lid->next();
504 ASSERT(!VLast.IsNull());
505 lid=aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes();
506 const SMDS_MeshNode * idLast = lid->next();
508 // --- edge internal IDNodes (relies on good order storage, not checked)
510 int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
512 UVPtStruct *uvslf = new UVPtStruct[nbPoints + 2];
515 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
517 map<double, const SMDS_MeshNode *> params;
518 SMDS_NodeIteratorPtr ite= aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes();
522 const SMDS_MeshNode * node = ite->next();
523 const SMDS_EdgePosition* epos =
524 static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
525 double param = epos->GetUParameter();
526 params[param] = node;
529 bool isForward = (((l - f) * (last - first)) > 0);
536 gp_Pnt2d p = C2d->Value(f); // first point = Vertex Forward
540 uvslf[0].node = idFirst;
541 //MESSAGE("__ f "<<f<<" "<<uvslf[0].x <<" "<<uvslf[0].y);
542 map < double, const SMDS_MeshNode* >::iterator itp = params.begin();
543 for (int i = 1; i <= nbPoints; i++) // nbPoints internal
545 double param = (*itp).first;
546 gp_Pnt2d p = C2d->Value(param);
549 uvslf[i].param = param;
550 uvslf[i].node = (*itp).second;
551 //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[i].x <<" "<<uvslf[i].y);
554 p = C2d->Value(l); // last point = Vertex Reversed
555 uvslf[nbPoints + 1].x = p.X();
556 uvslf[nbPoints + 1].y = p.Y();
557 uvslf[nbPoints + 1].param = l;
558 uvslf[nbPoints + 1].node = idLast;
559 //MESSAGE("__ l "<<l<<" "<<uvslf[nbPoints+1].x <<" "<<uvslf[nbPoints+1].y);
565 gp_Pnt2d p = C2d->Value(l); // first point = Vertex Reversed
569 uvslf[0].node = idLast;
570 //MESSAGE("__ l "<<l<<" "<<uvslf[0].x <<" "<<uvslf[0].y);
571 map < double, const SMDS_MeshNode* >::reverse_iterator itp = params.rbegin();
572 for (int j = nbPoints; j >= 1; j--) // nbPoints internal
574 double param = (*itp).first;
575 int i = nbPoints + 1 - j;
576 gp_Pnt2d p = C2d->Value(param);
579 uvslf[i].param = param;
580 uvslf[i].node = (*itp).second;
581 //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[i].x <<" "<<uvslf[i].y);
584 p = C2d->Value(f); // last point = Vertex Forward
585 uvslf[nbPoints + 1].x = p.X();
586 uvslf[nbPoints + 1].y = p.Y();
587 uvslf[nbPoints + 1].param = f;
588 uvslf[nbPoints + 1].node = idFirst;
589 //MESSAGE("__ f "<<f<<" "<<uvslf[nbPoints+1].x <<" "<<uvslf[nbPoints+1].y);
592 ASSERT(paramin != paramax);
593 for (int i = 0; i < nbPoints + 2; i++)
595 uvslf[i].normParam = (uvslf[i].param - paramin) / (paramax - paramin);
596 //SCRUTE(uvslf[i].normParam);
602 //=============================================================================
606 //=============================================================================
608 ostream & StdMeshers_Quadrangle_2D::SaveTo(ostream & save)
613 //=============================================================================
617 //=============================================================================
619 istream & StdMeshers_Quadrangle_2D::LoadFrom(istream & load)
624 //=============================================================================
628 //=============================================================================
630 ostream & operator <<(ostream & save, StdMeshers_Quadrangle_2D & hyp)
632 return hyp.SaveTo( save );
635 //=============================================================================
639 //=============================================================================
641 istream & operator >>(istream & load, StdMeshers_Quadrangle_2D & hyp)
643 return hyp.LoadFrom( load );