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"
34 #include "SMESH_subMesh.hxx"
36 #include "SMDS_MeshElement.hxx"
37 #include "SMDS_MeshNode.hxx"
38 #include "SMDS_EdgePosition.hxx"
39 #include "SMDS_FacePosition.hxx"
41 #include <BRep_Tool.hxx>
42 #include <BRepTools.hxx>
43 #include <BRepTools_WireExplorer.hxx>
45 #include <Geom_Surface.hxx>
46 #include <Geom_Curve.hxx>
47 #include <Geom2d_Curve.hxx>
48 #include <GeomAdaptor_Curve.hxx>
49 #include <GCPnts_UniformAbscissa.hxx>
51 #include <Precision.hxx>
52 #include <gp_Pnt2d.hxx>
53 #include <TColStd_ListIteratorOfListOfInteger.hxx>
55 #include "utilities.h"
56 #include "Utils_ExceptHandlers.hxx"
59 //=============================================================================
63 //=============================================================================
65 StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, SMESH_Gen* gen)
66 : SMESH_2D_Algo(hypId, studyId, gen)
68 MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D");
69 _name = "Quadrangle_2D";
70 // _shapeType = TopAbs_FACE;
71 _shapeType = (1 << TopAbs_FACE);
74 //=============================================================================
78 //=============================================================================
80 StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D()
82 MESSAGE("StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D");
85 //=============================================================================
89 //=============================================================================
91 bool StdMeshers_Quadrangle_2D::CheckHypothesis
93 const TopoDS_Shape& aShape,
94 SMESH_Hypothesis::Hypothesis_Status& aStatus)
96 //MESSAGE("StdMeshers_Quadrangle_2D::CheckHypothesis");
99 aStatus = SMESH_Hypothesis::HYP_OK;
106 //=============================================================================
110 //=============================================================================
112 bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh,
113 const TopoDS_Shape& aShape) throw (SALOME_Exception)
115 Unexpect aCatch(SalomeException);
116 //MESSAGE("StdMeshers_Quadrangle_2D::Compute");
117 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
118 aMesh.GetSubMesh(aShape);
120 FaceQuadStruct *quad = CheckAnd2Dcompute(aMesh, aShape);
124 // --- compute 3D values on points, store points & quadrangles
126 int nbdown = quad->nbPts[0];
127 int nbup = quad->nbPts[2];
128 // bool isDownOut = (nbdown > nbup);
129 // bool isUpOut = (nbdown < nbup);
131 int nbright = quad->nbPts[1];
132 int nbleft = quad->nbPts[3];
133 // bool isRightOut = (nbright > nbleft);
134 // bool isLeftOut = (nbright < nbleft);
136 int nbhoriz = Min(nbdown, nbup);
137 int nbvertic = Min(nbright, nbleft);
139 const TopoDS_Face& F = TopoDS::Face(aShape);
140 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
142 // internal mesh nodes
144 for (i = 1; i < nbhoriz - 1; i++) {
145 for (j = 1; j < nbvertic - 1; j++) {
146 int ij = j * nbhoriz + i;
147 double u = quad->uv_grid[ij].u;
148 double v = quad->uv_grid[ij].v;
149 gp_Pnt P = S->Value(u, v);
150 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
151 meshDS->SetNodeOnFace(node, F);
152 quad->uv_grid[ij].node = node;
153 SMDS_FacePosition* fpos =
154 dynamic_cast<SMDS_FacePosition*>(node->GetPosition().get());
155 fpos->SetUParameter(u);
156 fpos->SetVParameter(v);
163 // --.--.--.--.--.-- nbvertic
169 // ---.----.----.--- 0
170 // 0 > > > > > > > > nbhoriz
176 int iup = nbhoriz - 1;
177 if (quad->isEdgeOut[3]) { ilow++; } else { if (quad->isEdgeOut[1]) iup--; }
180 int jup = nbvertic - 1;
181 if (quad->isEdgeOut[0]) { jlow++; } else { if (quad->isEdgeOut[2]) jup--; }
183 // regular quadrangles
184 // bool isQuadForward = ( faceIsForward == quad->isEdgeForward[0]);
185 for (i = ilow; i < iup; i++) {
186 for (j = jlow; j < jup; j++) {
187 const SMDS_MeshNode *a, *b, *c, *d;
188 a = quad->uv_grid[j * nbhoriz + i].node;
189 b = quad->uv_grid[j * nbhoriz + i + 1].node;
190 c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node;
191 d = quad->uv_grid[(j + 1) * nbhoriz + i].node;
192 // if (isQuadForward) faceId = meshDS->AddFace(a,b,c,d);
193 // else faceId = meshDS->AddFace(a,d,c,b);
194 SMDS_MeshFace * face = meshDS->AddFace(a, b, c, d);
195 meshDS->SetMeshElementOnShape(face, F);
199 UVPtStruct *uv_e0 = quad->uv_edges[0];
200 UVPtStruct *uv_e1 = quad->uv_edges[1];
201 UVPtStruct *uv_e2 = quad->uv_edges[2];
202 UVPtStruct *uv_e3 = quad->uv_edges[3];
204 double eps = Precision::Confusion();
206 // Boundary quadrangles
208 if (quad->isEdgeOut[0]) {
211 // |___|___|___|___|___|___|
213 // |___|___|___|___|___|___|
215 // |___|___|___|___|___|___| __ first row of the regular grid
216 // . . . . . . . . . __ down edge nodes
218 // >->->->->->->->->->->->-> -- direction of processing
220 int g = 0; // number of last processed node in the regular grid
222 // number of last node of the down edge to be processed
223 int stop = nbdown - 1;
224 // if right edge is out, we will stop at a node, previous to the last one
225 if (quad->isEdgeOut[1]) stop--;
227 // for each node of the down edge find nearest node
228 // in the first row of the regular grid and link them
229 for (i = 0; i < stop; i++) {
230 const SMDS_MeshNode *a, *b, *c, *d;
232 b = uv_e0[i + 1].node;
233 gp_Pnt pb (b->X(), b->Y(), b->Z());
235 // find node c in the regular grid, which will be linked with node b
238 // right bound reached, link with the rightmost node
240 c = quad->uv_grid[nbhoriz + iup].node;
242 // find in the grid node c, nearest to the b
243 double mind = RealLast();
244 for (int k = g; k <= iup; k++) {
246 const SMDS_MeshNode *nk;
247 if (k < ilow) // this can be, if left edge is out
248 nk = uv_e3[1].node; // get node from the left edge
250 nk = quad->uv_grid[nbhoriz + k].node; // get one of middle nodes
252 gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
253 double dist = pb.Distance(pnk);
254 if (dist < mind - eps) {
264 if (near == g) { // make triangle
265 SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
266 meshDS->SetMeshElementOnShape(face, F);
267 } else { // make quadrangle
271 d = quad->uv_grid[nbhoriz + near - 1].node;
272 SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
273 meshDS->SetMeshElementOnShape(face, F);
275 // if node d is not at position g - make additional triangles
277 for (int k = near - 1; k > g; k--) {
278 c = quad->uv_grid[nbhoriz + k].node;
282 d = quad->uv_grid[nbhoriz + k - 1].node;
283 SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
284 meshDS->SetMeshElementOnShape(face, F);
291 if (quad->isEdgeOut[2]) {
294 // <-<-<-<-<-<-<-<-<-<-<-<-< -- direction of processing
296 // . . . . . . . . . __ up edge nodes
297 // ___ ___ ___ ___ ___ ___ __ first row of the regular grid
299 // |___|___|___|___|___|___|
301 // |___|___|___|___|___|___|
304 int g = nbhoriz - 1; // last processed node in the regular grid
307 // if left edge is out, we will stop at a second node
308 if (quad->isEdgeOut[3]) stop++;
310 // for each node of the up edge find nearest node
311 // in the first row of the regular grid and link them
312 for (i = nbup - 1; i > stop; i--) {
313 const SMDS_MeshNode *a, *b, *c, *d;
315 b = uv_e2[i - 1].node;
316 gp_Pnt pb (b->X(), b->Y(), b->Z());
318 // find node c in the grid, which will be linked with node b
320 if (i == stop + 1) { // left bound reached, link with the leftmost node
321 c = quad->uv_grid[nbhoriz*(nbvertic - 2) + ilow].node;
324 // find node c in the grid, nearest to the b
325 double mind = RealLast();
326 for (int k = g; k >= ilow; k--) {
327 const SMDS_MeshNode *nk;
329 nk = uv_e1[nbright - 2].node;
331 nk = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
332 gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
333 double dist = pb.Distance(pnk);
334 if (dist < mind - eps) {
344 if (near == g) { // make triangle
345 SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
346 meshDS->SetMeshElementOnShape(face, F);
347 } else { // make quadrangle
349 d = uv_e1[nbright - 2].node;
351 d = quad->uv_grid[nbhoriz*(nbvertic - 2) + near + 1].node;
352 SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
353 meshDS->SetMeshElementOnShape(face, F);
355 if (near + 1 < g) { // if d not is at g - make additional triangles
356 for (int k = near + 1; k < g; k++) {
357 c = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
359 d = uv_e1[nbright - 2].node;
361 d = quad->uv_grid[nbhoriz*(nbvertic - 2) + k + 1].node;
362 SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
363 meshDS->SetMeshElementOnShape(face, F);
372 // right or left boundary quadrangles
373 if (quad->isEdgeOut[1]) {
374 // MESSAGE("right edge is out");
375 int g = 0; // last processed node in the grid
376 int stop = nbright - 1;
377 if (quad->isEdgeOut[2]) stop--;
378 for (i = 0; i < stop; i++) {
379 const SMDS_MeshNode *a, *b, *c, *d;
381 b = uv_e1[i + 1].node;
382 gp_Pnt pb (b->X(), b->Y(), b->Z());
384 // find node c in the grid, nearest to the b
386 if (i == stop - 1) { // up bondary reached
387 c = quad->uv_grid[nbhoriz*(jup + 1) - 2].node;
390 double mind = RealLast();
391 for (int k = g; k <= jup; k++) {
392 const SMDS_MeshNode *nk;
394 nk = uv_e0[nbdown - 2].node;
396 nk = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
397 gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
398 double dist = pb.Distance(pnk);
399 if (dist < mind - eps) {
409 if (near == g) { // make triangle
410 SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
411 meshDS->SetMeshElementOnShape(face, F);
412 } else { // make quadrangle
414 d = uv_e0[nbdown - 2].node;
416 d = quad->uv_grid[nbhoriz*near - 2].node;
417 SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
418 meshDS->SetMeshElementOnShape(face, F);
420 if (near - 1 > g) { // if d not is at g - make additional triangles
421 for (int k = near - 1; k > g; k--) {
422 c = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
424 d = uv_e0[nbdown - 2].node;
426 d = quad->uv_grid[nbhoriz*k - 2].node;
427 SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
428 meshDS->SetMeshElementOnShape(face, F);
435 if (quad->isEdgeOut[3]) {
436 // MESSAGE("left edge is out");
437 int g = nbvertic - 1; // last processed node in the grid
439 if (quad->isEdgeOut[0]) stop++;
440 for (i = nbleft - 1; i > stop; i--) {
441 const SMDS_MeshNode *a, *b, *c, *d;
443 b = uv_e3[i - 1].node;
444 gp_Pnt pb (b->X(), b->Y(), b->Z());
446 // find node c in the grid, nearest to the b
448 if (i == stop + 1) { // down bondary reached
449 c = quad->uv_grid[nbhoriz*jlow + 1].node;
452 double mind = RealLast();
453 for (int k = g; k >= jlow; k--) {
454 const SMDS_MeshNode *nk;
458 nk = quad->uv_grid[nbhoriz*k + 1].node;
459 gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
460 double dist = pb.Distance(pnk);
461 if (dist < mind - eps) {
471 if (near == g) { // make triangle
472 SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
473 meshDS->SetMeshElementOnShape(face, F);
474 } else { // make quadrangle
478 d = quad->uv_grid[nbhoriz*(near + 1) + 1].node;
479 SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
480 meshDS->SetMeshElementOnShape(face, F);
482 if (near + 1 < g) { // if d not is at g - make additional triangles
483 for (int k = near + 1; k < g; k++) {
484 c = quad->uv_grid[nbhoriz*k + 1].node;
488 d = quad->uv_grid[nbhoriz*(k + 1) + 1].node;
489 SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
490 meshDS->SetMeshElementOnShape(face, F);
504 //=============================================================================
508 //=============================================================================
510 FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute
511 (SMESH_Mesh & aMesh, const TopoDS_Shape & aShape) throw(SALOME_Exception)
513 Unexpect aCatch(SalomeException);
515 //SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape);
517 const TopoDS_Face & F = TopoDS::Face(aShape);
519 // verify 1 wire only, with 4 edges
521 if (NumberOfWires(F) != 1)
523 INFOS("only 1 wire by face (quadrangles)");
526 const TopoDS_Wire& W = BRepTools::OuterWire(F);
527 BRepTools_WireExplorer wexp (W, F);
529 FaceQuadStruct *quad = new FaceQuadStruct;
530 for (int i = 0; i < 4; i++)
531 quad->uv_edges[i] = 0;
535 for (wexp.Init(W, F); wexp.More(); wexp.Next())
537 const TopoDS_Edge& E = wexp.Current();
538 int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
541 quad->edge[nbEdges] = E;
542 quad->nbPts[nbEdges] = nb + 2; // internal points + 2 extrema
549 INFOS("face must have 4 edges /quadrangles");
554 // set normalized grid on unit square in parametric domain
556 SetNormalizedGrid(aMesh, F, quad);
561 //=============================================================================
565 //=============================================================================
567 void StdMeshers_Quadrangle_2D::QuadDelete (FaceQuadStruct * quad)
569 //MESSAGE("StdMeshers_Quadrangle_2D::QuadDelete");
572 for (int i = 0; i < 4; i++)
574 if (quad->uv_edges[i])
575 delete [] quad->uv_edges[i];
576 quad->edge[i].Nullify();
579 delete [] quad->uv_grid;
584 //=============================================================================
588 //=============================================================================
590 void StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh,
591 const TopoDS_Shape& aShape,
592 FaceQuadStruct* quad) throw (SALOME_Exception)
594 Unexpect aCatch(SalomeException);
595 // Algorithme décrit dans "Génération automatique de maillages"
596 // P.L. GEORGE, MASSON, § 6.4.1 p. 84-85
597 // traitement dans le domaine paramétrique 2d u,v
598 // transport - projection sur le carré unité
600 // MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid");
601 const TopoDS_Face& F = TopoDS::Face(aShape);
603 // 1 --- find orientation of the 4 edges, by test on extrema
606 // |<----north-2-------^ a3 -------------> a2
608 // west-3 east-1 =right | |
612 // v----south-0--------> a0 -------------> a1
617 Handle(Geom2d_Curve) c2d[4];
620 for (int i = 0; i < 4; i++)
622 c2d[i] = BRep_Tool::CurveOnSurface(quad->edge[i], F,
623 quad->first[i], quad->last[i]);
624 pf[i] = c2d[i]->Value(quad->first[i]);
625 pl[i] = c2d[i]->Value(quad->last[i]);
626 quad->isEdgeForward[i] = false;
629 double l0f1 = pl[0].SquareDistance(pf[1]);
630 double l0l1 = pl[0].SquareDistance(pl[1]);
631 double f0f1 = pf[0].SquareDistance(pf[1]);
632 double f0l1 = pf[0].SquareDistance(pl[1]);
633 if ( Min( l0f1, l0l1 ) < Min ( f0f1, f0l1 ))
635 quad->isEdgeForward[0] = true;
637 double tmp = quad->first[0];
638 quad->first[0] = quad->last[0];
640 pf[0] = c2d[0]->Value(quad->first[0]);
641 pl[0] = c2d[0]->Value(quad->last[0]);
643 for (int i = 1; i < 4; i++)
645 l0l1 = pl[i - 1].SquareDistance(pl[i]);
646 l0f1 = pl[i - 1].SquareDistance(pf[i]);
647 quad->isEdgeForward[i] = ( l0f1 < l0l1 );
648 if (!quad->isEdgeForward[i])
650 double tmp = quad->first[i];
651 quad->first[i] = quad->last[i];
653 pf[i] = c2d[i]->Value(quad->first[i]);
654 pl[i] = c2d[i]->Value(quad->last[i]);
658 // 2 --- load 2d edge points (u,v) with orientation and value on unit square
661 for (int i = 0; i < 2; i++)
663 quad->uv_edges[i] = LoadEdgePoints(aMesh, F, quad->edge[i],
664 quad->first[i], quad->last[i]);
665 if (!quad->uv_edges[i]) loadOk = false;
668 for (int i = 2; i < 4; i++)
670 quad->uv_edges[i] = LoadEdgePoints(aMesh, F, quad->edge[i],
671 quad->last[i], quad->first[i]);
672 if (!quad->uv_edges[i]) loadOk = false;
677 INFOS("StdMeshers_Quadrangle_2D::SetNormalizedGrid - LoadEdgePoints failed");
682 // 3 --- 2D normalized values on unit square [0..1][0..1]
684 int nbhoriz = Min(quad->nbPts[0], quad->nbPts[2]);
685 int nbvertic = Min(quad->nbPts[1], quad->nbPts[3]);
687 quad->isEdgeOut[0] = (quad->nbPts[0] > quad->nbPts[2]);
688 quad->isEdgeOut[1] = (quad->nbPts[1] > quad->nbPts[3]);
689 quad->isEdgeOut[2] = (quad->nbPts[2] > quad->nbPts[0]);
690 quad->isEdgeOut[3] = (quad->nbPts[3] > quad->nbPts[1]);
692 quad->uv_grid = new UVPtStruct[nbvertic * nbhoriz];
694 UVPtStruct *uv_grid = quad->uv_grid;
695 UVPtStruct *uv_e0 = quad->uv_edges[0];
696 UVPtStruct *uv_e1 = quad->uv_edges[1];
697 UVPtStruct *uv_e2 = quad->uv_edges[2];
698 UVPtStruct *uv_e3 = quad->uv_edges[3];
700 // nodes Id on "in" edges
701 if (! quad->isEdgeOut[0]) {
703 for (int i = 0; i < nbhoriz; i++) { // down
704 int ij = j * nbhoriz + i;
705 uv_grid[ij].node = uv_e0[i].node;
708 if (! quad->isEdgeOut[1]) {
710 for (int j = 0; j < nbvertic; j++) { // right
711 int ij = j * nbhoriz + i;
712 uv_grid[ij].node = uv_e1[j].node;
715 if (! quad->isEdgeOut[2]) {
716 int j = nbvertic - 1;
717 for (int i = 0; i < nbhoriz; i++) { // up
718 int ij = j * nbhoriz + i;
719 uv_grid[ij].node = uv_e2[i].node;
722 if (! quad->isEdgeOut[3]) {
724 for (int j = 0; j < nbvertic; j++) { // left
725 int ij = j * nbhoriz + i;
726 uv_grid[ij].node = uv_e3[j].node;
730 // falsificate "out" edges
731 if (quad->isEdgeOut[0]) // down
732 uv_e0 = MakeEdgePoints
733 (aMesh, F, quad->edge[0], quad->first[0], quad->last[0], nbhoriz - 1);
734 else if (quad->isEdgeOut[2]) // up
735 uv_e2 = MakeEdgePoints
736 (aMesh, F, quad->edge[2], quad->last[2], quad->first[2], nbhoriz - 1);
738 if (quad->isEdgeOut[1]) // right
739 uv_e1 = MakeEdgePoints
740 (aMesh, F, quad->edge[1], quad->first[1], quad->last[1], nbvertic - 1);
741 else if (quad->isEdgeOut[3]) // left
742 uv_e3 = MakeEdgePoints
743 (aMesh, F, quad->edge[3], quad->last[3], quad->first[3], nbvertic - 1);
745 // normalized 2d values on grid
746 for (int i = 0; i < nbhoriz; i++)
748 for (int j = 0; j < nbvertic; j++)
750 int ij = j * nbhoriz + i;
751 // --- droite i cste : x = x0 + y(x1-x0)
752 double x0 = uv_e0[i].normParam; // bas - sud
753 double x1 = uv_e2[i].normParam; // haut - nord
754 // --- droite j cste : y = y0 + x(y1-y0)
755 double y0 = uv_e3[j].normParam; // gauche-ouest
756 double y1 = uv_e1[j].normParam; // droite - est
757 // --- intersection : x=x0+(y0+x(y1-y0))(x1-x0)
758 double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
759 double y = y0 + x * (y1 - y0);
762 //MESSAGE("-xy-01 "<<x0<<" "<<x1<<" "<<y0<<" "<<y1);
763 //MESSAGE("-xy-norm "<<i<<" "<<j<<" "<<x<<" "<<y);
767 // 4 --- projection on 2d domain (u,v)
773 for (int i = 0; i < nbhoriz; i++)
775 for (int j = 0; j < nbvertic; j++)
777 int ij = j * nbhoriz + i;
778 double x = uv_grid[ij].x;
779 double y = uv_grid[ij].y;
780 double param_0 = uv_e0[0].param + x * (uv_e0[nbhoriz - 1].param - uv_e0[0].param); // sud
781 double param_2 = uv_e2[0].param + x * (uv_e2[nbhoriz - 1].param - uv_e2[0].param); // nord
782 double param_1 = uv_e1[0].param + y * (uv_e1[nbvertic - 1].param - uv_e1[0].param); // est
783 double param_3 = uv_e3[0].param + y * (uv_e3[nbvertic - 1].param - uv_e3[0].param); // ouest
785 //MESSAGE("params "<<param_0<<" "<<param_1<<" "<<param_2<<" "<<param_3);
786 gp_Pnt2d p0 = c2d[0]->Value(param_0);
787 gp_Pnt2d p1 = c2d[1]->Value(param_1);
788 gp_Pnt2d p2 = c2d[2]->Value(param_2);
789 gp_Pnt2d p3 = c2d[3]->Value(param_3);
791 double u = (1 - y) * p0.X() + x * p1.X() + y * p2.X() + (1 - x) * p3.X();
792 double v = (1 - y) * p0.Y() + x * p1.Y() + y * p2.Y() + (1 - x) * p3.Y();
794 u -= (1 - x) * (1 - y) * a0.X() + x * (1 - y) * a1.X() +
795 x * y * a2.X() + (1 - x) * y * a3.X();
796 v -= (1 - x) * (1 - y) * a0.Y() + x * (1 - y) * a1.Y() +
797 x * y * a2.Y() + (1 - x) * y * a3.Y();
805 //=============================================================================
809 //=============================================================================
810 UVPtStruct* StdMeshers_Quadrangle_2D::LoadEdgePoints (SMESH_Mesh & aMesh,
811 const TopoDS_Face& F,
812 const TopoDS_Edge& E,
813 double first, double last)
816 //MESSAGE("StdMeshers_Quadrangle_2D::LoadEdgePoints");
818 //SMDS_Mesh* meshDS = aMesh.GetMeshDS();
820 // --- IDNodes of first and last Vertex
822 TopoDS_Vertex VFirst, VLast;
823 TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
825 ASSERT(!VFirst.IsNull());
826 SMDS_NodeIteratorPtr lid = aMesh.GetSubMesh(VFirst)->GetSubMeshDS()->GetNodes();
829 MESSAGE ( "NO NODE BUILT ON VERTEX" );
832 const SMDS_MeshNode* idFirst = lid->next();
834 ASSERT(!VLast.IsNull());
835 lid = aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes();
838 MESSAGE ( "NO NODE BUILT ON VERTEX" );
841 const SMDS_MeshNode* idLast = lid->next();
843 // --- edge internal IDNodes (relies on good order storage, not checked)
845 map<double, const SMDS_MeshNode *> params;
846 SMDS_NodeIteratorPtr ite = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes();
850 const SMDS_MeshNode* node = ite->next();
851 const SMDS_EdgePosition* epos =
852 static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
853 double param = epos->GetUParameter();
854 params[param] = node;
857 int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
858 if (nbPoints != params.size())
860 MESSAGE( "BAD NODE ON EDGE POSITIONS" );
863 UVPtStruct* uvslf = new UVPtStruct[nbPoints + 2];
866 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
868 bool isForward = (((l - f) * (last - first)) > 0);
875 gp_Pnt2d p = C2d->Value(f); // first point = Vertex Forward
879 uvslf[0].node = idFirst;
880 //MESSAGE("__ f "<<f<<" "<<uvslf[0].x <<" "<<uvslf[0].y);
881 map < double, const SMDS_MeshNode* >::iterator itp = params.begin();
882 for (int i = 1; i <= nbPoints; i++) // nbPoints internal
884 double param = (*itp).first;
885 gp_Pnt2d p = C2d->Value(param);
888 uvslf[i].param = param;
889 uvslf[i].node = (*itp).second;
890 //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[i].x <<" "<<uvslf[i].y);
893 p = C2d->Value(l); // last point = Vertex Reversed
894 uvslf[nbPoints + 1].x = p.X();
895 uvslf[nbPoints + 1].y = p.Y();
896 uvslf[nbPoints + 1].param = l;
897 uvslf[nbPoints + 1].node = idLast;
898 //MESSAGE("__ l "<<l<<" "<<uvslf[nbPoints+1].x <<" "<<uvslf[nbPoints+1].y);
903 gp_Pnt2d p = C2d->Value(l); // first point = Vertex Reversed
907 uvslf[0].node = idLast;
908 //MESSAGE("__ l "<<l<<" "<<uvslf[0].x <<" "<<uvslf[0].y);
909 map < double, const SMDS_MeshNode* >::reverse_iterator itp = params.rbegin();
911 for (int j = nbPoints; j >= 1; j--) // nbPoints internal
913 double param = (*itp).first;
914 int i = nbPoints + 1 - j;
915 gp_Pnt2d p = C2d->Value(param);
918 uvslf[i].param = param;
919 uvslf[i].node = (*itp).second;
920 //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[i].x <<" "<<uvslf[i].y);
923 p = C2d->Value(f); // last point = Vertex Forward
924 uvslf[nbPoints + 1].x = p.X();
925 uvslf[nbPoints + 1].y = p.Y();
926 uvslf[nbPoints + 1].param = f;
927 uvslf[nbPoints + 1].node = idFirst;
928 //MESSAGE("__ f "<<f<<" "<<uvslf[nbPoints+1].x <<" "<<uvslf[nbPoints+1].y);
931 ASSERT(paramin != paramax);
932 for (int i = 0; i < nbPoints + 2; i++)
934 uvslf[i].normParam = (uvslf[i].param - paramin) / (paramax - paramin);
940 //=============================================================================
944 //=============================================================================
945 UVPtStruct* StdMeshers_Quadrangle_2D::MakeEdgePoints (SMESH_Mesh & aMesh,
946 const TopoDS_Face& F,
947 const TopoDS_Edge& E,
948 double first, double last,
951 // MESSAGE("StdMeshers_Quadrangle_2D::MakeEdgePoints");
953 UVPtStruct* uvslf = new UVPtStruct[nb_segm + 1];
956 // --- edge internal points
958 Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, fi, li);
959 if (!Curve.IsNull()) {
961 GeomAdaptor_Curve C3d (Curve);
962 double length = EdgeLength(E);
963 double eltSize = length / nb_segm;
964 GCPnts_UniformAbscissa Discret (C3d, eltSize, fi, li);
965 if (!Discret.IsDone()) return false;
966 int NbPoints = Discret.NbPoints();
967 for (int i = 1; i <= NbPoints; i++) {
968 double param = Discret.Parameter(i);
969 params.push_back(param);
972 catch (Standard_Failure) {
978 // Edge is a degenerated Edge
979 BRep_Tool::Range(E, fi, li);
980 double du = (li - fi) / nb_segm;
981 for (int i = 1; i <= nb_segm + 1; i++)
983 double param = fi + (i - 1) * du;
984 params.push_back(param);
989 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
992 bool isForward = (((l - f) * (last - first)) > 0);
994 list<double>::iterator itU = params.begin();
995 for (int i = 0; i <= nb_segm; i++) // nbPoints internal
998 gp_Pnt2d p = C2d->Value(param);
1001 uvslf[i].param = param;
1002 uvslf[i].normParam = (param - f) / (l - f);
1006 list<double>::reverse_iterator itU = params.rbegin();
1007 for (int j = nb_segm; j >= 0; j--) // nbPoints internal
1009 double param = *itU;
1010 int i = nb_segm - j;
1011 gp_Pnt2d p = C2d->Value(param);
1014 uvslf[i].param = param;
1015 uvslf[i].normParam = (param - l) / (f - l);
1024 //=============================================================================
1028 //=============================================================================
1030 ostream & StdMeshers_Quadrangle_2D::SaveTo(ostream & save)
1035 //=============================================================================
1039 //=============================================================================
1041 istream & StdMeshers_Quadrangle_2D::LoadFrom(istream & load)
1046 //=============================================================================
1050 //=============================================================================
1052 ostream & operator <<(ostream & save, StdMeshers_Quadrangle_2D & hyp)
1054 return hyp.SaveTo( save );
1057 //=============================================================================
1061 //=============================================================================
1063 istream & operator >>(istream & load, StdMeshers_Quadrangle_2D & hyp)
1065 return hyp.LoadFrom( load );