Salome HOME
da6b6a2e4ebff1db3d8a6f701705963ff33e7329
[modules/smesh.git] / src / StdMeshers / StdMeshers_Quadrangle_2D.cxx
1 //  SMESH SMESH : implementaion of SMESH idl descriptions
2 //
3 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS 
5 // 
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. 
10 // 
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. 
15 // 
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 
19 // 
20 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org 
21 //
22 //
23 //
24 //  File   : StdMeshers_Quadrangle_2D.cxx
25 //           Moved here from SMESH_Quadrangle_2D.cxx
26 //  Author : Paul RASCLE, EDF
27 //  Module : SMESH
28 //  $Header$
29
30 using namespace std;
31 #include "StdMeshers_Quadrangle_2D.hxx"
32 #include "SMESH_Gen.hxx"
33 #include "SMESH_Mesh.hxx"
34 #include "SMESH_subMesh.hxx"
35
36 #include "SMDS_MeshElement.hxx"
37 #include "SMDS_MeshNode.hxx"
38 #include "SMDS_EdgePosition.hxx"
39 #include "SMDS_FacePosition.hxx"
40
41 #include <BRep_Tool.hxx>
42 #include <BRepTools.hxx>
43 #include <BRepTools_WireExplorer.hxx>
44
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>
50
51 #include <Precision.hxx>
52 #include <gp_Pnt2d.hxx>
53 #include <TColStd_ListIteratorOfListOfInteger.hxx>
54
55 #include "utilities.h"
56 #include "Utils_ExceptHandlers.hxx"
57
58
59 //=============================================================================
60 /*!
61  *  
62  */
63 //=============================================================================
64
65 StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, SMESH_Gen* gen)
66      : SMESH_2D_Algo(hypId, studyId, gen)
67 {
68   MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D");
69   _name = "Quadrangle_2D";
70   //  _shapeType = TopAbs_FACE;
71   _shapeType = (1 << TopAbs_FACE);
72 }
73
74 //=============================================================================
75 /*!
76  *  
77  */
78 //=============================================================================
79
80 StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D()
81 {
82   MESSAGE("StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D");
83 }
84
85 //=============================================================================
86 /*!
87  *  
88  */
89 //=============================================================================
90
91 bool StdMeshers_Quadrangle_2D::CheckHypothesis
92                          (SMESH_Mesh& aMesh,
93                           const TopoDS_Shape& aShape,
94                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
95 {
96   //MESSAGE("StdMeshers_Quadrangle_2D::CheckHypothesis");
97
98   bool isOk = true;
99   aStatus = SMESH_Hypothesis::HYP_OK;
100
101   // nothing to check
102
103   return isOk;
104 }
105
106 //=============================================================================
107 /*!
108  *  
109  */
110 //=============================================================================
111
112 bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh,
113                                         const TopoDS_Shape& aShape) throw (SALOME_Exception)
114 {
115   Unexpect aCatch(SalomeException);
116   //MESSAGE("StdMeshers_Quadrangle_2D::Compute");
117   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
118   aMesh.GetSubMesh(aShape);
119
120   FaceQuadStruct *quad = CheckAnd2Dcompute(aMesh, aShape);
121   if (!quad)
122     return false;
123
124   // --- compute 3D values on points, store points & quadrangles
125
126   int nbdown  = quad->nbPts[0];
127   int nbup    = quad->nbPts[2];
128 //  bool isDownOut = (nbdown > nbup);
129 //  bool isUpOut   = (nbdown < nbup);
130
131   int nbright = quad->nbPts[1];
132   int nbleft  = quad->nbPts[3];
133 //  bool isRightOut = (nbright > nbleft);
134 //  bool isLeftOut  = (nbright < nbleft);
135
136   int nbhoriz  = Min(nbdown, nbup);
137   int nbvertic = Min(nbright, nbleft);
138
139   int nbVertices = nbhoriz * nbvertic;
140   int nbQuad = (nbhoriz - 1) * (nbvertic - 1);
141   //SCRUTE(nbVertices);
142   //SCRUTE(nbQuad);
143
144   //   const TopoDS_Face& FF = TopoDS::Face(aShape);
145   //   bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
146   //   TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
147   const TopoDS_Face& F = TopoDS::Face(aShape);
148   bool faceIsForward = (F.Orientation() == TopAbs_FORWARD);
149   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
150
151   // internal mesh nodes
152   int i, j;
153   for (i = 1; i < nbhoriz - 1; i++) {
154     for (j = 1; j < nbvertic - 1; j++) {
155       int ij = j * nbhoriz + i;
156       double u = quad->uv_grid[ij].u;
157       double v = quad->uv_grid[ij].v;
158       gp_Pnt P = S->Value(u, v);
159       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
160       meshDS->SetNodeOnFace(node, F);
161       quad->uv_grid[ij].node = node;
162       SMDS_FacePosition* fpos =
163         dynamic_cast<SMDS_FacePosition*>(node->GetPosition().get());
164       fpos->SetUParameter(u);
165       fpos->SetVParameter(v);
166     }
167   }
168
169   // mesh faces
170
171   //             [2]
172   //      --.--.--.--.--.--  nbvertic
173   //     |                 | ^
174   //     |                 | ^
175   // [3] |                 | ^ j  [1]
176   //     |                 | ^
177   //     |                 | ^
178   //      ---.----.----.---  0
179   //     0 > > > > > > > > nbhoriz
180   //              i
181   //             [0]
182
183   i = 0;
184   int ilow = 0;
185   int iup = nbhoriz - 1;
186   if (quad->isEdgeOut[3]) { ilow++; } else { if (quad->isEdgeOut[1]) iup--; }
187
188   int jlow = 0;
189   int jup = nbvertic - 1;
190   if (quad->isEdgeOut[0]) { jlow++; } else { if (quad->isEdgeOut[2]) jup--; }
191
192   // regular quadrangles
193   //   bool isQuadForward = ( faceIsForward == quad->isEdgeForward[0]);
194   for (i = ilow; i < iup; i++) {
195     for (j = jlow; j < jup; j++) {
196       const SMDS_MeshNode *a, *b, *c, *d;
197       a = quad->uv_grid[j * nbhoriz + i].node;
198       b = quad->uv_grid[j * nbhoriz + i + 1].node;
199       c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node;
200       d = quad->uv_grid[(j + 1) * nbhoriz + i].node;
201       //  if (isQuadForward) faceId = meshDS->AddFace(a,b,c,d);
202       //  else faceId = meshDS->AddFace(a,d,c,b);
203       SMDS_MeshFace * face = meshDS->AddFace(a, b, c, d);
204       meshDS->SetMeshElementOnShape(face, F);
205     }
206   }
207
208   UVPtStruct *uv_e0 = quad->uv_edges[0];
209   UVPtStruct *uv_e1 = quad->uv_edges[1];
210   UVPtStruct *uv_e2 = quad->uv_edges[2];
211   UVPtStruct *uv_e3 = quad->uv_edges[3];
212
213   double eps = Precision::Confusion();
214
215   // Boundary quadrangles
216
217   if (quad->isEdgeOut[0]) {
218     // Down edge is out
219     // 
220     // |___|___|___|___|___|___|
221     // |   |   |   |   |   |   |
222     // |___|___|___|___|___|___|
223     // |   |   |   |   |   |   |
224     // |___|___|___|___|___|___| __ first row of the regular grid
225     // .  .  .  .  .  .  .  .  . __ down edge nodes
226     // 
227     // >->->->->->->->->->->->-> -- direction of processing
228
229     int g = 0; // number of last processed node in the regular grid
230
231     // number of last node of the down edge to be processed
232     int stop = nbdown - 1;
233     // if right edge is out, we will stop at a node, previous to the last one
234     if (quad->isEdgeOut[1]) stop--;
235
236     // for each node of the down edge find nearest node
237     // in the first row of the regular grid and link them
238     for (i = 0; i < stop; i++) {
239       const SMDS_MeshNode *a, *b, *c, *d;
240       a = uv_e0[i].node;
241       b = uv_e0[i + 1].node;
242       gp_Pnt pb (b->X(), b->Y(), b->Z());
243
244       // find node c in the regular grid, which will be linked with node b
245       int near = g;
246       if (i == stop - 1) {
247         // right bound reached, link with the rightmost node
248         near = iup;
249         c = quad->uv_grid[nbhoriz + iup].node;
250       } else {
251         // find in the grid node c, nearest to the b
252         double mind = RealLast();
253         for (int k = g; k <= iup; k++) {
254
255           const SMDS_MeshNode *nk;
256           if (k < ilow) // this can be, if left edge is out
257             nk = uv_e3[1].node; // get node from the left edge
258           else
259             nk = quad->uv_grid[nbhoriz + k].node; // get one of middle nodes
260
261           gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
262           double dist = pb.Distance(pnk);
263           if (dist < mind - eps) {
264             c = nk;
265             near = k;
266             mind = dist;
267           } else {
268             break;
269           }
270         }
271       }
272
273       if (near == g) { // make triangle
274         SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
275         meshDS->SetMeshElementOnShape(face, F);
276       } else { // make quadrangle
277         if (near - 1 < ilow)
278           d = uv_e3[1].node;
279         else
280           d = quad->uv_grid[nbhoriz + near - 1].node;
281         SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
282         meshDS->SetMeshElementOnShape(face, F);
283
284         // if node d is not at position g - make additional triangles
285         if (near - 1 > g) {
286           for (int k = near - 1; k > g; k--) {
287             c = quad->uv_grid[nbhoriz + k].node;
288             if (k - 1 < ilow)
289               d = uv_e3[1].node;
290             else
291               d = quad->uv_grid[nbhoriz + k - 1].node;
292             SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
293             meshDS->SetMeshElementOnShape(face, F);
294           }
295         }
296         g = near;
297       }
298     }
299   } else {
300     if (quad->isEdgeOut[2]) {
301       // Up edge is out
302       // 
303       // <-<-<-<-<-<-<-<-<-<-<-<-< -- direction of processing
304       // 
305       // .  .  .  .  .  .  .  .  . __ up edge nodes
306       //  ___ ___ ___ ___ ___ ___  __ first row of the regular grid
307       // |   |   |   |   |   |   |
308       // |___|___|___|___|___|___|
309       // |   |   |   |   |   |   |
310       // |___|___|___|___|___|___|
311       // |   |   |   |   |   |   |
312
313       int g = nbhoriz - 1; // last processed node in the regular grid
314
315       int stop = 0;
316       // if left edge is out, we will stop at a second node
317       if (quad->isEdgeOut[3]) stop++;
318
319       // for each node of the up edge find nearest node
320       // in the first row of the regular grid and link them
321       for (i = nbup - 1; i > stop; i--) {
322         const SMDS_MeshNode *a, *b, *c, *d;
323         a = uv_e2[i].node;
324         b = uv_e2[i - 1].node;
325         gp_Pnt pb (b->X(), b->Y(), b->Z());
326
327         // find node c in the grid, which will be linked with node b
328         int near = g;
329         if (i == stop + 1) { // left bound reached, link with the leftmost node
330           c = quad->uv_grid[nbhoriz*(nbvertic - 2) + ilow].node;
331           near = ilow;
332         } else {
333           // find node c in the grid, nearest to the b
334           double mind = RealLast();
335           for (int k = g; k >= ilow; k--) {
336             const SMDS_MeshNode *nk;
337             if (k > iup)
338               nk = uv_e1[nbright - 2].node;
339             else
340               nk = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
341             gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
342             double dist = pb.Distance(pnk);
343             if (dist < mind - eps) {
344               c = nk;
345               near = k;
346               mind = dist;
347             } else {
348               break;
349             }
350           }
351         }
352
353         if (near == g) { // make triangle
354           SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
355           meshDS->SetMeshElementOnShape(face, F);
356         } else { // make quadrangle
357           if (near + 1 > iup)
358             d = uv_e1[nbright - 2].node;
359           else
360             d = quad->uv_grid[nbhoriz*(nbvertic - 2) + near + 1].node;
361           SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
362           meshDS->SetMeshElementOnShape(face, F);
363
364           if (near + 1 < g) { // if d not is at g - make additional triangles
365             for (int k = near + 1; k < g; k++) {
366               c = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
367               if (k + 1 > iup)
368                 d = uv_e1[nbright - 2].node;
369               else
370                 d = quad->uv_grid[nbhoriz*(nbvertic - 2) + k + 1].node;
371               SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
372               meshDS->SetMeshElementOnShape(face, F);
373             }
374           }
375           g = near;
376         }
377       }
378     }
379   }
380
381   // right or left boundary quadrangles
382   if (quad->isEdgeOut[1]) {
383 //    MESSAGE("right edge is out");
384     int g = 0; // last processed node in the grid
385     int stop = nbright - 1;
386     if (quad->isEdgeOut[2]) stop--;
387     for (i = 0; i < stop; i++) {
388       const SMDS_MeshNode *a, *b, *c, *d;
389       a = uv_e1[i].node;
390       b = uv_e1[i + 1].node;
391       gp_Pnt pb (b->X(), b->Y(), b->Z());
392
393       // find node c in the grid, nearest to the b
394       int near = g;
395       if (i == stop - 1) { // up bondary reached
396         c = quad->uv_grid[nbhoriz*(jup + 1) - 2].node;
397         near = jup;
398       } else {
399         double mind = RealLast();
400         for (int k = g; k <= jup; k++) {
401           const SMDS_MeshNode *nk;
402           if (k < jlow)
403             nk = uv_e0[nbdown - 2].node;
404           else
405             nk = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
406           gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
407           double dist = pb.Distance(pnk);
408           if (dist < mind - eps) {
409             c = nk;
410             near = k;
411             mind = dist;
412           } else {
413             break;
414           }
415         }
416       }
417
418       if (near == g) { // make triangle
419         SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
420         meshDS->SetMeshElementOnShape(face, F);
421       } else { // make quadrangle
422         if (near - 1 < jlow)
423           d = uv_e0[nbdown - 2].node;
424         else
425           d = quad->uv_grid[nbhoriz*near - 2].node;
426         SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
427         meshDS->SetMeshElementOnShape(face, F);
428
429         if (near - 1 > g) { // if d not is at g - make additional triangles
430           for (int k = near - 1; k > g; k--) {
431             c = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
432             if (k - 1 < jlow)
433               d = uv_e0[nbdown - 2].node;
434             else
435               d = quad->uv_grid[nbhoriz*k - 2].node;
436             SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
437             meshDS->SetMeshElementOnShape(face, F);
438           }
439         }
440         g = near;
441       }
442     }
443   } else {
444     if (quad->isEdgeOut[3]) {
445 //      MESSAGE("left edge is out");
446       int g = nbvertic - 1; // last processed node in the grid
447       int stop = 0;
448       if (quad->isEdgeOut[0]) stop++;
449       for (i = nbleft - 1; i > stop; i--) {
450         const SMDS_MeshNode *a, *b, *c, *d;
451         a = uv_e3[i].node;
452         b = uv_e3[i - 1].node;
453         gp_Pnt pb (b->X(), b->Y(), b->Z());
454
455         // find node c in the grid, nearest to the b
456         int near = g;
457         if (i == stop + 1) { // down bondary reached
458           c = quad->uv_grid[nbhoriz*jlow + 1].node;
459           near = jlow;
460         } else {
461           double mind = RealLast();
462           for (int k = g; k >= jlow; k--) {
463             const SMDS_MeshNode *nk;
464             if (k > jup)
465               nk = uv_e2[1].node;
466             else
467               nk = quad->uv_grid[nbhoriz*k + 1].node;
468             gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
469             double dist = pb.Distance(pnk);
470             if (dist < mind - eps) {
471               c = nk;
472               near = k;
473               mind = dist;
474             } else {
475               break;
476             }
477           }
478         }
479
480         if (near == g) { // make triangle
481           SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
482           meshDS->SetMeshElementOnShape(face, F);
483         } else { // make quadrangle
484           if (near + 1 > jup)
485             d = uv_e2[1].node;
486           else
487             d = quad->uv_grid[nbhoriz*(near + 1) + 1].node;
488           SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
489           meshDS->SetMeshElementOnShape(face, F);
490
491           if (near + 1 < g) { // if d not is at g - make additional triangles
492             for (int k = near + 1; k < g; k++) {
493               c = quad->uv_grid[nbhoriz*k + 1].node;
494               if (k + 1 > jup)
495                 d = uv_e2[1].node;
496               else
497                 d = quad->uv_grid[nbhoriz*(k + 1) + 1].node;
498               SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
499               meshDS->SetMeshElementOnShape(face, F);
500             }
501           }
502           g = near;
503         }
504       }
505     }
506   }
507
508   QuadDelete(quad);
509   bool isOk = true;
510   return isOk;
511 }
512
513 //=============================================================================
514 /*!
515  *  
516  */
517 //=============================================================================
518
519 FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute
520   (SMESH_Mesh & aMesh, const TopoDS_Shape & aShape) throw(SALOME_Exception)
521 {
522   Unexpect aCatch(SalomeException);
523 //  MESSAGE("StdMeshers_Quadrangle_2D::CheckAnd2Dcompute");
524
525   SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape);
526
527   //   const TopoDS_Face& FF = TopoDS::Face(aShape);
528   //   bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
529   //   TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
530   const TopoDS_Face & F = TopoDS::Face(aShape);
531   bool faceIsForward = (F.Orientation() == TopAbs_FORWARD);
532
533   // verify 1 wire only, with 4 edges
534
535   if (NumberOfWires(F) != 1)
536   {
537     MESSAGE("only 1 wire by face (quadrangles)");
538     return 0;
539     //throw SALOME_Exception(LOCALIZED("only 1 wire by face (quadrangles)"));
540   }
541   //   const TopoDS_Wire WW = BRepTools::OuterWire(F);
542   //   TopoDS_Wire W = TopoDS::Wire(WW.Oriented(TopAbs_FORWARD));
543   const TopoDS_Wire& W = BRepTools::OuterWire(F);
544   BRepTools_WireExplorer wexp (W, F);
545
546   FaceQuadStruct *quad = new FaceQuadStruct;
547   for (int i = 0; i < 4; i++)
548     quad->uv_edges[i] = 0;
549   quad->uv_grid = 0;
550
551   int nbEdges = 0;
552   for (wexp.Init(W, F); wexp.More(); wexp.Next())
553   {
554     //       const TopoDS_Edge& EE = wexp.Current();
555     //       TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
556     const TopoDS_Edge& E = wexp.Current();
557     int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
558     if (nbEdges < 4)
559     {
560       quad->edge[nbEdges] = E;
561       quad->nbPts[nbEdges] = nb + 2; // internal points + 2 extrema
562     }
563     nbEdges++;
564   }
565
566   if (nbEdges != 4)
567   {
568     MESSAGE("face must have 4 edges /quadrangles");
569     QuadDelete(quad);
570     return 0;
571     //throw SALOME_Exception(LOCALIZED("face must have 4 edges /quadrangles"));
572   }
573
574 //  if (quad->nbPts[0] != quad->nbPts[2]) {
575 //    MESSAGE("different point number-opposed edge");
576 //    QuadDelete(quad);
577 //    return 0;
578 //    //throw SALOME_Exception(LOCALIZED("different point number-opposed edge"));
579 //  }
580 //
581 //  if (quad->nbPts[1] != quad->nbPts[3]) {
582 //    MESSAGE("different point number-opposed edge");
583 //    QuadDelete(quad);
584 //    return 0;
585 //    //throw SALOME_Exception(LOCALIZED("different point number-opposed edge"));
586 //  }
587
588   // set normalized grid on unit square in parametric domain
589
590   SetNormalizedGrid(aMesh, F, quad);
591
592   return quad;
593 }
594
595 //=============================================================================
596 /*!
597  *  
598  */
599 //=============================================================================
600
601 void StdMeshers_Quadrangle_2D::QuadDelete (FaceQuadStruct * quad)
602 {
603   //MESSAGE("StdMeshers_Quadrangle_2D::QuadDelete");
604   if (quad)
605   {
606     for (int i = 0; i < 4; i++)
607     {
608       if (quad->uv_edges[i])
609         delete [] quad->uv_edges[i];
610       quad->edge[i].Nullify();
611     }
612     if (quad->uv_grid)
613       delete [] quad->uv_grid;
614     delete quad;
615   }
616 }
617
618 //=============================================================================
619 /*!
620  *  
621  */
622 //=============================================================================
623
624 void StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh,
625                                                   const TopoDS_Shape& aShape,
626                                                   FaceQuadStruct* quad) throw (SALOME_Exception)
627 {
628   Unexpect aCatch(SalomeException);
629   // Algorithme décrit dans "Génération automatique de maillages"
630   // P.L. GEORGE, MASSON, § 6.4.1 p. 84-85
631   // traitement dans le domaine paramétrique 2d u,v
632   // transport - projection sur le carré unité
633
634 //  MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid");
635   const TopoDS_Face& F = TopoDS::Face(aShape);
636
637   // 1 --- find orientation of the 4 edges, by test on extrema
638
639   //      max             min                    0     x1     1
640   //     |<----north-2-------^                a3 -------------> a2
641   //     |                   |                   ^1          1^
642   //    west-3            east-1 =right          |            |
643   //     |                   |         ==>       |            |
644   //  y0 |                   | y1                |            |
645   //     |                   |                   |0          0|
646   //     v----south-0-------->                a0 -------------> a1
647   //      min             max                    0     x0     1
648   //             =down
649   //
650
651   Handle(Geom2d_Curve) c2d[4];
652   gp_Pnt2d pf[4];
653   gp_Pnt2d pl[4];
654   for (int i = 0; i < 4; i++)
655   {
656     c2d[i] = BRep_Tool::CurveOnSurface(quad->edge[i], F,
657                                        quad->first[i], quad->last[i]);
658     pf[i] = c2d[i]->Value(quad->first[i]);
659     pl[i] = c2d[i]->Value(quad->last[i]);
660     quad->isEdgeForward[i] = false;
661   }
662
663   double eps2d = 1.e-3; // *** utiliser plutot TopExp::CommonVertex, puis
664   // distances si piece fausse
665 //  int i = 0;
666   if ((pf[1].Distance(pl[0]) < eps2d) || (pl[1].Distance(pl[0]) < eps2d))
667   {
668     quad->isEdgeForward[0] = true;
669   } else {
670     double tmp = quad->first[0];
671     quad->first[0] = quad->last[0];
672     quad->last[0] = tmp;
673     pf[0] = c2d[0]->Value(quad->first[0]);
674     pl[0] = c2d[0]->Value(quad->last[0]);
675   }
676
677   for (int i = 1; i < 4; i++)
678   {
679     quad->isEdgeForward[i] = (pf[i].Distance(pl[i - 1]) < eps2d);
680     if (!quad->isEdgeForward[i])
681     {
682       double tmp = quad->first[i];
683       quad->first[i] = quad->last[i];
684       quad->last[i] = tmp;
685       pf[i] = c2d[i]->Value(quad->first[i]);
686       pl[i] = c2d[i]->Value(quad->last[i]);
687       //SCRUTE(pf[i].Distance(pl[i-1]));
688       ASSERT(pf[i].Distance(pl[i - 1]) < eps2d);
689     }
690   }
691   //SCRUTE(pf[0].Distance(pl[3]));
692   ASSERT(pf[0].Distance(pl[3]) < eps2d);
693
694 //   for (int i=0; i<4; i++)
695 //     {
696 //       SCRUTE(quad->isEdgeForward[i]);
697 //       MESSAGE(" -first "<<i<<" "<<pf[i].X()<<" "<<pf[i].Y());
698 //       MESSAGE(" -last  "<<i<<" "<<pl[i].X()<<" "<<pl[i].Y());
699 //     }
700
701   // 2 --- load 2d edge points (u,v) with orientation and value on unit square
702
703   bool loadOk = true;
704   for (int i = 0; i < 2; i++)
705   {
706     quad->uv_edges[i] = LoadEdgePoints(aMesh, F, quad->edge[i],
707                                        quad->first[i], quad->last[i]);
708     if (!quad->uv_edges[i]) loadOk = false;
709     //    quad->isEdgeForward[i]);
710   }
711
712   for (int i = 2; i < 4; i++)
713   {
714     quad->uv_edges[i] = LoadEdgePoints(aMesh, F, quad->edge[i],
715                                        quad->last[i], quad->first[i]);
716     if (!quad->uv_edges[i]) loadOk = false;
717     //    !quad->isEdgeForward[i]);
718   }
719
720   if (!loadOk)
721   {
722 //    MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid - LoadEdgePoints failed");
723     QuadDelete( quad );
724     quad = 0;
725     return;
726   }
727   // 3 --- 2D normalized values on unit square [0..1][0..1]
728
729 //  int nbdown = quad->nbPts[0];
730 //  int nbright = quad->nbPts[1];
731   int nbhoriz  = Min(quad->nbPts[0], quad->nbPts[2]);
732   int nbvertic = Min(quad->nbPts[1], quad->nbPts[3]);
733 //  MESSAGE("nbhoriz, nbvertic = " << nbhoriz << nbvertic);
734
735   quad->isEdgeOut[0] = (quad->nbPts[0] > quad->nbPts[2]);
736   quad->isEdgeOut[1] = (quad->nbPts[1] > quad->nbPts[3]);
737   quad->isEdgeOut[2] = (quad->nbPts[2] > quad->nbPts[0]);
738   quad->isEdgeOut[3] = (quad->nbPts[3] > quad->nbPts[1]);
739
740   quad->uv_grid = new UVPtStruct[nbvertic * nbhoriz];
741
742   UVPtStruct *uv_grid = quad->uv_grid;
743   UVPtStruct *uv_e0 = quad->uv_edges[0];
744   UVPtStruct *uv_e1 = quad->uv_edges[1];
745   UVPtStruct *uv_e2 = quad->uv_edges[2];
746   UVPtStruct *uv_e3 = quad->uv_edges[3];
747
748   // nodes Id on "in" edges
749   if (! quad->isEdgeOut[0]) {
750     int j = 0;
751     for (int i = 0; i < nbhoriz; i++) { // down
752       int ij = j * nbhoriz + i;
753       uv_grid[ij].node = uv_e0[i].node;
754     }
755   }
756   if (! quad->isEdgeOut[1]) {
757     int i = nbhoriz - 1;
758     for (int j = 0; j < nbvertic; j++) { // right
759       int ij = j * nbhoriz + i;
760       uv_grid[ij].node = uv_e1[j].node;
761     }
762   }
763   if (! quad->isEdgeOut[2]) {
764     int j = nbvertic - 1;
765     for (int i = 0; i < nbhoriz; i++) { // up
766       int ij = j * nbhoriz + i;
767       uv_grid[ij].node = uv_e2[i].node;
768     }
769   }
770   if (! quad->isEdgeOut[3]) {
771     int i = 0;
772     for (int j = 0; j < nbvertic; j++) { // left
773       int ij = j * nbhoriz + i;
774       uv_grid[ij].node = uv_e3[j].node;
775     }
776   }
777
778   // falsificate "out" edges
779   if (quad->isEdgeOut[0]) // down
780     uv_e0 = MakeEdgePoints
781       (aMesh, F, quad->edge[0], quad->first[0], quad->last[0], nbhoriz - 1);
782   else if (quad->isEdgeOut[2]) // up
783     uv_e2 = MakeEdgePoints
784       (aMesh, F, quad->edge[2], quad->last[2], quad->first[2], nbhoriz - 1);
785
786   if (quad->isEdgeOut[1]) // right
787     uv_e1 = MakeEdgePoints
788       (aMesh, F, quad->edge[1], quad->first[1], quad->last[1], nbvertic - 1);
789   else if (quad->isEdgeOut[3]) // left
790     uv_e3 = MakeEdgePoints
791       (aMesh, F, quad->edge[3], quad->last[3], quad->first[3], nbvertic - 1);
792
793   // normalized 2d values on grid
794   for (int i = 0; i < nbhoriz; i++)
795   {
796     for (int j = 0; j < nbvertic; j++)
797     {
798       int ij = j * nbhoriz + i;
799       // --- droite i cste : x = x0 + y(x1-x0)
800       double x0 = uv_e0[i].normParam;   // bas - sud
801       double x1 = uv_e2[i].normParam;   // haut - nord
802       // --- droite j cste : y = y0 + x(y1-y0)
803       double y0 = uv_e3[j].normParam;   // gauche-ouest
804       double y1 = uv_e1[j].normParam;   // droite - est
805       // --- intersection : x=x0+(y0+x(y1-y0))(x1-x0)
806       double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
807       double y = y0 + x * (y1 - y0);
808       uv_grid[ij].x = x;
809       uv_grid[ij].y = y;
810       //MESSAGE("-xy-01 "<<x0<<" "<<x1<<" "<<y0<<" "<<y1);
811       //MESSAGE("-xy-norm "<<i<<" "<<j<<" "<<x<<" "<<y);
812     }
813   }
814
815   // 4 --- projection on 2d domain (u,v)
816   gp_Pnt2d a0 = pf[0];
817   gp_Pnt2d a1 = pf[1];
818   gp_Pnt2d a2 = pf[2];
819   gp_Pnt2d a3 = pf[3];
820
821   for (int i = 0; i < nbhoriz; i++)
822   {
823     for (int j = 0; j < nbvertic; j++)
824     {
825       int ij = j * nbhoriz + i;
826       double x = uv_grid[ij].x;
827       double y = uv_grid[ij].y;
828       double param_0 = uv_e0[0].param + x * (uv_e0[nbhoriz - 1].param - uv_e0[0].param); // sud
829       double param_2 = uv_e2[0].param + x * (uv_e2[nbhoriz - 1].param - uv_e2[0].param); // nord
830       double param_1 = uv_e1[0].param + y * (uv_e1[nbvertic - 1].param - uv_e1[0].param); // est
831       double param_3 = uv_e3[0].param + y * (uv_e3[nbvertic - 1].param - uv_e3[0].param); // ouest
832
833       //MESSAGE("params "<<param_0<<" "<<param_1<<" "<<param_2<<" "<<param_3);
834       gp_Pnt2d p0 = c2d[0]->Value(param_0);
835       gp_Pnt2d p1 = c2d[1]->Value(param_1);
836       gp_Pnt2d p2 = c2d[2]->Value(param_2);
837       gp_Pnt2d p3 = c2d[3]->Value(param_3);
838
839       double u = (1 - y) * p0.X() + x * p1.X() + y * p2.X() + (1 - x) * p3.X();
840       double v = (1 - y) * p0.Y() + x * p1.Y() + y * p2.Y() + (1 - x) * p3.Y();
841
842       u -= (1 - x) * (1 - y) * a0.X() + x * (1 - y) * a1.X() +
843         x * y * a2.X() + (1 - x) * y * a3.X();
844       v -= (1 - x) * (1 - y) * a0.Y() + x * (1 - y) * a1.Y() +
845         x * y * a2.Y() + (1 - x) * y * a3.Y();
846
847       uv_grid[ij].u = u;
848       uv_grid[ij].v = v;
849
850       //MESSAGE("-uv- "<<i<<" "<<j<<" "<<uv_grid[ij].u<<" "<<uv_grid[ij].v);
851     }
852   }
853 }
854
855 //=============================================================================
856 /*!
857  *  LoadEdgePoints
858  */
859 //=============================================================================
860 UVPtStruct* StdMeshers_Quadrangle_2D::LoadEdgePoints (SMESH_Mesh & aMesh,
861                                                       const TopoDS_Face& F,
862                                                       const TopoDS_Edge& E,
863                                                       double first, double last)
864 //                        bool isForward)
865 {
866   //MESSAGE("StdMeshers_Quadrangle_2D::LoadEdgePoints");
867
868   SMDS_Mesh* meshDS = aMesh.GetMeshDS();
869
870   // --- IDNodes of first and last Vertex
871
872   TopoDS_Vertex VFirst, VLast;
873   TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
874
875   ASSERT(!VFirst.IsNull());
876   SMDS_NodeIteratorPtr lid = aMesh.GetSubMesh(VFirst)->GetSubMeshDS()->GetNodes();
877   if (!lid->more())
878   {
879     MESSAGE ( "NO NODE BUILT ON VERTEX" );
880     return 0;
881   }
882   const SMDS_MeshNode* idFirst = lid->next();
883
884   ASSERT(!VLast.IsNull());
885   lid = aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes();
886   if (!lid->more())
887   {
888     MESSAGE ( "NO NODE BUILT ON VERTEX" );
889     return 0;
890   }
891   const SMDS_MeshNode* idLast = lid->next();
892
893   // --- edge internal IDNodes (relies on good order storage, not checked)
894
895   map<double, const SMDS_MeshNode *> params;
896   SMDS_NodeIteratorPtr ite = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes();
897
898   while(ite->more())
899   {
900     const SMDS_MeshNode* node = ite->next();
901     const SMDS_EdgePosition* epos =
902       static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
903     double param = epos->GetUParameter();
904     params[param] = node;
905   }
906
907   int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
908   if (nbPoints != params.size())
909   {
910     MESSAGE( "BAD NODE ON EDGE POSITIONS" );
911     return 0;
912   }
913   UVPtStruct* uvslf = new UVPtStruct[nbPoints + 2];
914
915   double f, l;
916   Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
917
918   bool isForward = (((l - f) * (last - first)) > 0);
919   double paramin = 0;
920   double paramax = 0;
921   if (isForward)
922   {
923     paramin = f;
924     paramax = l;
925     gp_Pnt2d p = C2d->Value(f); // first point = Vertex Forward
926     uvslf[0].x = p.X();
927     uvslf[0].y = p.Y();
928     uvslf[0].param = f;
929     uvslf[0].node = idFirst;
930     //MESSAGE("__ f "<<f<<" "<<uvslf[0].x <<" "<<uvslf[0].y);
931     map < double, const SMDS_MeshNode* >::iterator itp = params.begin();
932     for (int i = 1; i <= nbPoints; i++) // nbPoints internal
933     {
934       double param = (*itp).first;
935       gp_Pnt2d p = C2d->Value(param);
936       uvslf[i].x = p.X();
937       uvslf[i].y = p.Y();
938       uvslf[i].param = param;
939       uvslf[i].node = (*itp).second;
940       //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[i].x <<" "<<uvslf[i].y);
941       itp++;
942     }
943     p = C2d->Value(l);          // last point = Vertex Reversed
944     uvslf[nbPoints + 1].x = p.X();
945     uvslf[nbPoints + 1].y = p.Y();
946     uvslf[nbPoints + 1].param = l;
947     uvslf[nbPoints + 1].node = idLast;
948     //MESSAGE("__ l "<<l<<" "<<uvslf[nbPoints+1].x <<" "<<uvslf[nbPoints+1].y);
949   } else
950   {
951     paramin = l;
952     paramax = f;
953     gp_Pnt2d p = C2d->Value(l); // first point = Vertex Reversed
954     uvslf[0].x = p.X();
955     uvslf[0].y = p.Y();
956     uvslf[0].param = l;
957     uvslf[0].node = idLast;
958     //MESSAGE("__ l "<<l<<" "<<uvslf[0].x <<" "<<uvslf[0].y);
959     map < double, const SMDS_MeshNode* >::reverse_iterator itp = params.rbegin();
960
961     for (int j = nbPoints; j >= 1; j--) // nbPoints internal
962     {
963       double param = (*itp).first;
964       int i = nbPoints + 1 - j;
965       gp_Pnt2d p = C2d->Value(param);
966       uvslf[i].x = p.X();
967       uvslf[i].y = p.Y();
968       uvslf[i].param = param;
969       uvslf[i].node = (*itp).second;
970       //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[i].x <<" "<<uvslf[i].y);
971       itp++;
972     }
973     p = C2d->Value(f);          // last point = Vertex Forward
974     uvslf[nbPoints + 1].x = p.X();
975     uvslf[nbPoints + 1].y = p.Y();
976     uvslf[nbPoints + 1].param = f;
977     uvslf[nbPoints + 1].node = idFirst;
978     //MESSAGE("__ f "<<f<<" "<<uvslf[nbPoints+1].x <<" "<<uvslf[nbPoints+1].y);
979   }
980
981   ASSERT(paramin != paramax);
982   for (int i = 0; i < nbPoints + 2; i++)
983   {
984     uvslf[i].normParam = (uvslf[i].param - paramin) / (paramax - paramin);
985     //SCRUTE(uvslf[i].normParam);
986   }
987
988   return uvslf;
989 }
990
991 //=============================================================================
992 /*!
993  *  MakeEdgePoints
994  */
995 //=============================================================================
996 UVPtStruct* StdMeshers_Quadrangle_2D::MakeEdgePoints (SMESH_Mesh & aMesh,
997                                                       const TopoDS_Face& F,
998                                                       const TopoDS_Edge& E,
999                                                       double first, double last,
1000                                                       int nb_segm)
1001 {
1002 //  MESSAGE("StdMeshers_Quadrangle_2D::MakeEdgePoints");
1003
1004   UVPtStruct* uvslf = new UVPtStruct[nb_segm + 1];
1005   list<double> params;
1006
1007   // --- edge internal points
1008   double fi, li;
1009   Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, fi, li);
1010   if (!Curve.IsNull()) {
1011     try {
1012       GeomAdaptor_Curve C3d (Curve);
1013       double length = EdgeLength(E);
1014       double eltSize = length / nb_segm;
1015       GCPnts_UniformAbscissa Discret (C3d, eltSize, fi, li);
1016       if (!Discret.IsDone()) return false;
1017       int NbPoints = Discret.NbPoints();
1018       for (int i = 1; i <= NbPoints; i++) {
1019         double param = Discret.Parameter(i);
1020         params.push_back(param);
1021       }
1022     }
1023     catch (Standard_Failure) {
1024       return 0;
1025     }
1026   }
1027   else
1028   {
1029     // Edge is a degenerated Edge
1030     BRep_Tool::Range(E, fi, li);
1031     double du = (li - fi) / nb_segm;
1032     for (int i = 1; i <= nb_segm + 1; i++)
1033     {
1034       double param = fi + (i - 1) * du;
1035       params.push_back(param);
1036     }
1037   }
1038
1039   double f, l;
1040   Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
1041   ASSERT(f != l);
1042
1043   bool isForward = (((l - f) * (last - first)) > 0);
1044   if (isForward) {
1045     list<double>::iterator itU = params.begin();
1046     for (int i = 0; i <= nb_segm; i++) // nbPoints internal
1047     {
1048       double param = *itU;
1049       gp_Pnt2d p = C2d->Value(param);
1050       uvslf[i].x = p.X();
1051       uvslf[i].y = p.Y();
1052       uvslf[i].param = param;
1053       uvslf[i].normParam = (param - f) / (l - f);
1054       itU++;
1055     }
1056   } else {
1057     list<double>::reverse_iterator itU = params.rbegin();
1058     for (int j = nb_segm; j >= 0; j--) // nbPoints internal
1059     {
1060       double param = *itU;
1061       int i = nb_segm - j;
1062       gp_Pnt2d p = C2d->Value(param);
1063       uvslf[i].x = p.X();
1064       uvslf[i].y = p.Y();
1065       uvslf[i].param = param;
1066       uvslf[i].normParam = (param - l) / (f - l);
1067       itU++;
1068     }
1069   }
1070
1071   return uvslf;
1072 }
1073
1074
1075 //=============================================================================
1076 /*!
1077  *  
1078  */
1079 //=============================================================================
1080
1081 ostream & StdMeshers_Quadrangle_2D::SaveTo(ostream & save)
1082 {
1083   return save;
1084 }
1085
1086 //=============================================================================
1087 /*!
1088  *  
1089  */
1090 //=============================================================================
1091
1092 istream & StdMeshers_Quadrangle_2D::LoadFrom(istream & load)
1093 {
1094   return load;
1095 }
1096
1097 //=============================================================================
1098 /*!
1099  *  
1100  */
1101 //=============================================================================
1102
1103 ostream & operator <<(ostream & save, StdMeshers_Quadrangle_2D & hyp)
1104 {
1105   return hyp.SaveTo( save );
1106 }
1107
1108 //=============================================================================
1109 /*!
1110  *  
1111  */
1112 //=============================================================================
1113
1114 istream & operator >>(istream & load, StdMeshers_Quadrangle_2D & hyp)
1115 {
1116   return hyp.LoadFrom( load );
1117 }