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 #include <TopExp.hxx>
  51
  52 #include <Precision.hxx>
  53 #include <gp_Pnt2d.hxx>
  54 #include <TColStd_ListIteratorOfListOfInteger.hxx>
  55 #include <TColStd_SequenceOfReal.hxx>
  56 #include <TColgp_SequenceOfXY.hxx>
  57
  58 #include "utilities.h"
  59 #include "Utils_ExceptHandlers.hxx"
  60
  61 #ifndef StdMeshers_Array2OfNode_HeaderFile
  62 #define StdMeshers_Array2OfNode_HeaderFile
  63 typedef const SMDS_MeshNode* SMDS_MeshNodePtr;
  64 #include <NCollection_DefineArray2.hxx>
  65 DEFINE_BASECOLLECTION (StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr)
  66 DEFINE_ARRAY2(StdMeshers_Array2OfNode,
  67               StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr)
  68 #endif
  69
  70
  71 //=============================================================================
  72 /*!
  73  *
  74  */
  75 //=============================================================================
  76
  77 StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, SMESH_Gen* gen)
  78      : SMESH_2D_Algo(hypId, studyId, gen)
  79 {
  80   MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D");
  81   _name = "Quadrangle_2D";
  82   _shapeType = (1 << TopAbs_FACE);
  83   _compatibleHypothesis.push_back("QuadranglePreference");
  84   myTool = 0;
  85 }
  86
  87 //=============================================================================
  88 /*!
  89  *
  90  */
  91 //=============================================================================
  92
  93 StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D()
  94 {
  95   MESSAGE("StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D");
  96   if ( myTool )
  97     delete myTool;
  98 }
  99
 100 //=============================================================================
 101 /*!
 102  *
 103  */
 104 //=============================================================================
 105
 106 bool StdMeshers_Quadrangle_2D::CheckHypothesis
 107                          (SMESH_Mesh&                          aMesh,
 108                           const TopoDS_Shape&                  aShape,
 109                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 110 {
 111   bool isOk = true;
 112   aStatus = SMESH_Hypothesis::HYP_OK;
 113
 114   // there is only one compatible Hypothesis so far
 115   const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
 116   myQuadranglePreference = hyps.size() > 0;
 117
 118   return isOk;
 119 }
 120
 121 //=============================================================================
 122 /*!
 123  *
 124  */
 125 //=============================================================================
 126
 127 bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh,
 128                                         const TopoDS_Shape& aShape) throw (SALOME_Exception)
 129 {
 130   Unexpect aCatch(SalomeException);
 131   //MESSAGE("StdMeshers_Quadrangle_2D::Compute");
 132   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 133   aMesh.GetSubMesh(aShape);
 134
 135   if ( !myTool )
 136     myTool = new SMESH_MesherHelper(aMesh);
 137   _quadraticMesh = myTool->IsQuadraticSubMesh(aShape);
 138
 139   //FaceQuadStruct *quad = CheckAnd2Dcompute(aMesh, aShape);
 140   FaceQuadStruct* quad = CheckNbEdges(aMesh, aShape);
 141
 142   if (!quad) {
 143     delete myTool; myTool = 0;
 144     return false;
 145   }
 146
 147   if(myQuadranglePreference) {
 148     int n1 = quad->nbPts[0];
 149     int n2 = quad->nbPts[1];
 150     int n3 = quad->nbPts[2];
 151     int n4 = quad->nbPts[3];
 152     int nfull = n1+n2+n3+n4;
 153     int ntmp = nfull/2;
 154     ntmp = ntmp*2;
 155     if( nfull==ntmp && ( (n1!=n3) || (n2!=n4) ) ) {
 156       // special path for using only quandrangle faces
 157       bool ok = ComputeQuadPref(aMesh, aShape, quad);
 158       delete myTool; myTool = 0;
 159       return ok;
 160     }
 161   }
 162
 163   // set normalized grid on unit square in parametric domain
 164   SetNormalizedGrid(aMesh, aShape, quad);
 165   if (!quad) {
 166     delete myTool; myTool = 0;
 167     return false;
 168   }
 169
 170   // --- compute 3D values on points, store points & quadrangles
 171
 172   int nbdown  = quad->nbPts[0];
 173   int nbup    = quad->nbPts[2];
 174
 175   int nbright = quad->nbPts[1];
 176   int nbleft  = quad->nbPts[3];
 177
 178   int nbhoriz  = Min(nbdown, nbup);
 179   int nbvertic = Min(nbright, nbleft);
 180
 181   const TopoDS_Face& F = TopoDS::Face(aShape);
 182   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
 183
 184   // internal mesh nodes
 185   int i, j, geomFaceID = meshDS->ShapeToIndex( F );
 186   for (i = 1; i < nbhoriz - 1; i++) {
 187     for (j = 1; j < nbvertic - 1; j++) {
 188       int ij = j * nbhoriz + i;
 189       double u = quad->uv_grid[ij].u;
 190       double v = quad->uv_grid[ij].v;
 191       gp_Pnt P = S->Value(u, v);
 192       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 193       meshDS->SetNodeOnFace(node, geomFaceID, u, v);
 194       quad->uv_grid[ij].node = node;
 195     }
 196   }
 197
 198   // mesh faces
 199
 200   //             [2]
 201   //      --.--.--.--.--.--  nbvertic
 202   //     |                 | ^
 203   //     |                 | ^
 204   // [3] |                 | ^ j  [1]
 205   //     |                 | ^
 206   //     |                 | ^
 207   //      ---.----.----.---  0
 208   //     0 > > > > > > > > nbhoriz
 209   //              i
 210   //             [0]
 211
 212   i = 0;
 213   int ilow = 0;
 214   int iup = nbhoriz - 1;
 215   if (quad->isEdgeOut[3]) { ilow++; } else { if (quad->isEdgeOut[1]) iup--; }
 216
 217   int jlow = 0;
 218   int jup = nbvertic - 1;
 219   if (quad->isEdgeOut[0]) { jlow++; } else { if (quad->isEdgeOut[2]) jup--; }
 220
 221   // regular quadrangles
 222   for (i = ilow; i < iup; i++) {
 223     for (j = jlow; j < jup; j++) {
 224       const SMDS_MeshNode *a, *b, *c, *d;
 225       a = quad->uv_grid[j * nbhoriz + i].node;
 226       b = quad->uv_grid[j * nbhoriz + i + 1].node;
 227       c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node;
 228       d = quad->uv_grid[(j + 1) * nbhoriz + i].node;
 229       //SMDS_MeshFace * face = meshDS->AddFace(a, b, c, d);
 230       SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
 231       meshDS->SetMeshElementOnShape(face, geomFaceID);
 232     }
 233   }
 234
 235   UVPtStruct *uv_e0 = quad->uv_edges[0];
 236   UVPtStruct *uv_e1 = quad->uv_edges[1];
 237   UVPtStruct *uv_e2 = quad->uv_edges[2];
 238   UVPtStruct *uv_e3 = quad->uv_edges[3];
 239
 240   double eps = Precision::Confusion();
 241
 242   // Boundary quadrangles
 243
 244   if (quad->isEdgeOut[0]) {
 245     // Down edge is out
 246     //
 247     // |___|___|___|___|___|___|
 248     // |   |   |   |   |   |   |
 249     // |___|___|___|___|___|___|
 250     // |   |   |   |   |   |   |
 251     // |___|___|___|___|___|___| __ first row of the regular grid
 252     // .  .  .  .  .  .  .  .  . __ down edge nodes
 253     //
 254     // >->->->->->->->->->->->-> -- direction of processing
 255
 256     int g = 0; // number of last processed node in the regular grid
 257
 258     // number of last node of the down edge to be processed
 259     int stop = nbdown - 1;
 260     // if right edge is out, we will stop at a node, previous to the last one
 261     if (quad->isEdgeOut[1]) stop--;
 262
 263     // for each node of the down edge find nearest node
 264     // in the first row of the regular grid and link them
 265     for (i = 0; i < stop; i++) {
 266       const SMDS_MeshNode *a, *b, *c, *d;
 267       a = uv_e0[i].node;
 268       b = uv_e0[i + 1].node;
 269       gp_Pnt pb (b->X(), b->Y(), b->Z());
 270
 271       // find node c in the regular grid, which will be linked with node b
 272       int near = g;
 273       if (i == stop - 1) {
 274         // right bound reached, link with the rightmost node
 275         near = iup;
 276         c = quad->uv_grid[nbhoriz + iup].node;
 277       }
 278       else {
 279         // find in the grid node c, nearest to the b
 280         double mind = RealLast();
 281         for (int k = g; k <= iup; k++) {
 282
 283           const SMDS_MeshNode *nk;
 284           if (k < ilow) // this can be, if left edge is out
 285             nk = uv_e3[1].node; // get node from the left edge
 286           else
 287             nk = quad->uv_grid[nbhoriz + k].node; // get one of middle nodes
 288
 289           gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 290           double dist = pb.Distance(pnk);
 291           if (dist < mind - eps) {
 292             c = nk;
 293             near = k;
 294             mind = dist;
 295           } else {
 296             break;
 297           }
 298         }
 299       }
 300
 301       if (near == g) { // make triangle
 302         //SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
 303         SMDS_MeshFace* face = myTool->AddFace(a, b, c);
 304         meshDS->SetMeshElementOnShape(face, geomFaceID);
 305       }
 306       else { // make quadrangle
 307         if (near - 1 < ilow)
 308           d = uv_e3[1].node;
 309         else
 310           d = quad->uv_grid[nbhoriz + near - 1].node;
 311         //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 312         SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
 313         meshDS->SetMeshElementOnShape(face, geomFaceID);
 314
 315         // if node d is not at position g - make additional triangles
 316         if (near - 1 > g) {
 317           for (int k = near - 1; k > g; k--) {
 318             c = quad->uv_grid[nbhoriz + k].node;
 319             if (k - 1 < ilow)
 320               d = uv_e3[1].node;
 321             else
 322               d = quad->uv_grid[nbhoriz + k - 1].node;
 323             //SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
 324             SMDS_MeshFace* face = myTool->AddFace(a, c, d);
 325             meshDS->SetMeshElementOnShape(face, geomFaceID);
 326           }
 327         }
 328         g = near;
 329       }
 330     }
 331   } else {
 332     if (quad->isEdgeOut[2]) {
 333       // Up edge is out
 334       //
 335       // <-<-<-<-<-<-<-<-<-<-<-<-< -- direction of processing
 336       //
 337       // .  .  .  .  .  .  .  .  . __ up edge nodes
 338       //  ___ ___ ___ ___ ___ ___  __ first row of the regular grid
 339       // |   |   |   |   |   |   |
 340       // |___|___|___|___|___|___|
 341       // |   |   |   |   |   |   |
 342       // |___|___|___|___|___|___|
 343       // |   |   |   |   |   |   |
 344
 345       int g = nbhoriz - 1; // last processed node in the regular grid
 346
 347       int stop = 0;
 348       // if left edge is out, we will stop at a second node
 349       if (quad->isEdgeOut[3]) stop++;
 350
 351       // for each node of the up edge find nearest node
 352       // in the first row of the regular grid and link them
 353       for (i = nbup - 1; i > stop; i--) {
 354         const SMDS_MeshNode *a, *b, *c, *d;
 355         a = uv_e2[i].node;
 356         b = uv_e2[i - 1].node;
 357         gp_Pnt pb (b->X(), b->Y(), b->Z());
 358
 359         // find node c in the grid, which will be linked with node b
 360         int near = g;
 361         if (i == stop + 1) { // left bound reached, link with the leftmost node
 362           c = quad->uv_grid[nbhoriz*(nbvertic - 2) + ilow].node;
 363           near = ilow;
 364         } else {
 365           // find node c in the grid, nearest to the b
 366           double mind = RealLast();
 367           for (int k = g; k >= ilow; k--) {
 368             const SMDS_MeshNode *nk;
 369             if (k > iup)
 370               nk = uv_e1[nbright - 2].node;
 371             else
 372               nk = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
 373             gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 374             double dist = pb.Distance(pnk);
 375             if (dist < mind - eps) {
 376               c = nk;
 377               near = k;
 378               mind = dist;
 379             } else {
 380               break;
 381             }
 382           }
 383         }
 384
 385         if (near == g) { // make triangle
 386           //SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
 387           SMDS_MeshFace* face = myTool->AddFace(a, b, c);
 388           meshDS->SetMeshElementOnShape(face, geomFaceID);
 389         }
 390         else { // make quadrangle
 391           if (near + 1 > iup)
 392             d = uv_e1[nbright - 2].node;
 393           else
 394             d = quad->uv_grid[nbhoriz*(nbvertic - 2) + near + 1].node;
 395           //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 396           SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
 397           meshDS->SetMeshElementOnShape(face, geomFaceID);
 398
 399           if (near + 1 < g) { // if d not is at g - make additional triangles
 400             for (int k = near + 1; k < g; k++) {
 401               c = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
 402               if (k + 1 > iup)
 403                 d = uv_e1[nbright - 2].node;
 404               else
 405                 d = quad->uv_grid[nbhoriz*(nbvertic - 2) + k + 1].node;
 406               //SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
 407               SMDS_MeshFace* face = myTool->AddFace(a, c, d);
 408               meshDS->SetMeshElementOnShape(face, geomFaceID);
 409             }
 410           }
 411           g = near;
 412         }
 413       }
 414     }
 415   }
 416
 417   // right or left boundary quadrangles
 418   if (quad->isEdgeOut[1]) {
 419 //    MESSAGE("right edge is out");
 420     int g = 0; // last processed node in the grid
 421     int stop = nbright - 1;
 422     if (quad->isEdgeOut[2]) stop--;
 423     for (i = 0; i < stop; i++) {
 424       const SMDS_MeshNode *a, *b, *c, *d;
 425       a = uv_e1[i].node;
 426       b = uv_e1[i + 1].node;
 427       gp_Pnt pb (b->X(), b->Y(), b->Z());
 428
 429       // find node c in the grid, nearest to the b
 430       int near = g;
 431       if (i == stop - 1) { // up bondary reached
 432         c = quad->uv_grid[nbhoriz*(jup + 1) - 2].node;
 433         near = jup;
 434       } else {
 435         double mind = RealLast();
 436         for (int k = g; k <= jup; k++) {
 437           const SMDS_MeshNode *nk;
 438           if (k < jlow)
 439             nk = uv_e0[nbdown - 2].node;
 440           else
 441             nk = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
 442           gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 443           double dist = pb.Distance(pnk);
 444           if (dist < mind - eps) {
 445             c = nk;
 446             near = k;
 447             mind = dist;
 448           } else {
 449             break;
 450           }
 451         }
 452       }
 453
 454       if (near == g) { // make triangle
 455         //SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
 456         SMDS_MeshFace* face = myTool->AddFace(a, b, c);
 457         meshDS->SetMeshElementOnShape(face, geomFaceID);
 458       }
 459       else { // make quadrangle
 460         if (near - 1 < jlow)
 461           d = uv_e0[nbdown - 2].node;
 462         else
 463           d = quad->uv_grid[nbhoriz*near - 2].node;
 464         //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 465         SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
 466         meshDS->SetMeshElementOnShape(face, geomFaceID);
 467
 468         if (near - 1 > g) { // if d not is at g - make additional triangles
 469           for (int k = near - 1; k > g; k--) {
 470             c = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
 471             if (k - 1 < jlow)
 472               d = uv_e0[nbdown - 2].node;
 473             else
 474               d = quad->uv_grid[nbhoriz*k - 2].node;
 475             //SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
 476             SMDS_MeshFace* face = myTool->AddFace(a, c, d);
 477             meshDS->SetMeshElementOnShape(face, geomFaceID);
 478           }
 479         }
 480         g = near;
 481       }
 482     }
 483   } else {
 484     if (quad->isEdgeOut[3]) {
 485 //      MESSAGE("left edge is out");
 486       int g = nbvertic - 1; // last processed node in the grid
 487       int stop = 0;
 488       if (quad->isEdgeOut[0]) stop++;
 489       for (i = nbleft - 1; i > stop; i--) {
 490         const SMDS_MeshNode *a, *b, *c, *d;
 491         a = uv_e3[i].node;
 492         b = uv_e3[i - 1].node;
 493         gp_Pnt pb (b->X(), b->Y(), b->Z());
 494
 495         // find node c in the grid, nearest to the b
 496         int near = g;
 497         if (i == stop + 1) { // down bondary reached
 498           c = quad->uv_grid[nbhoriz*jlow + 1].node;
 499           near = jlow;
 500         } else {
 501           double mind = RealLast();
 502           for (int k = g; k >= jlow; k--) {
 503             const SMDS_MeshNode *nk;
 504             if (k > jup)
 505               nk = uv_e2[1].node;
 506             else
 507               nk = quad->uv_grid[nbhoriz*k + 1].node;
 508             gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 509             double dist = pb.Distance(pnk);
 510             if (dist < mind - eps) {
 511               c = nk;
 512               near = k;
 513               mind = dist;
 514             } else {
 515               break;
 516             }
 517           }
 518         }
 519
 520         if (near == g) { // make triangle
 521           //SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
 522           SMDS_MeshFace* face = myTool->AddFace(a, b, c);
 523           meshDS->SetMeshElementOnShape(face, geomFaceID);
 524         }
 525         else { // make quadrangle
 526           if (near + 1 > jup)
 527             d = uv_e2[1].node;
 528           else
 529             d = quad->uv_grid[nbhoriz*(near + 1) + 1].node;
 530           //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 531           SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
 532           meshDS->SetMeshElementOnShape(face, geomFaceID);
 533
 534           if (near + 1 < g) { // if d not is at g - make additional triangles
 535             for (int k = near + 1; k < g; k++) {
 536               c = quad->uv_grid[nbhoriz*k + 1].node;
 537               if (k + 1 > jup)
 538                 d = uv_e2[1].node;
 539               else
 540                 d = quad->uv_grid[nbhoriz*(k + 1) + 1].node;
 541               //SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
 542               SMDS_MeshFace* face = myTool->AddFace(a, c, d);
 543               meshDS->SetMeshElementOnShape(face, geomFaceID);
 544             }
 545           }
 546           g = near;
 547         }
 548       }
 549     }
 550   }
 551
 552   QuadDelete(quad);
 553   delete myTool; myTool = 0;
 554
 555   bool isOk = true;
 556   return isOk;
 557 }
 558
 559
 560 //=============================================================================
 561 /*!
 562  *
 563  */
 564 //=============================================================================
 565
 566 FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMesh,
 567                                                        const TopoDS_Shape & aShape)
 568      throw(SALOME_Exception)
 569 {
 570   Unexpect aCatch(SalomeException);
 571
 572   const TopoDS_Face & F = TopoDS::Face(aShape);
 573
 574   // verify 1 wire only, with 4 edges
 575
 576   if (NumberOfWires(F) != 1) {
 577     INFOS("only 1 wire by face (quadrangles)");
 578     return 0;
 579   }
 580   const TopoDS_Wire& W = BRepTools::OuterWire(F);
 581   BRepTools_WireExplorer wexp (W, F);
 582
 583   FaceQuadStruct* quad = new FaceQuadStruct;
 584   for (int i = 0; i < 4; i++)
 585     quad->uv_edges[i] = 0;
 586   quad->uv_grid = 0;
 587
 588   int nbEdges = 0;
 589   for (wexp.Init(W, F); wexp.More(); wexp.Next()) {
 590     const TopoDS_Edge& E = wexp.Current();
 591     int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
 592     if (nbEdges < 4) {
 593       quad->edge[nbEdges] = E;
 594       if(!_quadraticMesh) {
 595         quad->nbPts[nbEdges] = nb + 2; // internal points + 2 extrema
 596       }
 597       else {
 598         int tmp = nb/2;
 599         quad->nbPts[nbEdges] = tmp + 2; // internal not medium points + 2 extrema
 600       }
 601     }
 602     nbEdges++;
 603   }
 604
 605   if (nbEdges != 4) {
 606     INFOS("face must have 4 edges /quadrangles");
 607     QuadDelete(quad);
 608     return 0;
 609   }
 610
 611   return quad;
 612 }
 613
 614 //=============================================================================
 615 /*!
 616  *  CheckAnd2Dcompute
 617  */
 618 //=============================================================================
 619
 620 FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute
 621                            (SMESH_Mesh & aMesh,
 622                             const TopoDS_Shape & aShape,
 623                             const bool CreateQuadratic) throw(SALOME_Exception)
 624 {
 625   Unexpect aCatch(SalomeException);
 626
 627   _quadraticMesh = CreateQuadratic;
 628
 629   FaceQuadStruct *quad = CheckNbEdges(aMesh, aShape);
 630
 631   if(!quad) return 0;
 632
 633   // set normalized grid on unit square in parametric domain
 634   SetNormalizedGrid(aMesh, aShape, quad);
 635
 636   return quad;
 637 }
 638
 639 //=============================================================================
 640 /*!
 641  *
 642  */
 643 //=============================================================================
 644
 645 void StdMeshers_Quadrangle_2D::QuadDelete (FaceQuadStruct * quad)
 646 {
 647   //MESSAGE("StdMeshers_Quadrangle_2D::QuadDelete");
 648   if (quad)
 649   {
 650     for (int i = 0; i < 4; i++)
 651     {
 652       if (quad->uv_edges[i])
 653         delete [] quad->uv_edges[i];
 654       quad->edge[i].Nullify();
 655     }
 656     if (quad->uv_grid)
 657       delete [] quad->uv_grid;
 658     delete quad;
 659   }
 660 }
 661
 662 //=============================================================================
 663 /*!
 664  *
 665  */
 666 //=============================================================================
 667
 668 void StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh,
 669                                                   const TopoDS_Shape& aShape,
 670                                                   FaceQuadStruct* quad) throw (SALOME_Exception)
 671 {
 672   Unexpect aCatch(SalomeException);
 673   // Algorithme décrit dans "Génération automatique de maillages"
 674   // P.L. GEORGE, MASSON, § 6.4.1 p. 84-85
 675   // traitement dans le domaine paramétrique 2d u,v
 676   // transport - projection sur le carré unité
 677
 678 //  MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid");
 679   const TopoDS_Face& F = TopoDS::Face(aShape);
 680
 681   // 1 --- find orientation of the 4 edges, by test on extrema
 682
 683   //      max             min                    0     x1     1
 684   //     |<----north-2-------^                a3 -------------> a2
 685   //     |                   |                   ^1          1^
 686   //    west-3            east-1 =right          |            |
 687   //     |                   |         ==>       |            |
 688   //  y0 |                   | y1                |            |
 689   //     |                   |                   |0          0|
 690   //     v----south-0-------->                a0 -------------> a1
 691   //      min             max                    0     x0     1
 692   //             =down
 693   //
 694
 695   Handle(Geom2d_Curve) c2d[4];
 696   gp_Pnt2d pf[4];
 697   gp_Pnt2d pl[4];
 698   for (int i = 0; i < 4; i++)
 699   {
 700     c2d[i] = BRep_Tool::CurveOnSurface(quad->edge[i], F,
 701                                        quad->first[i], quad->last[i]);
 702     pf[i] = c2d[i]->Value(quad->first[i]);
 703     pl[i] = c2d[i]->Value(quad->last[i]);
 704     quad->isEdgeForward[i] = false;
 705   }
 706
 707   double l0f1 = pl[0].SquareDistance(pf[1]);
 708   double l0l1 = pl[0].SquareDistance(pl[1]);
 709   double f0f1 = pf[0].SquareDistance(pf[1]);
 710   double f0l1 = pf[0].SquareDistance(pl[1]);
 711   if ( Min( l0f1, l0l1 ) < Min ( f0f1, f0l1 ))
 712   {
 713     quad->isEdgeForward[0] = true;
 714   } else {
 715     double tmp = quad->first[0];
 716     quad->first[0] = quad->last[0];
 717     quad->last[0] = tmp;
 718     pf[0] = c2d[0]->Value(quad->first[0]);
 719     pl[0] = c2d[0]->Value(quad->last[0]);
 720   }
 721   for (int i = 1; i < 4; i++)
 722   {
 723     l0l1 = pl[i - 1].SquareDistance(pl[i]);
 724     l0f1 = pl[i - 1].SquareDistance(pf[i]);
 725     quad->isEdgeForward[i] = ( l0f1 < l0l1 );
 726     if (!quad->isEdgeForward[i])
 727     {
 728       double tmp = quad->first[i];
 729       quad->first[i] = quad->last[i];
 730       quad->last[i] = tmp;
 731       pf[i] = c2d[i]->Value(quad->first[i]);
 732       pl[i] = c2d[i]->Value(quad->last[i]);
 733     }
 734   }
 735
 736   // 2 --- load 2d edge points (u,v) with orientation and value on unit square
 737
 738   bool loadOk = true;
 739   for (int i = 0; i < 2; i++)
 740   {
 741     quad->uv_edges[i] = LoadEdgePoints(aMesh, F, quad->edge[i],
 742                                        quad->first[i], quad->last[i]);
 743     if (!quad->uv_edges[i]) loadOk = false;
 744   }
 745
 746   for (int i = 2; i < 4; i++)
 747   {
 748     quad->uv_edges[i] = LoadEdgePoints(aMesh, F, quad->edge[i],
 749                                        quad->last[i], quad->first[i]);
 750     if (!quad->uv_edges[i]) loadOk = false;
 751   }
 752
 753   if (!loadOk)
 754   {
 755     INFOS("StdMeshers_Quadrangle_2D::SetNormalizedGrid - LoadEdgePoints failed");
 756     QuadDelete( quad );
 757     quad = 0;
 758     return;
 759   }
 760   // 3 --- 2D normalized values on unit square [0..1][0..1]
 761
 762   int nbhoriz  = Min(quad->nbPts[0], quad->nbPts[2]);
 763   int nbvertic = Min(quad->nbPts[1], quad->nbPts[3]);
 764
 765   quad->isEdgeOut[0] = (quad->nbPts[0] > quad->nbPts[2]);
 766   quad->isEdgeOut[1] = (quad->nbPts[1] > quad->nbPts[3]);
 767   quad->isEdgeOut[2] = (quad->nbPts[2] > quad->nbPts[0]);
 768   quad->isEdgeOut[3] = (quad->nbPts[3] > quad->nbPts[1]);
 769
 770   quad->uv_grid = new UVPtStruct[nbvertic * nbhoriz];
 771
 772   UVPtStruct *uv_grid = quad->uv_grid;
 773   UVPtStruct *uv_e0 = quad->uv_edges[0];
 774   UVPtStruct *uv_e1 = quad->uv_edges[1];
 775   UVPtStruct *uv_e2 = quad->uv_edges[2];
 776   UVPtStruct *uv_e3 = quad->uv_edges[3];
 777
 778   // nodes Id on "in" edges
 779   if (! quad->isEdgeOut[0]) {
 780     int j = 0;
 781     for (int i = 0; i < nbhoriz; i++) { // down
 782       int ij = j * nbhoriz + i;
 783       uv_grid[ij].node = uv_e0[i].node;
 784     }
 785   }
 786   if (! quad->isEdgeOut[1]) {
 787     int i = nbhoriz - 1;
 788     for (int j = 0; j < nbvertic; j++) { // right
 789       int ij = j * nbhoriz + i;
 790       uv_grid[ij].node = uv_e1[j].node;
 791     }
 792   }
 793   if (! quad->isEdgeOut[2]) {
 794     int j = nbvertic - 1;
 795     for (int i = 0; i < nbhoriz; i++) { // up
 796       int ij = j * nbhoriz + i;
 797       uv_grid[ij].node = uv_e2[i].node;
 798     }
 799   }
 800   if (! quad->isEdgeOut[3]) {
 801     int i = 0;
 802     for (int j = 0; j < nbvertic; j++) { // left
 803       int ij = j * nbhoriz + i;
 804       uv_grid[ij].node = uv_e3[j].node;
 805     }
 806   }
 807
 808   // falsificate "out" edges
 809   if (quad->isEdgeOut[0]) // down
 810     uv_e0 = MakeEdgePoints
 811       (aMesh, F, quad->edge[0], quad->first[0], quad->last[0], nbhoriz - 1);
 812   else if (quad->isEdgeOut[2]) // up
 813     uv_e2 = MakeEdgePoints
 814       (aMesh, F, quad->edge[2], quad->last[2], quad->first[2], nbhoriz - 1);
 815
 816   if (quad->isEdgeOut[1]) // right
 817     uv_e1 = MakeEdgePoints
 818       (aMesh, F, quad->edge[1], quad->first[1], quad->last[1], nbvertic - 1);
 819   else if (quad->isEdgeOut[3]) // left
 820     uv_e3 = MakeEdgePoints
 821       (aMesh, F, quad->edge[3], quad->last[3], quad->first[3], nbvertic - 1);
 822
 823   // normalized 2d values on grid
 824   for (int i = 0; i < nbhoriz; i++)
 825   {
 826     for (int j = 0; j < nbvertic; j++)
 827     {
 828       int ij = j * nbhoriz + i;
 829       // --- droite i cste : x = x0 + y(x1-x0)
 830       double x0 = uv_e0[i].normParam;   // bas - sud
 831       double x1 = uv_e2[i].normParam;   // haut - nord
 832       // --- droite j cste : y = y0 + x(y1-y0)
 833       double y0 = uv_e3[j].normParam;   // gauche-ouest
 834       double y1 = uv_e1[j].normParam;   // droite - est
 835       // --- intersection : x=x0+(y0+x(y1-y0))(x1-x0)
 836       double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
 837       double y = y0 + x * (y1 - y0);
 838       uv_grid[ij].x = x;
 839       uv_grid[ij].y = y;
 840       //MESSAGE("-xy-01 "<<x0<<" "<<x1<<" "<<y0<<" "<<y1);
 841       //MESSAGE("-xy-norm "<<i<<" "<<j<<" "<<x<<" "<<y);
 842     }
 843   }
 844
 845   // 4 --- projection on 2d domain (u,v)
 846   gp_Pnt2d a0 = pf[0];
 847   gp_Pnt2d a1 = pf[1];
 848   gp_Pnt2d a2 = pf[2];
 849   gp_Pnt2d a3 = pf[3];
 850
 851   for (int i = 0; i < nbhoriz; i++)
 852   {
 853     for (int j = 0; j < nbvertic; j++)
 854     {
 855       int ij = j * nbhoriz + i;
 856       double x = uv_grid[ij].x;
 857       double y = uv_grid[ij].y;
 858       double param_0 = uv_e0[0].param + x * (uv_e0[nbhoriz - 1].param - uv_e0[0].param); // sud
 859       double param_2 = uv_e2[0].param + x * (uv_e2[nbhoriz - 1].param - uv_e2[0].param); // nord
 860       double param_1 = uv_e1[0].param + y * (uv_e1[nbvertic - 1].param - uv_e1[0].param); // est
 861       double param_3 = uv_e3[0].param + y * (uv_e3[nbvertic - 1].param - uv_e3[0].param); // ouest
 862
 863       //MESSAGE("params "<<param_0<<" "<<param_1<<" "<<param_2<<" "<<param_3);
 864       gp_Pnt2d p0 = c2d[0]->Value(param_0);
 865       gp_Pnt2d p1 = c2d[1]->Value(param_1);
 866       gp_Pnt2d p2 = c2d[2]->Value(param_2);
 867       gp_Pnt2d p3 = c2d[3]->Value(param_3);
 868
 869       double u = (1 - y) * p0.X() + x * p1.X() + y * p2.X() + (1 - x) * p3.X();
 870       double v = (1 - y) * p0.Y() + x * p1.Y() + y * p2.Y() + (1 - x) * p3.Y();
 871
 872       u -= (1 - x) * (1 - y) * a0.X() + x * (1 - y) * a1.X() +
 873         x * y * a2.X() + (1 - x) * y * a3.X();
 874       v -= (1 - x) * (1 - y) * a0.Y() + x * (1 - y) * a1.Y() +
 875         x * y * a2.Y() + (1 - x) * y * a3.Y();
 876
 877       uv_grid[ij].u = u;
 878       uv_grid[ij].v = v;
 879     }
 880   }
 881 }
 882
 883
 884 //=======================================================================
 885 //function : ShiftQuad
 886 //purpose  : auxilary function for ComputeQuadPref
 887 //=======================================================================
 888 static void ShiftQuad(FaceQuadStruct* quad, const int num, bool WisF)
 889 {
 890   if(num>3) return;
 891   int i;
 892   for(i=1; i<=num; i++) {
 893     int nbPts3 = quad->nbPts[0];
 894     quad->nbPts[0] = quad->nbPts[1];
 895     quad->nbPts[1] = quad->nbPts[2];
 896     quad->nbPts[2] = quad->nbPts[3];
 897     quad->nbPts[3] = nbPts3;
 898     TopoDS_Edge edge3 = quad->edge[0];
 899     quad->edge[0] = quad->edge[1];
 900     quad->edge[1] = quad->edge[2];
 901     quad->edge[2] = quad->edge[3];
 902     quad->edge[3] = edge3;
 903     double first3 = quad->first[0];
 904     quad->first[0] = quad->first[1];
 905     quad->first[1] = quad->first[2];
 906     quad->first[2] = quad->first[3];
 907     quad->first[3] = first3;
 908     double last3 = quad->last[0];
 909     quad->last[0] = quad->last[1];
 910     quad->last[1] = quad->last[2];
 911     quad->last[2] = quad->last[3];
 912     quad->last[3] = last3;
 913     bool isEdgeForward3 = quad->isEdgeForward[0];
 914     quad->isEdgeForward[0] = quad->isEdgeForward[1];
 915     quad->isEdgeForward[1] = quad->isEdgeForward[2];
 916     quad->isEdgeForward[2] = quad->isEdgeForward[3];
 917     quad->isEdgeForward[3] = isEdgeForward3;
 918     bool isEdgeOut3 = quad->isEdgeOut[0];
 919     quad->isEdgeOut[0] = quad->isEdgeOut[1];
 920     quad->isEdgeOut[1] = quad->isEdgeOut[2];
 921     quad->isEdgeOut[2] = quad->isEdgeOut[3];
 922     quad->isEdgeOut[3] = isEdgeOut3;
 923     UVPtStruct* uv_edges3 = quad->uv_edges[0];
 924     quad->uv_edges[0] = quad->uv_edges[1];
 925     quad->uv_edges[1] = quad->uv_edges[2];
 926     quad->uv_edges[2] = quad->uv_edges[3];
 927     quad->uv_edges[3] = uv_edges3;
 928   }
 929   if(!WisF) {
 930     // replacement left and right edges
 931     int nbPts3 = quad->nbPts[1];
 932     quad->nbPts[1] = quad->nbPts[3];
 933     quad->nbPts[3] = nbPts3;
 934     TopoDS_Edge edge3 = quad->edge[1];
 935     quad->edge[1] = quad->edge[3];
 936     quad->edge[3] = edge3;
 937     double first3 = quad->first[1];
 938     quad->first[1] = quad->first[3];
 939     quad->first[3] = first3;
 940     double last3 = quad->last[1];
 941     quad->last[1] = quad->last[2];
 942     quad->last[3] = last3;
 943     bool isEdgeForward3 = quad->isEdgeForward[1];
 944     quad->isEdgeForward[1] = quad->isEdgeForward[3];
 945     quad->isEdgeForward[3] = isEdgeForward3;
 946     bool isEdgeOut3 = quad->isEdgeOut[1];
 947     quad->isEdgeOut[1] = quad->isEdgeOut[3];
 948     quad->isEdgeOut[3] = isEdgeOut3;
 949     UVPtStruct* uv_edges3 = quad->uv_edges[1];
 950     quad->uv_edges[1] = quad->uv_edges[3];
 951     quad->uv_edges[3] = uv_edges3;
 952   }
 953 }
 954
 955
 956 //=======================================================================
 957 //function : CalcUV
 958 //purpose  : auxilary function for ComputeQuadPref
 959 //=======================================================================
 960 static gp_XY CalcUV(double x0, double x1, double y0, double y1,
 961                     FaceQuadStruct* quad,
 962                     const gp_Pnt2d& a0, const gp_Pnt2d& a1,
 963                     const gp_Pnt2d& a2, const gp_Pnt2d& a3,
 964                     const Handle(Geom2d_Curve)& c2db,
 965                     const Handle(Geom2d_Curve)& c2dr,
 966                     const Handle(Geom2d_Curve)& c2dt,
 967                     const Handle(Geom2d_Curve)& c2dl)
 968 {
 969   int nb = quad->nbPts[0];
 970   int nr = quad->nbPts[1];
 971   int nt = quad->nbPts[2];
 972   int nl = quad->nbPts[3];
 973
 974   UVPtStruct* uv_eb = quad->uv_edges[0];
 975   UVPtStruct* uv_er = quad->uv_edges[1];
 976   UVPtStruct* uv_et = quad->uv_edges[2];
 977   UVPtStruct* uv_el = quad->uv_edges[3];
 978
 979   double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
 980   double y = y0 + x * (y1 - y0);
 981
 982   double param_b = uv_eb[0].param + x * (uv_eb[nb-1].param - uv_eb[0].param);
 983   double param_t = uv_et[0].param + x * (uv_et[nt-1].param - uv_et[0].param);
 984   double param_r = uv_er[0].param + y * (uv_er[nr-1].param - uv_er[0].param);
 985   double param_l = uv_el[0].param + y * (uv_el[nl-1].param - uv_el[0].param);
 986
 987   gp_Pnt2d p0 = c2db->Value(param_b);
 988   gp_Pnt2d p1 = c2dr->Value(param_r);
 989   gp_Pnt2d p2 = c2dt->Value(param_t);
 990   gp_Pnt2d p3 = c2dl->Value(param_l);
 991
 992   double u = (1 - y) * p0.X() + x * p1.X() + y * p2.X() + (1 - x) * p3.X();
 993   double v = (1 - y) * p0.Y() + x * p1.Y() + y * p2.Y() + (1 - x) * p3.Y();
 994
 995   u -= (1 - x) * (1 - y) * a0.X() + x * (1 - y) * a1.X() +
 996     x * y * a2.X() + (1 - x) * y * a3.X();
 997   v -= (1 - x) * (1 - y) * a0.Y() + x * (1 - y) * a1.Y() +
 998     x * y * a2.Y() + (1 - x) * y * a3.Y();
 999
1000   //cout<<"x0="<<x0<<" x1="<<x1<<" y0="<<y0<<" y1="<<y1<<endl;
1001   //cout<<"x="<<x<<" y="<<y<<endl;
1002   //cout<<"param_b="<<param_b<<" param_t="<<param_t<<" param_r="<<param_r<<" param_l="<<param_l<<endl;
1003   //cout<<"u="<<u<<" v="<<v<<endl;
1004
1005   return gp_XY(u,v);
1006 }
1007
1008
1009 //=======================================================================
1010 //function : ComputeQuadPref
1011 //purpose  :
1012 //=======================================================================
1013 /*!
1014  * Special function for creation only quandrangle faces
1015  */
1016 bool StdMeshers_Quadrangle_2D::ComputeQuadPref
1017                           (SMESH_Mesh & aMesh,
1018                            const TopoDS_Shape& aShape,
1019                            FaceQuadStruct* quad) throw (SALOME_Exception)
1020 {
1021   Unexpect aCatch(SalomeException);
1022
1023   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
1024   const TopoDS_Face& F = TopoDS::Face(aShape);
1025   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
1026   const TopoDS_Wire& W = BRepTools::OuterWire(F);
1027   bool WisF = false;
1028   if(W.Orientation()==TopAbs_FORWARD)
1029     WisF = true;
1030   //if(WisF) cout<<"W is FORWARD"<<endl;
1031   //else cout<<"W is REVERSED"<<endl;
1032   bool FisF = (F.Orientation()==TopAbs_FORWARD);
1033   if(!FisF) WisF = !WisF;
1034   int i,j,geomFaceID = meshDS->ShapeToIndex( F );
1035
1036   int nb = quad->nbPts[0];
1037   int nr = quad->nbPts[1];
1038   int nt = quad->nbPts[2];
1039   int nl = quad->nbPts[3];
1040   int dh = abs(nb-nt);
1041   int dv = abs(nr-nl);
1042
1043   if( dh>=dv ) {
1044     if( nt>nb ) {
1045       // it is a base case => not shift quad but me be replacement is need
1046       ShiftQuad(quad,0,WisF);
1047     }
1048     else {
1049       // we have to shift quad on 2
1050       ShiftQuad(quad,2,WisF);
1051     }
1052   }
1053   else {
1054     if( nr>nl ) {
1055       // we have to shift quad on 3
1056       ShiftQuad(quad,3,WisF);
1057     }
1058     else {
1059       // we have to shift quad on 1
1060       ShiftQuad(quad,1,WisF);
1061     }
1062   }
1063
1064   nb = quad->nbPts[0];
1065   nr = quad->nbPts[1];
1066   nt = quad->nbPts[2];
1067   nl = quad->nbPts[3];
1068   dh = abs(nb-nt);
1069   dv = abs(nr-nl);
1070   int nbh  = Max(nb,nt);
1071   int nbv = Max(nr,nl);
1072   int addh = 0;
1073   int addv = 0;
1074
1075   // orientation of face and 3 main domain for future faces
1076   //       0   top    1
1077   //      1------------1
1078   //       |   |  |   |
1079   //       |   |  |   |
1080   //       | L |  | R |
1081   //  left |   |  |   | rigth
1082   //       |  /    \  |
1083   //       | /  C   \ |
1084   //       |/        \|
1085   //      0------------0
1086   //       0  bottom  1
1087
1088   if(dh>dv) {
1089     addv = (dh-dv)/2;
1090     nbv = nbv + addv;
1091   }
1092   else { // dv>=dh
1093     addh = (dv-dh)/2;
1094     nbh = nbh + addh;
1095   }
1096
1097   Handle(Geom2d_Curve) c2d[4];
1098   for(i=0; i<4; i++) {
1099     c2d[i] = BRep_Tool::CurveOnSurface(quad->edge[i], F,
1100                                        quad->first[i], quad->last[i]);
1101   }
1102
1103   bool loadOk = true;
1104   for(i=0; i<2; i++) {
1105     quad->uv_edges[i] = LoadEdgePoints2(aMesh, F, quad->edge[i], false);
1106     if(!quad->uv_edges[i]) loadOk = false;
1107   }
1108   for(i=2; i<4; i++) {
1109     quad->uv_edges[i] = LoadEdgePoints2(aMesh, F, quad->edge[i], true);
1110     if (!quad->uv_edges[i]) loadOk = false;
1111   }
1112   if (!loadOk) {
1113     INFOS("StdMeshers_Quadrangle_2D::ComputeQuadPref - LoadEdgePoints failed");
1114     QuadDelete( quad );
1115     quad = 0;
1116     return false;
1117   }
1118
1119   UVPtStruct* uv_eb = quad->uv_edges[0];
1120   UVPtStruct* uv_er = quad->uv_edges[1];
1121   UVPtStruct* uv_et = quad->uv_edges[2];
1122   UVPtStruct* uv_el = quad->uv_edges[3];
1123
1124   // arrays for normalized params
1125   //cout<<"Dump B:"<<endl;
1126   TColStd_SequenceOfReal npb, npr, npt, npl;
1127   for(i=0; i<nb; i++) {
1128     npb.Append(uv_eb[i].normParam);
1129     //cout<<"i="<<i<<" par="<<uv_eb[i].param<<" npar="<<uv_eb[i].normParam;
1130     //const SMDS_MeshNode* N = uv_eb[i].node;
1131     //cout<<" node("<<N->X()<<","<<N->Y()<<","<<N->Z()<<")"<<endl;
1132   }
1133   for(i=0; i<nr; i++) {
1134     npr.Append(uv_er[i].normParam);
1135   }
1136   for(i=0; i<nt; i++) {
1137     npt.Append(uv_et[i].normParam);
1138   }
1139   for(i=0; i<nl; i++) {
1140     npl.Append(uv_el[i].normParam);
1141   }
1142
1143   // we have to add few values of params to right and left
1144   // insert them after first param
1145   // insert to right
1146   int dr = nbv - nr;
1147   double dpr = (npr.Value(2) - npr.Value(1))/(dr+1);
1148   for(i=1; i<=dr; i++) {
1149     npr.InsertAfter(1,npr.Value(2)-dpr);
1150   }
1151   // insert to left
1152   int dl = nbv - nl;
1153   dpr = (npl.Value(2) - npl.Value(1))/(dl+1);
1154   for(i=1; i<=dl; i++) {
1155     npl.InsertAfter(1,npl.Value(2)-dpr);
1156   }
1157   //cout<<"npb:";
1158   //for(i=1; i<=npb.Length(); i++) {
1159   //  cout<<" "<<npb.Value(i);
1160   //}
1161   //cout<<endl;
1162
1163   gp_Pnt2d a[4];
1164   c2d[0]->D0(uv_eb[0].param,a[0]);
1165   c2d[0]->D0(uv_eb[nb-1].param,a[1]);
1166   c2d[2]->D0(uv_et[nt-1].param,a[2]);
1167   c2d[2]->D0(uv_et[0].param,a[3]);
1168   //cout<<" a[0]("<<a[0].X()<<","<<a[0].Y()<<")"<<" a[1]("<<a[1].X()<<","<<a[1].Y()<<")"
1169   //    <<" a[2]("<<a[2].X()<<","<<a[2].Y()<<")"<<" a[3]("<<a[3].X()<<","<<a[3].Y()<<")"<<endl;
1170
1171   int nnn = Min(nr,nl);
1172   // auxilary sequence of XY for creation nodes
1173   // in the bottom part of central domain
1174   // it's length must be == nbv-nnn-1
1175   TColgp_SequenceOfXY UVL;
1176   TColgp_SequenceOfXY UVR;
1177
1178   // step1: create faces for left domain
1179   StdMeshers_Array2OfNode NodesL(1,dl+1,1,nl);
1180   // add left nodes
1181   for(j=1; j<=nl; j++)
1182     NodesL.SetValue(1,j,uv_el[j-1].node);
1183   if(dl>0) {
1184     // add top nodes
1185     for(i=1; i<=dl; i++)
1186       NodesL.SetValue(i+1,nl,uv_et[i].node);
1187     // create and add needed nodes
1188     TColgp_SequenceOfXY UVtmp;
1189     for(i=1; i<=dl; i++) {
1190       double x0 = npt.Value(i+1);
1191       double x1 = x0;
1192       // diagonal node
1193       double y0 = npl.Value(i+1);
1194       double y1 = npr.Value(i+1);
1195       gp_XY UV = CalcUV(x0, x1, y0, y1, quad, a[0], a[1], a[2], a[3],
1196                         c2d[0], c2d[1], c2d[2], c2d[3]);
1197       gp_Pnt P = S->Value(UV.X(),UV.Y());
1198       SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1199       meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1200       NodesL.SetValue(i+1,1,N);
1201       if(UVL.Length()<nbv-nnn-1) UVL.Append(UV);
1202       // internal nodes
1203       for(j=2; j<nl; j++) {
1204         double y0 = npl.Value(dl+j);
1205         double y1 = npr.Value(dl+j);
1206         gp_XY UV = CalcUV(x0, x1, y0, y1, quad, a[0], a[1], a[2], a[3],
1207                           c2d[0], c2d[1], c2d[2], c2d[3]);
1208         gp_Pnt P = S->Value(UV.X(),UV.Y());
1209         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1210         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1211         NodesL.SetValue(i+1,j,N);
1212         if( i==dl ) UVtmp.Append(UV);
1213       }
1214     }
1215     for(i=1; i<=UVtmp.Length() && UVL.Length()<nbv-nnn-1; i++) {
1216       UVL.Append(UVtmp.Value(i));
1217     }
1218     //cout<<"Dump NodesL:"<<endl;
1219     //for(i=1; i<=dl+1; i++) {
1220     //  cout<<"i="<<i;
1221     //  for(j=1; j<=nl; j++) {
1222     //    cout<<" ("<<NodesL.Value(i,j)->X()<<","<<NodesL.Value(i,j)->Y()<<","<<NodesL.Value(i,j)->Z()<<")";
1223     //  }
1224     //  cout<<endl;
1225     //}
1226     // create faces
1227     for(i=1; i<=dl; i++) {
1228       for(j=1; j<nl; j++) {
1229         if(WisF) {
1230           SMDS_MeshFace* F =
1231             myTool->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j),
1232                             NodesL.Value(i+1,j+1), NodesL.Value(i,j+1));
1233           meshDS->SetMeshElementOnShape(F, geomFaceID);
1234         }
1235         else {
1236           SMDS_MeshFace* F =
1237             myTool->AddFace(NodesL.Value(i,j), NodesL.Value(i,j+1),
1238                             NodesL.Value(i+1,j+1), NodesL.Value(i+1,j));
1239           meshDS->SetMeshElementOnShape(F, geomFaceID);
1240         }
1241       }
1242     }
1243   }
1244   else {
1245     // fill UVL using c2d
1246     for(i=1; i<npl.Length() && UVL.Length()<nbv-nnn-1; i++) {
1247       gp_Pnt2d p2d;
1248       c2d[3]->D0(uv_el[i].param,p2d);
1249       UVL.Append(p2d.XY());
1250     }
1251   }
1252
1253   // step2: create faces for right domain
1254   StdMeshers_Array2OfNode NodesR(1,dr+1,1,nr);
1255   // add right nodes
1256   for(j=1; j<=nr; j++)
1257     NodesR.SetValue(1,j,uv_er[nr-j].node);
1258   if(dr>0) {
1259     // add top nodes
1260     for(i=1; i<=dr; i++)
1261       NodesR.SetValue(i+1,1,uv_et[nt-1-i].node);
1262     // create and add needed nodes
1263     TColgp_SequenceOfXY UVtmp;
1264     for(i=1; i<=dr; i++) {
1265       double x0 = npt.Value(nt-i);
1266       double x1 = x0;
1267       // diagonal node
1268       double y0 = npl.Value(i+1);
1269       double y1 = npr.Value(i+1);
1270       gp_XY UV = CalcUV(x0, x1, y0, y1, quad, a[0], a[1], a[2], a[3],
1271                         c2d[0], c2d[1], c2d[2], c2d[3]);
1272       gp_Pnt P = S->Value(UV.X(),UV.Y());
1273       SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1274       meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1275       NodesR.SetValue(i+1,nr,N);
1276       if(UVR.Length()<nbv-nnn-1) UVR.Append(UV);
1277       // internal nodes
1278       for(j=2; j<nr; j++) {
1279         double y0 = npl.Value(nbv-j+1);
1280         double y1 = npr.Value(nbv-j+1);
1281         gp_XY UV = CalcUV(x0, x1, y0, y1, quad, a[0], a[1], a[2], a[3],
1282                           c2d[0], c2d[1], c2d[2], c2d[3]);
1283         gp_Pnt P = S->Value(UV.X(),UV.Y());
1284         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1285         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1286         NodesR.SetValue(i+1,j,N);
1287         if( i==dr ) UVtmp.Prepend(UV);
1288       }
1289     }
1290     for(i=1; i<=UVtmp.Length() && UVR.Length()<nbv-nnn-1; i++) {
1291       UVR.Append(UVtmp.Value(i));
1292     }
1293     // create faces
1294     for(i=1; i<=dr; i++) {
1295       for(j=1; j<nr; j++) {
1296         if(WisF) {
1297           SMDS_MeshFace* F =
1298             myTool->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j),
1299                             NodesR.Value(i+1,j+1), NodesR.Value(i,j+1));
1300           meshDS->SetMeshElementOnShape(F, geomFaceID);
1301         }
1302         else {
1303           SMDS_MeshFace* F =
1304             myTool->AddFace(NodesR.Value(i,j), NodesR.Value(i,j+1),
1305                             NodesR.Value(i+1,j+1), NodesR.Value(i+1,j));
1306           meshDS->SetMeshElementOnShape(F, geomFaceID);
1307         }
1308       }
1309     }
1310   }
1311   else {
1312     // fill UVR using c2d
1313     for(i=1; i<npr.Length() && UVR.Length()<nbv-nnn-1; i++) {
1314       gp_Pnt2d p2d;
1315       c2d[1]->D0(uv_er[i].param,p2d);
1316       UVR.Append(p2d.XY());
1317     }
1318   }
1319
1320   // step3: create faces for central domain
1321   StdMeshers_Array2OfNode NodesC(1,nb,1,nbv);
1322   // add first string using NodesL
1323   for(i=1; i<=dl+1; i++)
1324     NodesC.SetValue(1,i,NodesL(i,1));
1325   for(i=2; i<=nl; i++)
1326     NodesC.SetValue(1,dl+i,NodesL(dl+1,i));
1327   // add last string using NodesR
1328   for(i=1; i<=dr+1; i++)
1329     NodesC.SetValue(nb,i,NodesR(i,nr));
1330   for(i=1; i<nr; i++)
1331     NodesC.SetValue(nb,dr+i+1,NodesR(dr+1,nr-i));
1332   // add top nodes (last columns)
1333   for(i=dl+2; i<nbh-dr; i++)
1334     NodesC.SetValue(i-dl,nbv,uv_et[i-1].node);
1335   // add bottom nodes (first columns)
1336   for(i=2; i<nb; i++) {
1337     NodesC.SetValue(i,1,uv_eb[i-1].node);
1338     gp_Pnt2d p2d;
1339     c2d[0]->D0(uv_eb[i-1].param,p2d);
1340   }
1341   // create and add needed nodes
1342   // add linear layers
1343   for(i=2; i<nb; i++) {
1344     double x0 = npt.Value(dl+i);
1345     double x1 = x0;
1346     for(j=1; j<nnn; j++) {
1347       double y0 = npl.Value(nbv-nnn+j);
1348       double y1 = npr.Value(nbv-nnn+j);
1349       gp_XY UV = CalcUV(x0, x1, y0, y1, quad, a[0], a[1], a[2], a[3],
1350                         c2d[0], c2d[1], c2d[2], c2d[3]);
1351       gp_Pnt P = S->Value(UV.X(),UV.Y());
1352       SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1353       meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1354       NodesC.SetValue(i,nbv-nnn+j,N);
1355     }
1356   }
1357   // add diagonal layers
1358   //cout<<"UVL.Length()="<<UVL.Length()<<" UVR.Length()="<<UVR.Length()<<endl;
1359   //cout<<"Dump UVL:"<<endl;
1360   //for(i=1; i<=UVL.Length(); i++) {
1361   //  cout<<" ("<<UVL.Value(i).X()<<","<<UVL.Value(i).Y()<<")";
1362   //}
1363   //cout<<endl;
1364   for(i=1; i<nbv-nnn; i++) {
1365     double du = UVR.Value(i).X() - UVL.Value(i).X();
1366     double dv = UVR.Value(i).Y() - UVL.Value(i).Y();
1367     for(j=2; j<nb; j++) {
1368       double u = UVL.Value(i).X() + du*npb.Value(j);
1369       double v = UVL.Value(i).Y() + dv*npb.Value(j);
1370       gp_Pnt P = S->Value(u,v);
1371       SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1372       meshDS->SetNodeOnFace(N, geomFaceID, u, v);
1373       NodesC.SetValue(j,i+1,N);
1374     }
1375   }
1376   // create faces
1377   for(i=1; i<nb; i++) {
1378     for(j=1; j<nbv; j++) {
1379       if(WisF) {
1380         SMDS_MeshFace* F =
1381           myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
1382                           NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
1383         meshDS->SetMeshElementOnShape(F, geomFaceID);
1384       }
1385       else {
1386         SMDS_MeshFace* F =
1387           myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1),
1388                           NodesC.Value(i+1,j+1), NodesC.Value(i+1,j));
1389         meshDS->SetMeshElementOnShape(F, geomFaceID);
1390       }
1391     }
1392   }
1393
1394   QuadDelete(quad);
1395   bool isOk = true;
1396   return isOk;
1397 }
1398
1399
1400 //=============================================================================
1401 /*!
1402  *  LoadEdgePoints2
1403  */
1404 //=============================================================================
1405 UVPtStruct* StdMeshers_Quadrangle_2D::LoadEdgePoints2 (SMESH_Mesh & aMesh,
1406                                                        const TopoDS_Face& F,
1407                                                        const TopoDS_Edge& E,
1408                                                        bool IsReverse)
1409 {
1410   //MESSAGE("StdMeshers_Quadrangle_2D::LoadEdgePoints");
1411   // --- IDNodes of first and last Vertex
1412   TopoDS_Vertex VFirst, VLast;
1413   TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
1414
1415   ASSERT(!VFirst.IsNull());
1416   SMDS_NodeIteratorPtr lid = aMesh.GetSubMesh(VFirst)->GetSubMeshDS()->GetNodes();
1417   if (!lid->more()) {
1418     MESSAGE ( "NO NODE BUILT ON VERTEX" );
1419     return 0;
1420   }
1421   const SMDS_MeshNode* idFirst = lid->next();
1422
1423   ASSERT(!VLast.IsNull());
1424   lid = aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes();
1425   if (!lid->more()) {
1426     MESSAGE ( "NO NODE BUILT ON VERTEX" );
1427     return 0;
1428   }
1429   const SMDS_MeshNode* idLast = lid->next();
1430
1431   // --- edge internal IDNodes (relies on good order storage, not checked)
1432
1433   map<double, const SMDS_MeshNode *> params;
1434   SMDS_NodeIteratorPtr ite = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes();
1435   int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
1436
1437   if(!_quadraticMesh) {
1438     while(ite->more()) {
1439       const SMDS_MeshNode* node = ite->next();
1440       const SMDS_EdgePosition* epos =
1441         static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
1442       double param = epos->GetUParameter();
1443       params[param] = node;
1444     }
1445   }
1446   else {
1447     vector<const SMDS_MeshNode*> nodes(nbPoints+2);
1448     nodes[0] = idFirst;
1449     nodes[nbPoints+1] = idLast;
1450     nbPoints = nbPoints/2;
1451     int nn = 1;
1452     while(ite->more()) {
1453       const SMDS_MeshNode* node = ite->next();
1454       nodes[nn++] = node;
1455       // check if node is medium
1456       bool IsMedium = false;
1457       SMDS_ElemIteratorPtr itn = node->GetInverseElementIterator();
1458       while (itn->more()) {
1459         const SMDS_MeshElement* elem = itn->next();
1460         if ( elem->GetType() != SMDSAbs_Edge )
1461           continue;
1462         if(elem->IsMediumNode(node)) {
1463           IsMedium = true;
1464           break;
1465         }
1466       }
1467       if(IsMedium)
1468         continue;
1469       const SMDS_EdgePosition* epos =
1470         static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
1471       double param = epos->GetUParameter();
1472       params[param] = node;
1473     }
1474   }
1475
1476   if (nbPoints != params.size()) {
1477     MESSAGE( "BAD NODE ON EDGE POSITIONS" );
1478     return 0;
1479   }
1480   UVPtStruct* uvslf = new UVPtStruct[nbPoints + 2];
1481
1482   double f, l;
1483   Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
1484
1485   const TopoDS_Wire& W = BRepTools::OuterWire(F);
1486   bool FisF = (F.Orientation()==TopAbs_FORWARD);
1487   bool WisF = (W.Orientation()==TopAbs_FORWARD);
1488   bool isForward = (E.Orientation()==TopAbs_FORWARD);
1489   //if(isForward) cout<<"E is FORWARD"<<endl;
1490   //else cout<<"E is REVERSED"<<endl;
1491   if(!WisF) isForward = !isForward;
1492   if(!FisF) isForward = !isForward;
1493   //bool isForward = !(E.Orientation()==TopAbs_FORWARD);
1494   if(IsReverse) isForward = !isForward;
1495   double paramin = 0;
1496   double paramax = 0;
1497   if (isForward) {
1498     paramin = f;
1499     paramax = l;
1500     gp_Pnt2d p = C2d->Value(f); // first point = Vertex Forward
1501     uvslf[0].x = p.X();
1502     uvslf[0].y = p.Y();
1503     uvslf[0].param = f;
1504     uvslf[0].node = idFirst;
1505     //MESSAGE("__ f "<<f<<" "<<uvslf[0].x <<" "<<uvslf[0].y);
1506     map < double, const SMDS_MeshNode* >::iterator itp = params.begin();
1507     for (int i = 1; i <= nbPoints; i++) { // nbPoints internal
1508       double param = (*itp).first;
1509       gp_Pnt2d p = C2d->Value(param);
1510       uvslf[i].x = p.X();
1511       uvslf[i].y = p.Y();
1512       uvslf[i].param = param;
1513       uvslf[i].node = (*itp).second;
1514       //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[i].x <<" "<<uvslf[i].y);
1515       itp++;
1516     }
1517     p = C2d->Value(l);          // last point = Vertex Reversed
1518     uvslf[nbPoints + 1].x = p.X();
1519     uvslf[nbPoints + 1].y = p.Y();
1520     uvslf[nbPoints + 1].param = l;
1521     uvslf[nbPoints + 1].node = idLast;
1522     //MESSAGE("__ l "<<l<<" "<<uvslf[nbPoints+1].x <<" "<<uvslf[nbPoints+1].y);
1523   }
1524   else {
1525     paramin = l;
1526     paramax = f;
1527     gp_Pnt2d p = C2d->Value(l); // first point = Vertex Reversed
1528     uvslf[0].x = p.X();
1529     uvslf[0].y = p.Y();
1530     uvslf[0].param = l;
1531     uvslf[0].node = idLast;
1532     //MESSAGE("__ l "<<l<<" "<<uvslf[0].x <<" "<<uvslf[0].y);
1533     map < double, const SMDS_MeshNode* >::reverse_iterator itp = params.rbegin();
1534     for (int j = nbPoints; j >= 1; j--) { // nbPoints internal
1535       double param = (*itp).first;
1536       int i = nbPoints + 1 - j;
1537       gp_Pnt2d p = C2d->Value(param);
1538       uvslf[i].x = p.X();
1539       uvslf[i].y = p.Y();
1540       uvslf[i].param = param;
1541       uvslf[i].node = (*itp).second;
1542       //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[i].x <<" "<<uvslf[i].y);
1543       itp++;
1544     }
1545     p = C2d->Value(f);          // last point = Vertex Forward
1546     uvslf[nbPoints + 1].x = p.X();
1547     uvslf[nbPoints + 1].y = p.Y();
1548     uvslf[nbPoints + 1].param = f;
1549     uvslf[nbPoints + 1].node = idFirst;
1550     //MESSAGE("__ f "<<f<<" "<<uvslf[nbPoints+1].x <<" "<<uvslf[nbPoints+1].y);
1551   }
1552
1553   ASSERT(paramin != paramax);
1554   for (int i = 0; i < nbPoints + 2; i++) {
1555     uvslf[i].normParam = (uvslf[i].param - paramin) / (paramax - paramin);
1556   }
1557
1558   return uvslf;
1559 }
1560
1561
1562 //=============================================================================
1563 /*!
1564  *  LoadEdgePoints
1565  */
1566 //=============================================================================
1567 UVPtStruct* StdMeshers_Quadrangle_2D::LoadEdgePoints (SMESH_Mesh & aMesh,
1568                                                       const TopoDS_Face& F,
1569                                                       const TopoDS_Edge& E,
1570                                                       double first, double last)
1571 //                        bool isForward)
1572 {
1573   //MESSAGE("StdMeshers_Quadrangle_2D::LoadEdgePoints");
1574
1575   // --- IDNodes of first and last Vertex
1576
1577   TopoDS_Vertex VFirst, VLast;
1578   TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
1579
1580   ASSERT(!VFirst.IsNull());
1581   SMDS_NodeIteratorPtr lid = aMesh.GetSubMesh(VFirst)->GetSubMeshDS()->GetNodes();
1582   if (!lid->more())
1583   {
1584     MESSAGE ( "NO NODE BUILT ON VERTEX" );
1585     return 0;
1586   }
1587   const SMDS_MeshNode* idFirst = lid->next();
1588
1589   ASSERT(!VLast.IsNull());
1590   lid = aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes();
1591   if (!lid->more())
1592   {
1593     MESSAGE ( "NO NODE BUILT ON VERTEX" );
1594     return 0;
1595   }
1596   const SMDS_MeshNode* idLast = lid->next();
1597
1598   // --- edge internal IDNodes (relies on good order storage, not checked)
1599
1600 //  if(_quadraticMesh) {
1601     // fill myNLinkNodeMap
1602 //    SMDS_ElemIteratorPtr iter = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetElements();
1603 //    while(iter->more()) {
1604 //      const SMDS_MeshElement* elem = iter->next();
1605 //      SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
1606 //      const SMDS_MeshNode* n1 = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
1607 //      const SMDS_MeshNode* n2 = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
1608 //      const SMDS_MeshNode* n3 = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
1609 //      NLink link(( n1 < n2 ? n1 : n2 ), ( n1 < n2 ? n2 : n1 ));
1610 //      myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n3));
1611 //      myNLinkNodeMap[link] = n3;
1612 //    }
1613 //  }
1614
1615   map<double, const SMDS_MeshNode *> params;
1616   SMDS_NodeIteratorPtr ite = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes();
1617   int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
1618
1619   if(!_quadraticMesh) {
1620     while(ite->more()) {
1621       const SMDS_MeshNode* node = ite->next();
1622       const SMDS_EdgePosition* epos =
1623         static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
1624       double param = epos->GetUParameter();
1625       params[param] = node;
1626     }
1627   }
1628   else {
1629     nbPoints = nbPoints/2;
1630     while(ite->more()) {
1631       const SMDS_MeshNode* node = ite->next();
1632       // check if node is medium
1633       bool IsMedium = false;
1634       SMDS_ElemIteratorPtr itn = node->GetInverseElementIterator();
1635       while (itn->more()) {
1636         const SMDS_MeshElement* elem = itn->next();
1637         if ( elem->GetType() != SMDSAbs_Edge )
1638           continue;
1639         if(elem->IsMediumNode(node)) {
1640           IsMedium = true;
1641           break;
1642         }
1643       }
1644       if(IsMedium)
1645         continue;
1646       const SMDS_EdgePosition* epos =
1647         static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
1648       double param = epos->GetUParameter();
1649       params[param] = node;
1650     }
1651   }
1652
1653   if (nbPoints != params.size()) {
1654     MESSAGE( "BAD NODE ON EDGE POSITIONS" );
1655     return 0;
1656   }
1657   UVPtStruct* uvslf = new UVPtStruct[nbPoints + 2];
1658
1659   double f, l;
1660   Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
1661
1662   bool isForward = (((l - f) * (last - first)) > 0);
1663   double paramin = 0;
1664   double paramax = 0;
1665   if (isForward)
1666   {
1667     paramin = f;
1668     paramax = l;
1669     gp_Pnt2d p = C2d->Value(f); // first point = Vertex Forward
1670     uvslf[0].x = p.X();
1671     uvslf[0].y = p.Y();
1672     uvslf[0].param = f;
1673     uvslf[0].node = idFirst;
1674     //MESSAGE("__ f "<<f<<" "<<uvslf[0].x <<" "<<uvslf[0].y);
1675     map < double, const SMDS_MeshNode* >::iterator itp = params.begin();
1676     for (int i = 1; i <= nbPoints; i++) // nbPoints internal
1677     {
1678       double param = (*itp).first;
1679       gp_Pnt2d p = C2d->Value(param);
1680       uvslf[i].x = p.X();
1681       uvslf[i].y = p.Y();
1682       uvslf[i].param = param;
1683       uvslf[i].node = (*itp).second;
1684       //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[i].x <<" "<<uvslf[i].y);
1685       itp++;
1686     }
1687     p = C2d->Value(l);          // last point = Vertex Reversed
1688     uvslf[nbPoints + 1].x = p.X();
1689     uvslf[nbPoints + 1].y = p.Y();
1690     uvslf[nbPoints + 1].param = l;
1691     uvslf[nbPoints + 1].node = idLast;
1692     //MESSAGE("__ l "<<l<<" "<<uvslf[nbPoints+1].x <<" "<<uvslf[nbPoints+1].y);
1693   } else
1694   {
1695     paramin = l;
1696     paramax = f;
1697     gp_Pnt2d p = C2d->Value(l); // first point = Vertex Reversed
1698     uvslf[0].x = p.X();
1699     uvslf[0].y = p.Y();
1700     uvslf[0].param = l;
1701     uvslf[0].node = idLast;
1702     //MESSAGE("__ l "<<l<<" "<<uvslf[0].x <<" "<<uvslf[0].y);
1703     map < double, const SMDS_MeshNode* >::reverse_iterator itp = params.rbegin();
1704
1705     for (int j = nbPoints; j >= 1; j--) // nbPoints internal
1706     {
1707       double param = (*itp).first;
1708       int i = nbPoints + 1 - j;
1709       gp_Pnt2d p = C2d->Value(param);
1710       uvslf[i].x = p.X();
1711       uvslf[i].y = p.Y();
1712       uvslf[i].param = param;
1713       uvslf[i].node = (*itp).second;
1714       //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[i].x <<" "<<uvslf[i].y);
1715       itp++;
1716     }
1717     p = C2d->Value(f);          // last point = Vertex Forward
1718     uvslf[nbPoints + 1].x = p.X();
1719     uvslf[nbPoints + 1].y = p.Y();
1720     uvslf[nbPoints + 1].param = f;
1721     uvslf[nbPoints + 1].node = idFirst;
1722     //MESSAGE("__ f "<<f<<" "<<uvslf[nbPoints+1].x <<" "<<uvslf[nbPoints+1].y);
1723   }
1724
1725   ASSERT(paramin != paramax);
1726   for (int i = 0; i < nbPoints + 2; i++)
1727   {
1728     uvslf[i].normParam = (uvslf[i].param - paramin) / (paramax - paramin);
1729   }
1730
1731   return uvslf;
1732 }
1733
1734 //=============================================================================
1735 /*!
1736  *  MakeEdgePoints
1737  */
1738 //=============================================================================
1739 UVPtStruct* StdMeshers_Quadrangle_2D::MakeEdgePoints (SMESH_Mesh & aMesh,
1740                                                       const TopoDS_Face& F,
1741                                                       const TopoDS_Edge& E,
1742                                                       double first, double last,
1743                                                       int nb_segm)
1744 {
1745 //  MESSAGE("StdMeshers_Quadrangle_2D::MakeEdgePoints");
1746
1747   UVPtStruct* uvslf = new UVPtStruct[nb_segm + 1];
1748   list<double> params;
1749
1750   // --- edge internal points
1751   double fi, li;
1752   Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, fi, li);
1753   if (!Curve.IsNull()) {
1754     try {
1755       GeomAdaptor_Curve C3d (Curve);
1756       double length = EdgeLength(E);
1757       double eltSize = length / nb_segm;
1758       GCPnts_UniformAbscissa Discret (C3d, eltSize, fi, li);
1759       if (!Discret.IsDone()) return false;
1760       int NbPoints = Discret.NbPoints();
1761       for (int i = 1; i <= NbPoints; i++) {
1762         double param = Discret.Parameter(i);
1763         params.push_back(param);
1764       }
1765     }
1766     catch (Standard_Failure) {
1767       return 0;
1768     }
1769   }
1770   else
1771   {
1772     // Edge is a degenerated Edge
1773     BRep_Tool::Range(E, fi, li);
1774     double du = (li - fi) / nb_segm;
1775     for (int i = 1; i <= nb_segm + 1; i++)
1776     {
1777       double param = fi + (i - 1) * du;
1778       params.push_back(param);
1779     }
1780   }
1781
1782   double f, l;
1783   Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
1784   ASSERT(f != l);
1785
1786   bool isForward = (((l - f) * (last - first)) > 0);
1787   if (isForward) {
1788     list<double>::iterator itU = params.begin();
1789     for (int i = 0; i <= nb_segm; i++) // nbPoints internal
1790     {
1791       double param = *itU;
1792       gp_Pnt2d p = C2d->Value(param);
1793       uvslf[i].x = p.X();
1794       uvslf[i].y = p.Y();
1795       uvslf[i].param = param;
1796       uvslf[i].normParam = (param - f) / (l - f);
1797       itU++;
1798     }
1799   } else {
1800     list<double>::reverse_iterator itU = params.rbegin();
1801     for (int j = nb_segm; j >= 0; j--) // nbPoints internal
1802     {
1803       double param = *itU;
1804       int i = nb_segm - j;
1805       gp_Pnt2d p = C2d->Value(param);
1806       uvslf[i].x = p.X();
1807       uvslf[i].y = p.Y();
1808       uvslf[i].param = param;
1809       uvslf[i].normParam = (param - l) / (f - l);
1810       itU++;
1811     }
1812   }
1813
1814   return uvslf;
1815 }
1816
1817
1818 //=============================================================================
1819 /*!
1820  *
1821  */
1822 //=============================================================================
1823
1824 ostream & StdMeshers_Quadrangle_2D::SaveTo(ostream & save)
1825 {
1826   return save;
1827 }
1828
1829 //=============================================================================
1830 /*!
1831  *
1832  */
1833 //=============================================================================
1834
1835 istream & StdMeshers_Quadrangle_2D::LoadFrom(istream & load)
1836 {
1837   return load;
1838 }
1839
1840 //=============================================================================
1841 /*!
1842  *
1843  */
1844 //=============================================================================
1845
1846 ostream & operator <<(ostream & save, StdMeshers_Quadrangle_2D & hyp)
1847 {
1848   return hyp.SaveTo( save );
1849 }
1850
1851 //=============================================================================
1852 /*!
1853  *
1854  */
1855 //=============================================================================
1856
1857 istream & operator >>(istream & load, StdMeshers_Quadrangle_2D & hyp)
1858 {
1859   return hyp.LoadFrom( load );
1860 }