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>
44 #include <Geom_Surface.hxx>
45 #include <Geom_Curve.hxx>
46 #include <Geom2d_Curve.hxx>
47 #include <GeomAdaptor_Curve.hxx>
48 #include <GCPnts_UniformAbscissa.hxx>
50 #include <Precision.hxx>
51 #include <gp_Pnt2d.hxx>
52 #include <TColStd_ListIteratorOfListOfInteger.hxx>
54 #include "utilities.h"
55 #include "Utils_ExceptHandlers.hxx"
58 //=============================================================================
62 //=============================================================================
64 StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, SMESH_Gen* gen)
65 : SMESH_2D_Algo(hypId, studyId, gen)
67 MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D");
68 _name = "Quadrangle_2D";
69 // _shapeType = TopAbs_FACE;
70 _shapeType = (1 << TopAbs_FACE);
73 //=============================================================================
77 //=============================================================================
79 StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D()
81 MESSAGE("StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D");
84 //=============================================================================
88 //=============================================================================
90 bool StdMeshers_Quadrangle_2D::CheckHypothesis
92 const TopoDS_Shape& aShape,
93 SMESH_Hypothesis::Hypothesis_Status& aStatus)
95 //MESSAGE("StdMeshers_Quadrangle_2D::CheckHypothesis");
98 aStatus = SMESH_Hypothesis::HYP_OK;
105 //=============================================================================
109 //=============================================================================
111 bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh,
112 const TopoDS_Shape& aShape) throw (SALOME_Exception)
114 Unexpect aCatch(SalomeException);
115 //MESSAGE("StdMeshers_Quadrangle_2D::Compute");
116 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
117 SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape);
119 FaceQuadStruct *quad = CheckAnd2Dcompute(aMesh, aShape);
123 // --- compute 3D values on points, store points & quadrangles
125 int nbdown = quad->nbPts[0];
126 int nbup = quad->nbPts[2];
127 // bool isDownOut = (nbdown > nbup);
128 // bool isUpOut = (nbdown < nbup);
130 int nbright = quad->nbPts[1];
131 int nbleft = quad->nbPts[3];
132 // bool isRightOut = (nbright > nbleft);
133 // bool isLeftOut = (nbright < nbleft);
135 int nbhoriz = Min(nbdown, nbup);
136 int nbvertic = Min(nbright, nbleft);
138 int nbVertices = nbhoriz * nbvertic;
139 int nbQuad = (nbhoriz - 1) * (nbvertic - 1);
140 //SCRUTE(nbVertices);
143 // const TopoDS_Face& FF = TopoDS::Face(aShape);
144 // bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
145 // TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
146 const TopoDS_Face& F = TopoDS::Face(aShape);
147 bool faceIsForward = (F.Orientation() == TopAbs_FORWARD);
148 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
150 // internal mesh nodes
152 for (i = 1; i < nbhoriz - 1; i++) {
153 for (j = 1; j < nbvertic - 1; j++) {
154 int ij = j * nbhoriz + i;
155 double u = quad->uv_grid[ij].u;
156 double v = quad->uv_grid[ij].v;
157 gp_Pnt P = S->Value(u, v);
158 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
159 meshDS->SetNodeOnFace(node, F);
160 quad->uv_grid[ij].node = node;
161 SMDS_FacePosition* fpos =
162 dynamic_cast<SMDS_FacePosition*>(node->GetPosition().get());
163 fpos->SetUParameter(u);
164 fpos->SetVParameter(v);
171 // --.--.--.--.--.-- nbvertic
177 // ---.----.----.--- 0
178 // 0 > > > > > > > > nbhoriz
184 int iup = nbhoriz - 1;
185 if (quad->isEdgeOut[3]) { ilow++; } else { if (quad->isEdgeOut[1]) iup--; }
188 int jup = nbvertic - 1;
189 if (quad->isEdgeOut[0]) { jlow++; } else { if (quad->isEdgeOut[2]) jup--; }
191 // regular quadrangles
192 // bool isQuadForward = ( faceIsForward == quad->isEdgeForward[0]);
193 for (i = ilow; i < iup; i++) {
194 for (j = jlow; j < jup; j++) {
195 const SMDS_MeshNode *a, *b, *c, *d;
196 a = quad->uv_grid[j * nbhoriz + i].node;
197 b = quad->uv_grid[j * nbhoriz + i + 1].node;
198 c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node;
199 d = quad->uv_grid[(j + 1) * nbhoriz + i].node;
200 // if (isQuadForward) faceId = meshDS->AddFace(a,b,c,d);
201 // else faceId = meshDS->AddFace(a,d,c,b);
202 SMDS_MeshFace * face = meshDS->AddFace(a, b, c, d);
203 meshDS->SetMeshElementOnShape(face, F);
207 UVPtStruct *uv_e0 = quad->uv_edges[0];
208 UVPtStruct *uv_e1 = quad->uv_edges[1];
209 UVPtStruct *uv_e2 = quad->uv_edges[2];
210 UVPtStruct *uv_e3 = quad->uv_edges[3];
212 double eps = Precision::Confusion();
214 // Boundary quadrangles
216 if (quad->isEdgeOut[0]) {
219 // |___|___|___|___|___|___|
221 // |___|___|___|___|___|___|
223 // |___|___|___|___|___|___| __ first row of the regular grid
224 // . . . . . . . . . __ down edge nodes
226 // >->->->->->->->->->->->-> -- direction of processing
228 int g = 0; // number of last processed node in the regular grid
230 // number of last node of the down edge to be processed
231 int stop = nbdown - 1;
232 // if right edge is out, we will stop at a node, previous to the last one
233 if (quad->isEdgeOut[1]) stop--;
235 // for each node of the down edge find nearest node
236 // in the first row of the regular grid and link them
237 for (i = 0; i < stop; i++) {
238 const SMDS_MeshNode *a, *b, *c, *d;
240 b = uv_e0[i + 1].node;
241 gp_Pnt pb (b->X(), b->Y(), b->Z());
243 // find node c in the regular grid, which will be linked with node b
246 // right bound reached, link with the rightmost node
248 c = quad->uv_grid[nbhoriz + iup].node;
250 // find in the grid node c, nearest to the b
251 double mind = RealLast();
252 for (int k = g; k <= iup; k++) {
254 const SMDS_MeshNode *nk;
255 if (k < ilow) // this can be, if left edge is out
256 nk = uv_e3[1].node; // get node from the left edge
258 nk = quad->uv_grid[nbhoriz + k].node; // get one of middle nodes
260 gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
261 double dist = pb.Distance(pnk);
262 if (dist < mind - eps) {
272 if (near == g) { // make triangle
273 SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
274 meshDS->SetMeshElementOnShape(face, F);
275 } else { // make quadrangle
279 d = quad->uv_grid[nbhoriz + near - 1].node;
280 SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
281 meshDS->SetMeshElementOnShape(face, F);
283 // if node d is not at position g - make additional triangles
285 for (int k = near - 1; k > g; k--) {
286 c = quad->uv_grid[nbhoriz + k].node;
290 d = quad->uv_grid[nbhoriz + k - 1].node;
291 SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
292 meshDS->SetMeshElementOnShape(face, F);
299 if (quad->isEdgeOut[2]) {
302 // <-<-<-<-<-<-<-<-<-<-<-<-< -- direction of processing
304 // . . . . . . . . . __ up edge nodes
305 // ___ ___ ___ ___ ___ ___ __ first row of the regular grid
307 // |___|___|___|___|___|___|
309 // |___|___|___|___|___|___|
312 int g = nbhoriz - 1; // last processed node in the regular grid
315 // if left edge is out, we will stop at a second node
316 if (quad->isEdgeOut[3]) stop++;
318 // for each node of the up edge find nearest node
319 // in the first row of the regular grid and link them
320 for (i = nbup - 1; i > stop; i--) {
321 const SMDS_MeshNode *a, *b, *c, *d;
323 b = uv_e2[i - 1].node;
324 gp_Pnt pb (b->X(), b->Y(), b->Z());
326 // find node c in the grid, which will be linked with node b
328 if (i == stop + 1) { // left bound reached, link with the leftmost node
329 c = quad->uv_grid[nbhoriz*(nbvertic - 2) + ilow].node;
332 // find node c in the grid, nearest to the b
333 double mind = RealLast();
334 for (int k = g; k >= ilow; k--) {
335 const SMDS_MeshNode *nk;
337 nk = uv_e1[nbright - 2].node;
339 nk = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
340 gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
341 double dist = pb.Distance(pnk);
342 if (dist < mind - eps) {
352 if (near == g) { // make triangle
353 SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
354 meshDS->SetMeshElementOnShape(face, F);
355 } else { // make quadrangle
357 d = uv_e1[nbright - 2].node;
359 d = quad->uv_grid[nbhoriz*(nbvertic - 2) + near + 1].node;
360 SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
361 meshDS->SetMeshElementOnShape(face, F);
363 if (near + 1 < g) { // if d not is at g - make additional triangles
364 for (int k = near + 1; k < g; k++) {
365 c = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
367 d = uv_e1[nbright - 2].node;
369 d = quad->uv_grid[nbhoriz*(nbvertic - 2) + k + 1].node;
370 SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
371 meshDS->SetMeshElementOnShape(face, F);
380 // right or left boundary quadrangles
381 if (quad->isEdgeOut[1]) {
382 // MESSAGE("right edge is out");
383 int g = 0; // last processed node in the grid
384 int stop = nbright - 1;
385 if (quad->isEdgeOut[2]) stop--;
386 for (i = 0; i < stop; i++) {
387 const SMDS_MeshNode *a, *b, *c, *d;
389 b = uv_e1[i + 1].node;
390 gp_Pnt pb (b->X(), b->Y(), b->Z());
392 // find node c in the grid, nearest to the b
394 if (i == stop - 1) { // up bondary reached
395 c = quad->uv_grid[nbhoriz*(jup + 1) - 2].node;
398 double mind = RealLast();
399 for (int k = g; k <= jup; k++) {
400 const SMDS_MeshNode *nk;
402 nk = uv_e0[nbdown - 2].node;
404 nk = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
405 gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
406 double dist = pb.Distance(pnk);
407 if (dist < mind - eps) {
417 if (near == g) { // make triangle
418 SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
419 meshDS->SetMeshElementOnShape(face, F);
420 } else { // make quadrangle
422 d = uv_e0[nbdown - 2].node;
424 d = quad->uv_grid[nbhoriz*near - 2].node;
425 SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
426 meshDS->SetMeshElementOnShape(face, F);
428 if (near - 1 > g) { // if d not is at g - make additional triangles
429 for (int k = near - 1; k > g; k--) {
430 c = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
432 d = uv_e0[nbdown - 2].node;
434 d = quad->uv_grid[nbhoriz*k - 2].node;
435 SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
436 meshDS->SetMeshElementOnShape(face, F);
443 if (quad->isEdgeOut[3]) {
444 // MESSAGE("left edge is out");
445 int g = nbvertic - 1; // last processed node in the grid
447 if (quad->isEdgeOut[0]) stop++;
448 for (i = nbleft - 1; i > stop; i--) {
449 const SMDS_MeshNode *a, *b, *c, *d;
451 b = uv_e3[i - 1].node;
452 gp_Pnt pb (b->X(), b->Y(), b->Z());
454 // find node c in the grid, nearest to the b
456 if (i == stop + 1) { // down bondary reached
457 c = quad->uv_grid[nbhoriz*jlow + 1].node;
460 double mind = RealLast();
461 for (int k = g; k >= jlow; k--) {
462 const SMDS_MeshNode *nk;
466 nk = quad->uv_grid[nbhoriz*k + 1].node;
467 gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
468 double dist = pb.Distance(pnk);
469 if (dist < mind - eps) {
479 if (near == g) { // make triangle
480 SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
481 meshDS->SetMeshElementOnShape(face, F);
482 } else { // make quadrangle
486 d = quad->uv_grid[nbhoriz*(near + 1) + 1].node;
487 SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
488 meshDS->SetMeshElementOnShape(face, F);
490 if (near + 1 < g) { // if d not is at g - make additional triangles
491 for (int k = near + 1; k < g; k++) {
492 c = quad->uv_grid[nbhoriz*k + 1].node;
496 d = quad->uv_grid[nbhoriz*(k + 1) + 1].node;
497 SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
498 meshDS->SetMeshElementOnShape(face, F);
512 //=============================================================================
516 //=============================================================================
518 FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute
519 (SMESH_Mesh & aMesh, const TopoDS_Shape & aShape) throw(SALOME_Exception)
521 Unexpect aCatch(SalomeException);
522 // MESSAGE("StdMeshers_Quadrangle_2D::CheckAnd2Dcompute");
524 SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape);
526 // const TopoDS_Face& FF = TopoDS::Face(aShape);
527 // bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
528 // TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
529 const TopoDS_Face & F = TopoDS::Face(aShape);
530 bool faceIsForward = (F.Orientation() == TopAbs_FORWARD);
532 // verify 1 wire only, with 4 edges
534 if (NumberOfWires(F) != 1)
536 MESSAGE("only 1 wire by face (quadrangles)");
538 //throw SALOME_Exception(LOCALIZED("only 1 wire by face (quadrangles)"));
540 // const TopoDS_Wire WW = BRepTools::OuterWire(F);
541 // TopoDS_Wire W = TopoDS::Wire(WW.Oriented(TopAbs_FORWARD));
542 const TopoDS_Wire& W = BRepTools::OuterWire(F);
543 BRepTools_WireExplorer wexp (W, F);
545 FaceQuadStruct *quad = new FaceQuadStruct;
546 for (int i = 0; i < 4; i++)
547 quad->uv_edges[i] = 0;
551 for (wexp.Init(W, F); wexp.More(); wexp.Next())
553 // const TopoDS_Edge& EE = wexp.Current();
554 // TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
555 const TopoDS_Edge& E = wexp.Current();
556 int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
559 quad->edge[nbEdges] = E;
560 quad->nbPts[nbEdges] = nb + 2; // internal points + 2 extrema
567 MESSAGE("face must have 4 edges /quadrangles");
570 //throw SALOME_Exception(LOCALIZED("face must have 4 edges /quadrangles"));
573 // if (quad->nbPts[0] != quad->nbPts[2]) {
574 // MESSAGE("different point number-opposed edge");
577 // //throw SALOME_Exception(LOCALIZED("different point number-opposed edge"));
580 // if (quad->nbPts[1] != quad->nbPts[3]) {
581 // MESSAGE("different point number-opposed edge");
584 // //throw SALOME_Exception(LOCALIZED("different point number-opposed edge"));
587 // set normalized grid on unit square in parametric domain
589 SetNormalizedGrid(aMesh, F, quad);
594 //=============================================================================
598 //=============================================================================
600 void StdMeshers_Quadrangle_2D::QuadDelete (FaceQuadStruct * quad)
602 //MESSAGE("StdMeshers_Quadrangle_2D::QuadDelete");
605 for (int i = 0; i < 4; i++)
607 if (quad->uv_edges[i])
608 delete [] quad->uv_edges[i];
609 quad->edge[i].Nullify();
612 delete [] quad->uv_grid;
617 //=============================================================================
621 //=============================================================================
623 void StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh,
624 const TopoDS_Shape& aShape,
625 FaceQuadStruct* quad) throw (SALOME_Exception)
627 Unexpect aCatch(SalomeException);
628 // Algorithme décrit dans "Génération automatique de maillages"
629 // P.L. GEORGE, MASSON, § 6.4.1 p. 84-85
630 // traitement dans le domaine paramétrique 2d u,v
631 // transport - projection sur le carré unité
633 // MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid");
634 const TopoDS_Face& F = TopoDS::Face(aShape);
636 // 1 --- find orientation of the 4 edges, by test on extrema
639 // |<----north-2-------^ a3 -------------> a2
641 // west-3 east-1 =right | |
645 // v----south-0--------> a0 -------------> a1
650 Handle(Geom2d_Curve) c2d[4];
653 for (int i = 0; i < 4; i++)
655 c2d[i] = BRep_Tool::CurveOnSurface(quad->edge[i], F,
656 quad->first[i], quad->last[i]);
657 pf[i] = c2d[i]->Value(quad->first[i]);
658 pl[i] = c2d[i]->Value(quad->last[i]);
659 quad->isEdgeForward[i] = false;
662 double eps2d = 1.e-3; // *** utiliser plutot TopExp::CommonVertex, puis
663 // distances si piece fausse
665 if ((pf[1].Distance(pl[0]) < eps2d) || (pl[1].Distance(pl[0]) < eps2d))
667 quad->isEdgeForward[0] = true;
669 double tmp = quad->first[0];
670 quad->first[0] = quad->last[0];
672 pf[0] = c2d[0]->Value(quad->first[0]);
673 pl[0] = c2d[0]->Value(quad->last[0]);
676 for (int i = 1; i < 4; i++)
678 quad->isEdgeForward[i] = (pf[i].Distance(pl[i - 1]) < eps2d);
679 if (!quad->isEdgeForward[i])
681 double tmp = quad->first[i];
682 quad->first[i] = quad->last[i];
684 pf[i] = c2d[i]->Value(quad->first[i]);
685 pl[i] = c2d[i]->Value(quad->last[i]);
686 //SCRUTE(pf[i].Distance(pl[i-1]));
687 ASSERT(pf[i].Distance(pl[i - 1]) < eps2d);
690 //SCRUTE(pf[0].Distance(pl[3]));
691 ASSERT(pf[0].Distance(pl[3]) < eps2d);
693 // for (int i=0; i<4; i++)
695 // SCRUTE(quad->isEdgeForward[i]);
696 // MESSAGE(" -first "<<i<<" "<<pf[i].X()<<" "<<pf[i].Y());
697 // MESSAGE(" -last "<<i<<" "<<pl[i].X()<<" "<<pl[i].Y());
700 // 2 --- load 2d edge points (u,v) with orientation and value on unit square
703 for (int i = 0; i < 2; i++)
705 quad->uv_edges[i] = LoadEdgePoints(aMesh, F, quad->edge[i],
706 quad->first[i], quad->last[i]);
707 if (!quad->uv_edges[i]) loadOk = false;
708 // quad->isEdgeForward[i]);
711 for (int i = 2; i < 4; i++)
713 quad->uv_edges[i] = LoadEdgePoints(aMesh, F, quad->edge[i],
714 quad->last[i], quad->first[i]);
715 if (!quad->uv_edges[i]) loadOk = false;
716 // !quad->isEdgeForward[i]);
721 // MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid - LoadEdgePoints failed");
726 // 3 --- 2D normalized values on unit square [0..1][0..1]
728 // int nbdown = quad->nbPts[0];
729 // int nbright = quad->nbPts[1];
730 int nbhoriz = Min(quad->nbPts[0], quad->nbPts[2]);
731 int nbvertic = Min(quad->nbPts[1], quad->nbPts[3]);
732 // MESSAGE("nbhoriz, nbvertic = " << nbhoriz << nbvertic);
734 quad->isEdgeOut[0] = (quad->nbPts[0] > quad->nbPts[2]);
735 quad->isEdgeOut[1] = (quad->nbPts[1] > quad->nbPts[3]);
736 quad->isEdgeOut[2] = (quad->nbPts[2] > quad->nbPts[0]);
737 quad->isEdgeOut[3] = (quad->nbPts[3] > quad->nbPts[1]);
739 quad->uv_grid = new UVPtStruct[nbvertic * nbhoriz];
741 UVPtStruct *uv_grid = quad->uv_grid;
742 UVPtStruct *uv_e0 = quad->uv_edges[0];
743 UVPtStruct *uv_e1 = quad->uv_edges[1];
744 UVPtStruct *uv_e2 = quad->uv_edges[2];
745 UVPtStruct *uv_e3 = quad->uv_edges[3];
747 // nodes Id on "in" edges
748 if (! quad->isEdgeOut[0]) {
750 for (int i = 0; i < nbhoriz; i++) { // down
751 int ij = j * nbhoriz + i;
752 uv_grid[ij].node = uv_e0[i].node;
755 if (! quad->isEdgeOut[1]) {
757 for (int j = 0; j < nbvertic; j++) { // right
758 int ij = j * nbhoriz + i;
759 uv_grid[ij].node = uv_e1[j].node;
762 if (! quad->isEdgeOut[2]) {
763 int j = nbvertic - 1;
764 for (int i = 0; i < nbhoriz; i++) { // up
765 int ij = j * nbhoriz + i;
766 uv_grid[ij].node = uv_e2[i].node;
769 if (! quad->isEdgeOut[3]) {
771 for (int j = 0; j < nbvertic; j++) { // left
772 int ij = j * nbhoriz + i;
773 uv_grid[ij].node = uv_e3[j].node;
777 // falsificate "out" edges
778 if (quad->isEdgeOut[0]) // down
779 uv_e0 = MakeEdgePoints
780 (aMesh, F, quad->edge[0], quad->first[0], quad->last[0], nbhoriz - 1);
781 else if (quad->isEdgeOut[2]) // up
782 uv_e2 = MakeEdgePoints
783 (aMesh, F, quad->edge[2], quad->last[2], quad->first[2], nbhoriz - 1);
785 if (quad->isEdgeOut[1]) // right
786 uv_e1 = MakeEdgePoints
787 (aMesh, F, quad->edge[1], quad->first[1], quad->last[1], nbvertic - 1);
788 else if (quad->isEdgeOut[3]) // left
789 uv_e3 = MakeEdgePoints
790 (aMesh, F, quad->edge[3], quad->last[3], quad->first[3], nbvertic - 1);
792 // normalized 2d values on grid
793 for (int i = 0; i < nbhoriz; i++)
795 for (int j = 0; j < nbvertic; j++)
797 int ij = j * nbhoriz + i;
798 // --- droite i cste : x = x0 + y(x1-x0)
799 double x0 = uv_e0[i].normParam; // bas - sud
800 double x1 = uv_e2[i].normParam; // haut - nord
801 // --- droite j cste : y = y0 + x(y1-y0)
802 double y0 = uv_e3[j].normParam; // gauche-ouest
803 double y1 = uv_e1[j].normParam; // droite - est
804 // --- intersection : x=x0+(y0+x(y1-y0))(x1-x0)
805 double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
806 double y = y0 + x * (y1 - y0);
809 //MESSAGE("-xy-01 "<<x0<<" "<<x1<<" "<<y0<<" "<<y1);
810 //MESSAGE("-xy-norm "<<i<<" "<<j<<" "<<x<<" "<<y);
814 // 4 --- projection on 2d domain (u,v)
820 for (int i = 0; i < nbhoriz; i++)
822 for (int j = 0; j < nbvertic; j++)
824 int ij = j * nbhoriz + i;
825 double x = uv_grid[ij].x;
826 double y = uv_grid[ij].y;
827 double param_0 = uv_e0[0].param + x * (uv_e0[nbhoriz - 1].param - uv_e0[0].param); // sud
828 double param_2 = uv_e2[0].param + x * (uv_e2[nbhoriz - 1].param - uv_e2[0].param); // nord
829 double param_1 = uv_e1[0].param + y * (uv_e1[nbvertic - 1].param - uv_e1[0].param); // est
830 double param_3 = uv_e3[0].param + y * (uv_e3[nbvertic - 1].param - uv_e3[0].param); // ouest
832 //MESSAGE("params "<<param_0<<" "<<param_1<<" "<<param_2<<" "<<param_3);
833 gp_Pnt2d p0 = c2d[0]->Value(param_0);
834 gp_Pnt2d p1 = c2d[1]->Value(param_1);
835 gp_Pnt2d p2 = c2d[2]->Value(param_2);
836 gp_Pnt2d p3 = c2d[3]->Value(param_3);
838 double u = (1 - y) * p0.X() + x * p1.X() + y * p2.X() + (1 - x) * p3.X();
839 double v = (1 - y) * p0.Y() + x * p1.Y() + y * p2.Y() + (1 - x) * p3.Y();
841 u -= (1 - x) * (1 - y) * a0.X() + x * (1 - y) * a1.X() +
842 x * y * a2.X() + (1 - x) * y * a3.X();
843 v -= (1 - x) * (1 - y) * a0.Y() + x * (1 - y) * a1.Y() +
844 x * y * a2.Y() + (1 - x) * y * a3.Y();
849 //MESSAGE("-uv- "<<i<<" "<<j<<" "<<uv_grid[ij].u<<" "<<uv_grid[ij].v);
854 //=============================================================================
858 //=============================================================================
859 UVPtStruct* StdMeshers_Quadrangle_2D::LoadEdgePoints (SMESH_Mesh & aMesh,
860 const TopoDS_Face& F,
861 const TopoDS_Edge& E,
862 double first, double last)
865 //MESSAGE("StdMeshers_Quadrangle_2D::LoadEdgePoints");
867 SMDS_Mesh* meshDS = aMesh.GetMeshDS();
869 // --- IDNodes of first and last Vertex
871 TopoDS_Vertex VFirst, VLast;
872 TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
874 ASSERT(!VFirst.IsNull());
875 SMDS_NodeIteratorPtr lid = aMesh.GetSubMesh(VFirst)->GetSubMeshDS()->GetNodes();
878 MESSAGE ( "NO NODE BUILT ON VERTEX" );
881 const SMDS_MeshNode* idFirst = lid->next();
883 ASSERT(!VLast.IsNull());
884 lid = aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes();
887 MESSAGE ( "NO NODE BUILT ON VERTEX" );
890 const SMDS_MeshNode* idLast = lid->next();
892 // --- edge internal IDNodes (relies on good order storage, not checked)
894 map<double, const SMDS_MeshNode *> params;
895 SMDS_NodeIteratorPtr ite = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes();
899 const SMDS_MeshNode* node = ite->next();
900 const SMDS_EdgePosition* epos =
901 static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
902 double param = epos->GetUParameter();
903 params[param] = node;
906 int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
907 if (nbPoints != params.size())
909 MESSAGE( "BAD NODE ON EDGE POSITIONS" );
912 UVPtStruct* uvslf = new UVPtStruct[nbPoints + 2];
915 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
917 bool isForward = (((l - f) * (last - first)) > 0);
924 gp_Pnt2d p = C2d->Value(f); // first point = Vertex Forward
928 uvslf[0].node = idFirst;
929 //MESSAGE("__ f "<<f<<" "<<uvslf[0].x <<" "<<uvslf[0].y);
930 map < double, const SMDS_MeshNode* >::iterator itp = params.begin();
931 for (int i = 1; i <= nbPoints; i++) // nbPoints internal
933 double param = (*itp).first;
934 gp_Pnt2d p = C2d->Value(param);
937 uvslf[i].param = param;
938 uvslf[i].node = (*itp).second;
939 //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[i].x <<" "<<uvslf[i].y);
942 p = C2d->Value(l); // last point = Vertex Reversed
943 uvslf[nbPoints + 1].x = p.X();
944 uvslf[nbPoints + 1].y = p.Y();
945 uvslf[nbPoints + 1].param = l;
946 uvslf[nbPoints + 1].node = idLast;
947 //MESSAGE("__ l "<<l<<" "<<uvslf[nbPoints+1].x <<" "<<uvslf[nbPoints+1].y);
952 gp_Pnt2d p = C2d->Value(l); // first point = Vertex Reversed
956 uvslf[0].node = idLast;
957 //MESSAGE("__ l "<<l<<" "<<uvslf[0].x <<" "<<uvslf[0].y);
958 map < double, const SMDS_MeshNode* >::reverse_iterator itp = params.rbegin();
960 for (int j = nbPoints; j >= 1; j--) // nbPoints internal
962 double param = (*itp).first;
963 int i = nbPoints + 1 - j;
964 gp_Pnt2d p = C2d->Value(param);
967 uvslf[i].param = param;
968 uvslf[i].node = (*itp).second;
969 //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[i].x <<" "<<uvslf[i].y);
972 p = C2d->Value(f); // last point = Vertex Forward
973 uvslf[nbPoints + 1].x = p.X();
974 uvslf[nbPoints + 1].y = p.Y();
975 uvslf[nbPoints + 1].param = f;
976 uvslf[nbPoints + 1].node = idFirst;
977 //MESSAGE("__ f "<<f<<" "<<uvslf[nbPoints+1].x <<" "<<uvslf[nbPoints+1].y);
980 ASSERT(paramin != paramax);
981 for (int i = 0; i < nbPoints + 2; i++)
983 uvslf[i].normParam = (uvslf[i].param - paramin) / (paramax - paramin);
984 //SCRUTE(uvslf[i].normParam);
990 //=============================================================================
994 //=============================================================================
995 UVPtStruct* StdMeshers_Quadrangle_2D::MakeEdgePoints (SMESH_Mesh & aMesh,
996 const TopoDS_Face& F,
997 const TopoDS_Edge& E,
998 double first, double last,
1001 // MESSAGE("StdMeshers_Quadrangle_2D::MakeEdgePoints");
1003 UVPtStruct* uvslf = new UVPtStruct[nb_segm + 1];
1004 list<double> params;
1006 // --- edge internal points
1008 Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, fi, li);
1009 if (!Curve.IsNull()) {
1011 GeomAdaptor_Curve C3d (Curve);
1012 double length = EdgeLength(E);
1013 double eltSize = length / nb_segm;
1014 GCPnts_UniformAbscissa Discret (C3d, eltSize, fi, li);
1015 if (!Discret.IsDone()) return false;
1016 int NbPoints = Discret.NbPoints();
1017 for (int i = 1; i <= NbPoints; i++) {
1018 double param = Discret.Parameter(i);
1019 params.push_back(param);
1022 catch (Standard_Failure) {
1028 // Edge is a degenerated Edge
1029 BRep_Tool::Range(E, fi, li);
1030 double du = (li - fi) / nb_segm;
1031 for (int i = 1; i <= nb_segm + 1; i++)
1033 double param = fi + (i - 1) * du;
1034 params.push_back(param);
1039 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
1042 bool isForward = (((l - f) * (last - first)) > 0);
1044 list<double>::iterator itU = params.begin();
1045 for (int i = 0; i <= nb_segm; i++) // nbPoints internal
1047 double param = *itU;
1048 gp_Pnt2d p = C2d->Value(param);
1051 uvslf[i].param = param;
1052 uvslf[i].normParam = (param - f) / (l - f);
1056 list<double>::reverse_iterator itU = params.rbegin();
1057 for (int j = nb_segm; j >= 0; j--) // nbPoints internal
1059 double param = *itU;
1060 int i = nb_segm - j;
1061 gp_Pnt2d p = C2d->Value(param);
1064 uvslf[i].param = param;
1065 uvslf[i].normParam = (param - l) / (f - l);
1074 //=============================================================================
1078 //=============================================================================
1080 ostream & StdMeshers_Quadrangle_2D::SaveTo(ostream & save)
1085 //=============================================================================
1089 //=============================================================================
1091 istream & StdMeshers_Quadrangle_2D::LoadFrom(istream & load)
1096 //=============================================================================
1100 //=============================================================================
1102 ostream & operator <<(ostream & save, StdMeshers_Quadrangle_2D & hyp)
1104 return hyp.SaveTo( save );
1107 //=============================================================================
1111 //=============================================================================
1113 istream & operator >>(istream & load, StdMeshers_Quadrangle_2D & hyp)
1115 return hyp.LoadFrom( load );