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