src/StdMeshers/StdMeshers_Quadrangle_2D.cxx

   1 //  Copyright (C) 2007-2010  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 //  Copyright (C) 2003-2007  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.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21
  22 //  File   : StdMeshers_Quadrangle_2D.cxx
  23 //  Author : Paul RASCLE, EDF
  24 //  Module : SMESH
  25
  26 #include "StdMeshers_Quadrangle_2D.hxx"
  27
  28 #include "StdMeshers_FaceSide.hxx"
  29
  30 #include "StdMeshers_QuadrangleParams.hxx"
  31
  32 #include "SMESH_Gen.hxx"
  33 #include "SMESH_Mesh.hxx"
  34 #include "SMESH_subMesh.hxx"
  35 #include "SMESH_MesherHelper.hxx"
  36 #include "SMESH_Block.hxx"
  37 #include "SMESH_Comment.hxx"
  38
  39 #include "SMDS_MeshElement.hxx"
  40 #include "SMDS_MeshNode.hxx"
  41 #include "SMDS_EdgePosition.hxx"
  42 #include "SMDS_FacePosition.hxx"
  43
  44 #include <BRep_Tool.hxx>
  45 #include <Geom_Surface.hxx>
  46 #include <NCollection_DefineArray2.hxx>
  47 #include <Precision.hxx>
  48 #include <TColStd_SequenceOfReal.hxx>
  49 #include <TColStd_SequenceOfInteger.hxx>
  50 #include <TColgp_SequenceOfXY.hxx>
  51 #include <TopExp.hxx>
  52 #include <TopExp_Explorer.hxx>
  53 #include <TopTools_ListIteratorOfListOfShape.hxx>
  54 #include <TopTools_MapOfShape.hxx>
  55 #include <TopoDS.hxx>
  56
  57 #include "utilities.h"
  58 #include "Utils_ExceptHandlers.hxx"
  59
  60 #ifndef StdMeshers_Array2OfNode_HeaderFile
  61 #define StdMeshers_Array2OfNode_HeaderFile
  62 typedef const SMDS_MeshNode* SMDS_MeshNodePtr;
  63 DEFINE_BASECOLLECTION (StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr)
  64 DEFINE_ARRAY2(StdMeshers_Array2OfNode,
  65               StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr)
  66 #endif
  67
  68 using namespace std;
  69
  70 typedef gp_XY gp_UV;
  71 typedef SMESH_Comment TComm;
  72
  73 //=============================================================================
  74 /*!
  75  *
  76  */
  77 //=============================================================================
  78
  79 StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId,
  80                                                     SMESH_Gen* gen)
  81      : SMESH_2D_Algo(hypId, studyId, gen)
  82 {
  83   MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D");
  84   _name = "Quadrangle_2D";
  85   _shapeType = (1 << TopAbs_FACE);
  86   _compatibleHypothesis.push_back("QuadrangleParams");
  87   _compatibleHypothesis.push_back("QuadranglePreference");
  88   _compatibleHypothesis.push_back("TrianglePreference");
  89   myHelper = 0;
  90 }
  91
  92 //=============================================================================
  93 /*!
  94  *
  95  */
  96 //=============================================================================
  97
  98 StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D()
  99 {
 100   MESSAGE("StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D");
 101 }
 102
 103 //=============================================================================
 104 /*!
 105  *
 106  */
 107 //=============================================================================
 108
 109 bool StdMeshers_Quadrangle_2D::CheckHypothesis
 110                          (SMESH_Mesh&                          aMesh,
 111                           const TopoDS_Shape&                  aShape,
 112                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 113 {
 114   bool isOk = true;
 115   aStatus = SMESH_Hypothesis::HYP_OK;
 116
 117   const list <const SMESHDS_Hypothesis * >& hyps =
 118     GetUsedHypothesis(aMesh, aShape, false);
 119   const SMESHDS_Hypothesis * aHyp = 0;
 120
 121   myTriaVertexID = -1;
 122   myQuadType = QUAD_STANDARD;
 123   myQuadranglePreference = false;
 124   myTrianglePreference = false;
 125
 126   bool isFirstParams = true;
 127
 128   // First assigned hypothesis (if any) is processed now
 129   if (hyps.size() > 0) {
 130     aHyp = hyps.front();
 131     if (strcmp("QuadrangleParams", aHyp->GetName()) == 0) {
 132       const StdMeshers_QuadrangleParams* aHyp1 =
 133         (const StdMeshers_QuadrangleParams*)aHyp;
 134       myTriaVertexID = aHyp1->GetTriaVertex();
 135       myQuadType = aHyp1->GetQuadType();
 136       if (myQuadType == QUAD_QUADRANGLE_PREF ||
 137           myQuadType == QUAD_QUADRANGLE_PREF_REVERSED)
 138         myQuadranglePreference = true;
 139       else if (myQuadType == QUAD_TRIANGLE_PREF)
 140         myTrianglePreference = true;
 141     }
 142     else if (strcmp("QuadranglePreference", aHyp->GetName()) == 0) {
 143       isFirstParams = false;
 144       myQuadranglePreference = true;
 145     }
 146     else if (strcmp("TrianglePreference", aHyp->GetName()) == 0){
 147       isFirstParams = false;
 148       myTrianglePreference = true;
 149     }
 150     else {
 151       isFirstParams = false;
 152     }
 153   }
 154
 155   // Second(last) assigned hypothesis (if any) is processed now
 156   if (hyps.size() > 1) {
 157     aHyp = hyps.back();
 158     if (isFirstParams) {
 159       if (strcmp("QuadranglePreference", aHyp->GetName()) == 0) {
 160         myQuadranglePreference = true;
 161         myTrianglePreference = false;
 162         myQuadType = QUAD_STANDARD;
 163       }
 164       else if (strcmp("TrianglePreference", aHyp->GetName()) == 0){
 165         myQuadranglePreference = false;
 166         myTrianglePreference = true;
 167         myQuadType = QUAD_STANDARD;
 168       }
 169     }
 170     else {
 171       const StdMeshers_QuadrangleParams* aHyp2 =
 172         (const StdMeshers_QuadrangleParams*)aHyp;
 173       myTriaVertexID = aHyp2->GetTriaVertex();
 174
 175       if (!myQuadranglePreference && !myTrianglePreference) { // priority of hypos
 176         myQuadType = aHyp2->GetQuadType();
 177         if (myQuadType == QUAD_QUADRANGLE_PREF ||
 178             myQuadType == QUAD_QUADRANGLE_PREF_REVERSED)
 179           myQuadranglePreference = true;
 180         else if (myQuadType == QUAD_TRIANGLE_PREF)
 181           myTrianglePreference = true;
 182       }
 183     }
 184   }
 185
 186   return isOk;
 187 }
 188
 189 //=============================================================================
 190 /*!
 191  *
 192  */
 193 //=============================================================================
 194
 195 bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh,
 196                                         const TopoDS_Shape& aShape)// throw (SALOME_Exception)
 197 {
 198   // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside
 199   //Unexpect aCatchSalomeException);
 200
 201   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 202   aMesh.GetSubMesh(aShape);
 203
 204   SMESH_MesherHelper helper (aMesh);
 205   myHelper = &helper;
 206
 207   _quadraticMesh = myHelper->IsQuadraticSubMesh(aShape);
 208   myNeedSmooth = false;
 209
 210   FaceQuadStruct *quad = CheckNbEdges(aMesh, aShape);
 211   std::auto_ptr<FaceQuadStruct> quadDeleter (quad); // to delete quad at exit from Compute()
 212   if (!quad)
 213     return false;
 214
 215   if (myQuadranglePreference) {
 216     int n1 = quad->side[0]->NbPoints();
 217     int n2 = quad->side[1]->NbPoints();
 218     int n3 = quad->side[2]->NbPoints();
 219     int n4 = quad->side[3]->NbPoints();
 220     int nfull = n1+n2+n3+n4;
 221     int ntmp = nfull/2;
 222     ntmp = ntmp*2;
 223     if (nfull == ntmp && ((n1 != n3) || (n2 != n4))) {
 224       // special path for using only quandrangle faces
 225       bool ok = ComputeQuadPref(aMesh, aShape, quad);
 226       if ( ok && myNeedSmooth )
 227         Smooth( quad );
 228       return ok;
 229     }
 230   }
 231   else if (myQuadType == QUAD_REDUCED) {
 232     int n1 = quad->side[0]->NbPoints();
 233     int n2 = quad->side[1]->NbPoints();
 234     int n3 = quad->side[2]->NbPoints();
 235     int n4 = quad->side[3]->NbPoints();
 236     int n13 = n1 - n3;
 237     int n24 = n2 - n4;
 238     int n13tmp = n13/2; n13tmp = n13tmp*2;
 239     int n24tmp = n24/2; n24tmp = n24tmp*2;
 240     if ((n1 == n3 && n2 != n4 && n24tmp == n24) ||
 241         (n2 == n4 && n1 != n3 && n13tmp == n13)) {
 242       bool ok = ComputeReduced(aMesh, aShape, quad);
 243       if ( ok && myNeedSmooth )
 244         Smooth( quad );
 245       return ok;
 246     }
 247   }
 248
 249   // set normalized grid on unit square in parametric domain
 250
 251   if (!SetNormalizedGrid(aMesh, aShape, quad))
 252     return false;
 253
 254   // --- compute 3D values on points, store points & quadrangles
 255
 256   int nbdown  = quad->side[0]->NbPoints();
 257   int nbup    = quad->side[2]->NbPoints();
 258
 259   int nbright = quad->side[1]->NbPoints();
 260   int nbleft  = quad->side[3]->NbPoints();
 261
 262   int nbhoriz  = Min(nbdown, nbup);
 263   int nbvertic = Min(nbright, nbleft);
 264
 265   const TopoDS_Face& F = TopoDS::Face(aShape);
 266   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
 267
 268   // internal mesh nodes
 269   int i, j, geomFaceID = meshDS->ShapeToIndex(F);
 270   for (i = 1; i < nbhoriz - 1; i++) {
 271     for (j = 1; j < nbvertic - 1; j++) {
 272       int ij = j * nbhoriz + i;
 273       double u = quad->uv_grid[ij].u;
 274       double v = quad->uv_grid[ij].v;
 275       gp_Pnt P = S->Value(u, v);
 276       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 277       meshDS->SetNodeOnFace(node, geomFaceID, u, v);
 278       quad->uv_grid[ij].node = node;
 279     }
 280   }
 281
 282   // mesh faces
 283
 284   //             [2]
 285   //      --.--.--.--.--.--  nbvertic
 286   //     |                 | ^
 287   //     |                 | ^
 288   // [3] |                 | ^ j  [1]
 289   //     |                 | ^
 290   //     |                 | ^
 291   //      ---.----.----.---  0
 292   //     0 > > > > > > > > nbhoriz
 293   //              i
 294   //             [0]
 295
 296   i = 0;
 297   int ilow = 0;
 298   int iup = nbhoriz - 1;
 299   if (quad->isEdgeOut[3]) { ilow++; } else { if (quad->isEdgeOut[1]) iup--; }
 300
 301   int jlow = 0;
 302   int jup = nbvertic - 1;
 303   if (quad->isEdgeOut[0]) { jlow++; } else { if (quad->isEdgeOut[2]) jup--; }
 304
 305   // regular quadrangles
 306   for (i = ilow; i < iup; i++) {
 307     for (j = jlow; j < jup; j++) {
 308       const SMDS_MeshNode *a, *b, *c, *d;
 309       a = quad->uv_grid[j * nbhoriz + i].node;
 310       b = quad->uv_grid[j * nbhoriz + i + 1].node;
 311       c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node;
 312       d = quad->uv_grid[(j + 1) * nbhoriz + i].node;
 313       SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
 314       if (face) {
 315         meshDS->SetMeshElementOnShape(face, geomFaceID);
 316       }
 317     }
 318   }
 319
 320   const vector<UVPtStruct>& uv_e0 = quad->side[0]->GetUVPtStruct(true,0);
 321   const vector<UVPtStruct>& uv_e1 = quad->side[1]->GetUVPtStruct(false,1);
 322   const vector<UVPtStruct>& uv_e2 = quad->side[2]->GetUVPtStruct(true,1);
 323   const vector<UVPtStruct>& uv_e3 = quad->side[3]->GetUVPtStruct(false,0);
 324
 325   if (uv_e0.empty() || uv_e1.empty() || uv_e2.empty() || uv_e3.empty())
 326     return error(COMPERR_BAD_INPUT_MESH);
 327
 328   double eps = Precision::Confusion();
 329
 330   // Boundary quadrangles
 331
 332   if (quad->isEdgeOut[0]) {
 333     // Down edge is out
 334     //
 335     // |___|___|___|___|___|___|
 336     // |   |   |   |   |   |   |
 337     // |___|___|___|___|___|___|
 338     // |   |   |   |   |   |   |
 339     // |___|___|___|___|___|___| __ first row of the regular grid
 340     // .  .  .  .  .  .  .  .  . __ down edge nodes
 341     //
 342     // >->->->->->->->->->->->-> -- direction of processing
 343
 344     int g = 0; // number of last processed node in the regular grid
 345
 346     // number of last node of the down edge to be processed
 347     int stop = nbdown - 1;
 348     // if right edge is out, we will stop at a node, previous to the last one
 349     if (quad->isEdgeOut[1]) stop--;
 350
 351     // for each node of the down edge find nearest node
 352     // in the first row of the regular grid and link them
 353     for (i = 0; i < stop; i++) {
 354       const SMDS_MeshNode *a, *b, *c, *d;
 355       a = uv_e0[i].node;
 356       b = uv_e0[i + 1].node;
 357       gp_Pnt pb (b->X(), b->Y(), b->Z());
 358
 359       // find node c in the regular grid, which will be linked with node b
 360       int near = g;
 361       if (i == stop - 1) {
 362         // right bound reached, link with the rightmost node
 363         near = iup;
 364         c = quad->uv_grid[nbhoriz + iup].node;
 365       }
 366       else {
 367         // find in the grid node c, nearest to the b
 368         double mind = RealLast();
 369         for (int k = g; k <= iup; k++) {
 370
 371           const SMDS_MeshNode *nk;
 372           if (k < ilow) // this can be, if left edge is out
 373             nk = uv_e3[1].node; // get node from the left edge
 374           else
 375             nk = quad->uv_grid[nbhoriz + k].node; // get one of middle nodes
 376
 377           gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 378           double dist = pb.Distance(pnk);
 379           if (dist < mind - eps) {
 380             c = nk;
 381             near = k;
 382             mind = dist;
 383           } else {
 384             break;
 385           }
 386         }
 387       }
 388
 389       if (near == g) { // make triangle
 390         SMDS_MeshFace* face = myHelper->AddFace(a, b, c);
 391         if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 392       }
 393       else { // make quadrangle
 394         if (near - 1 < ilow)
 395           d = uv_e3[1].node;
 396         else
 397           d = quad->uv_grid[nbhoriz + near - 1].node;
 398         //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 399
 400         if (!myTrianglePreference){
 401           SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
 402           if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 403         }
 404         else {
 405           SplitQuad(meshDS, geomFaceID, a, b, c, d);
 406         }
 407
 408         // if node d is not at position g - make additional triangles
 409         if (near - 1 > g) {
 410           for (int k = near - 1; k > g; k--) {
 411             c = quad->uv_grid[nbhoriz + k].node;
 412             if (k - 1 < ilow)
 413               d = uv_e3[1].node;
 414             else
 415               d = quad->uv_grid[nbhoriz + k - 1].node;
 416             SMDS_MeshFace* face = myHelper->AddFace(a, c, d);
 417             if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 418           }
 419         }
 420         g = near;
 421       }
 422     }
 423   } else {
 424     if (quad->isEdgeOut[2]) {
 425       // Up edge is out
 426       //
 427       // <-<-<-<-<-<-<-<-<-<-<-<-< -- direction of processing
 428       //
 429       // .  .  .  .  .  .  .  .  . __ up edge nodes
 430       //  ___ ___ ___ ___ ___ ___  __ first row of the regular grid
 431       // |   |   |   |   |   |   |
 432       // |___|___|___|___|___|___|
 433       // |   |   |   |   |   |   |
 434       // |___|___|___|___|___|___|
 435       // |   |   |   |   |   |   |
 436
 437       int g = nbhoriz - 1; // last processed node in the regular grid
 438
 439       int stop = 0;
 440       // if left edge is out, we will stop at a second node
 441       if (quad->isEdgeOut[3]) stop++;
 442
 443       // for each node of the up edge find nearest node
 444       // in the first row of the regular grid and link them
 445       for (i = nbup - 1; i > stop; i--) {
 446         const SMDS_MeshNode *a, *b, *c, *d;
 447         a = uv_e2[i].node;
 448         b = uv_e2[i - 1].node;
 449         gp_Pnt pb (b->X(), b->Y(), b->Z());
 450
 451         // find node c in the grid, which will be linked with node b
 452         int near = g;
 453         if (i == stop + 1) { // left bound reached, link with the leftmost node
 454           c = quad->uv_grid[nbhoriz*(nbvertic - 2) + ilow].node;
 455           near = ilow;
 456         } else {
 457           // find node c in the grid, nearest to the b
 458           double mind = RealLast();
 459           for (int k = g; k >= ilow; k--) {
 460             const SMDS_MeshNode *nk;
 461             if (k > iup)
 462               nk = uv_e1[nbright - 2].node;
 463             else
 464               nk = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
 465             gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 466             double dist = pb.Distance(pnk);
 467             if (dist < mind - eps) {
 468               c = nk;
 469               near = k;
 470               mind = dist;
 471             } else {
 472               break;
 473             }
 474           }
 475         }
 476
 477         if (near == g) { // make triangle
 478           SMDS_MeshFace* face = myHelper->AddFace(a, b, c);
 479           if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 480         }
 481         else { // make quadrangle
 482           if (near + 1 > iup)
 483             d = uv_e1[nbright - 2].node;
 484           else
 485             d = quad->uv_grid[nbhoriz*(nbvertic - 2) + near + 1].node;
 486           //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 487           if (!myTrianglePreference){
 488             SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
 489             if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 490           }
 491           else {
 492             SplitQuad(meshDS, geomFaceID, a, b, c, d);
 493           }
 494
 495           if (near + 1 < g) { // if d not is at g - make additional triangles
 496             for (int k = near + 1; k < g; k++) {
 497               c = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
 498               if (k + 1 > iup)
 499                 d = uv_e1[nbright - 2].node;
 500               else
 501                 d = quad->uv_grid[nbhoriz*(nbvertic - 2) + k + 1].node;
 502               SMDS_MeshFace* face = myHelper->AddFace(a, c, d);
 503               if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 504             }
 505           }
 506           g = near;
 507         }
 508       }
 509     }
 510   }
 511
 512   // right or left boundary quadrangles
 513   if (quad->isEdgeOut[1]) {
 514 //    MESSAGE("right edge is out");
 515     int g = 0; // last processed node in the grid
 516     int stop = nbright - 1;
 517     if (quad->isEdgeOut[2]) stop--;
 518     for (i = 0; i < stop; i++) {
 519       const SMDS_MeshNode *a, *b, *c, *d;
 520       a = uv_e1[i].node;
 521       b = uv_e1[i + 1].node;
 522       gp_Pnt pb (b->X(), b->Y(), b->Z());
 523
 524       // find node c in the grid, nearest to the b
 525       int near = g;
 526       if (i == stop - 1) { // up bondary reached
 527         c = quad->uv_grid[nbhoriz*(jup + 1) - 2].node;
 528         near = jup;
 529       } else {
 530         double mind = RealLast();
 531         for (int k = g; k <= jup; k++) {
 532           const SMDS_MeshNode *nk;
 533           if (k < jlow)
 534             nk = uv_e0[nbdown - 2].node;
 535           else
 536             nk = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
 537           gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 538           double dist = pb.Distance(pnk);
 539           if (dist < mind - eps) {
 540             c = nk;
 541             near = k;
 542             mind = dist;
 543           } else {
 544             break;
 545           }
 546         }
 547       }
 548
 549       if (near == g) { // make triangle
 550         SMDS_MeshFace* face = myHelper->AddFace(a, b, c);
 551         if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 552       }
 553       else { // make quadrangle
 554         if (near - 1 < jlow)
 555           d = uv_e0[nbdown - 2].node;
 556         else
 557           d = quad->uv_grid[nbhoriz*near - 2].node;
 558         //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 559
 560         if (!myTrianglePreference){
 561           SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
 562           if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 563         }
 564         else {
 565           SplitQuad(meshDS, geomFaceID, a, b, c, d);
 566         }
 567
 568         if (near - 1 > g) { // if d not is at g - make additional triangles
 569           for (int k = near - 1; k > g; k--) {
 570             c = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
 571             if (k - 1 < jlow)
 572               d = uv_e0[nbdown - 2].node;
 573             else
 574               d = quad->uv_grid[nbhoriz*k - 2].node;
 575             SMDS_MeshFace* face = myHelper->AddFace(a, c, d);
 576             if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 577           }
 578         }
 579         g = near;
 580       }
 581     }
 582   } else {
 583     if (quad->isEdgeOut[3]) {
 584 //      MESSAGE("left edge is out");
 585       int g = nbvertic - 1; // last processed node in the grid
 586       int stop = 0;
 587       if (quad->isEdgeOut[0]) stop++;
 588       for (i = nbleft - 1; i > stop; i--) {
 589         const SMDS_MeshNode *a, *b, *c, *d;
 590         a = uv_e3[i].node;
 591         b = uv_e3[i - 1].node;
 592         gp_Pnt pb (b->X(), b->Y(), b->Z());
 593
 594         // find node c in the grid, nearest to the b
 595         int near = g;
 596         if (i == stop + 1) { // down bondary reached
 597           c = quad->uv_grid[nbhoriz*jlow + 1].node;
 598           near = jlow;
 599         } else {
 600           double mind = RealLast();
 601           for (int k = g; k >= jlow; k--) {
 602             const SMDS_MeshNode *nk;
 603             if (k > jup)
 604               nk = uv_e2[1].node;
 605             else
 606               nk = quad->uv_grid[nbhoriz*k + 1].node;
 607             gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 608             double dist = pb.Distance(pnk);
 609             if (dist < mind - eps) {
 610               c = nk;
 611               near = k;
 612               mind = dist;
 613             } else {
 614               break;
 615             }
 616           }
 617         }
 618
 619         if (near == g) { // make triangle
 620           SMDS_MeshFace* face = myHelper->AddFace(a, b, c);
 621           if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 622         }
 623         else { // make quadrangle
 624           if (near + 1 > jup)
 625             d = uv_e2[1].node;
 626           else
 627             d = quad->uv_grid[nbhoriz*(near + 1) + 1].node;
 628           //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 629           if (!myTrianglePreference){
 630             SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
 631             if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 632           }
 633           else {
 634             SplitQuad(meshDS, geomFaceID, a, b, c, d);
 635           }
 636
 637           if (near + 1 < g) { // if d not is at g - make additional triangles
 638             for (int k = near + 1; k < g; k++) {
 639               c = quad->uv_grid[nbhoriz*k + 1].node;
 640               if (k + 1 > jup)
 641                 d = uv_e2[1].node;
 642               else
 643                 d = quad->uv_grid[nbhoriz*(k + 1) + 1].node;
 644               SMDS_MeshFace* face = myHelper->AddFace(a, c, d);
 645               if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 646             }
 647           }
 648           g = near;
 649         }
 650       }
 651     }
 652   }
 653
 654   if ( myNeedSmooth )
 655     Smooth( quad );
 656
 657   bool isOk = true;
 658   return isOk;
 659 }
 660
 661
 662 //=============================================================================
 663 /*!
 664  *  Evaluate
 665  */
 666 //=============================================================================
 667
 668 bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh& aMesh,
 669                                         const TopoDS_Shape& aShape,
 670                                         MapShapeNbElems& aResMap)
 671
 672 {
 673   aMesh.GetSubMesh(aShape);
 674
 675   std::vector<int> aNbNodes(4);
 676   bool IsQuadratic = false;
 677   if (!CheckNbEdgesForEvaluate(aMesh, aShape, aResMap, aNbNodes, IsQuadratic)) {
 678     std::vector<int> aResVec(SMDSEntity_Last);
 679     for (int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
 680     SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 681     aResMap.insert(std::make_pair(sm,aResVec));
 682     SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 683     smError.reset(new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
 684     return false;
 685   }
 686
 687   if (myQuadranglePreference) {
 688     int n1 = aNbNodes[0];
 689     int n2 = aNbNodes[1];
 690     int n3 = aNbNodes[2];
 691     int n4 = aNbNodes[3];
 692     int nfull = n1+n2+n3+n4;
 693     int ntmp = nfull/2;
 694     ntmp = ntmp*2;
 695     if (nfull==ntmp && ((n1!=n3) || (n2!=n4))) {
 696       // special path for using only quandrangle faces
 697       return EvaluateQuadPref(aMesh, aShape, aNbNodes, aResMap, IsQuadratic);
 698       //return true;
 699     }
 700   }
 701
 702   int nbdown  = aNbNodes[0];
 703   int nbup    = aNbNodes[2];
 704
 705   int nbright = aNbNodes[1];
 706   int nbleft  = aNbNodes[3];
 707
 708   int nbhoriz  = Min(nbdown, nbup);
 709   int nbvertic = Min(nbright, nbleft);
 710
 711   int dh = Max(nbdown, nbup) - nbhoriz;
 712   int dv = Max(nbright, nbleft) - nbvertic;
 713
 714   //int kdh = 0;
 715   //if (dh>0) kdh = 1;
 716   //int kdv = 0;
 717   //if (dv>0) kdv = 1;
 718
 719   int nbNodes = (nbhoriz-2)*(nbvertic-2);
 720   //int nbFaces3 = dh + dv + kdh*(nbvertic-1)*2 + kdv*(nbhoriz-1)*2;
 721   int nbFaces3 = dh + dv;
 722   //if (kdh==1 && kdv==1) nbFaces3 -= 2;
 723   //if (dh>0 && dv>0) nbFaces3 -= 2;
 724   //int nbFaces4 = (nbhoriz-1-kdh)*(nbvertic-1-kdv);
 725   int nbFaces4 = (nbhoriz-1)*(nbvertic-1);
 726
 727   std::vector<int> aVec(SMDSEntity_Last);
 728   for (int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i] = 0;
 729   if (IsQuadratic) {
 730     aVec[SMDSEntity_Quad_Triangle] = nbFaces3;
 731     aVec[SMDSEntity_Quad_Quadrangle] = nbFaces4;
 732     int nbbndedges = nbdown + nbup + nbright + nbleft -4;
 733     int nbintedges = (nbFaces4*4 + nbFaces3*3 - nbbndedges) / 2;
 734     aVec[SMDSEntity_Node] = nbNodes + nbintedges;
 735     if (aNbNodes.size()==5) {
 736       aVec[SMDSEntity_Quad_Triangle] = nbFaces3 + aNbNodes[3] -1;
 737       aVec[SMDSEntity_Quad_Quadrangle] = nbFaces4 - aNbNodes[3] +1;
 738     }
 739   }
 740   else {
 741     aVec[SMDSEntity_Node] = nbNodes;
 742     aVec[SMDSEntity_Triangle] = nbFaces3;
 743     aVec[SMDSEntity_Quadrangle] = nbFaces4;
 744     if (aNbNodes.size()==5) {
 745       aVec[SMDSEntity_Triangle] = nbFaces3 + aNbNodes[3] - 1;
 746       aVec[SMDSEntity_Quadrangle] = nbFaces4 - aNbNodes[3] + 1;
 747     }
 748   }
 749   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 750   aResMap.insert(std::make_pair(sm,aVec));
 751
 752   return true;
 753 }
 754
 755
 756 //================================================================================
 757 /*!
 758  * \brief Return true if only two given edges meat at their common vertex
 759  */
 760 //================================================================================
 761
 762 static bool twoEdgesMeatAtVertex(const TopoDS_Edge& e1,
 763                                  const TopoDS_Edge& e2,
 764                                  SMESH_Mesh &       mesh)
 765 {
 766   TopoDS_Vertex v;
 767   if (!TopExp::CommonVertex(e1, e2, v))
 768     return false;
 769   TopTools_ListIteratorOfListOfShape ancestIt(mesh.GetAncestors(v));
 770   for (; ancestIt.More() ; ancestIt.Next())
 771     if (ancestIt.Value().ShapeType() == TopAbs_EDGE)
 772       if (!e1.IsSame(ancestIt.Value()) && !e2.IsSame(ancestIt.Value()))
 773         return false;
 774   return true;
 775 }
 776
 777 //=============================================================================
 778 /*!
 779  *
 780  */
 781 //=============================================================================
 782
 783 FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh &         aMesh,
 784                                                        const TopoDS_Shape & aShape)
 785   //throw(SALOME_Exception)
 786 {
 787   TopoDS_Face F = TopoDS::Face(aShape);
 788   if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD );
 789   const bool ignoreMediumNodes = _quadraticMesh;
 790
 791   // verify 1 wire only, with 4 edges
 792   TopoDS_Vertex V;
 793   list< TopoDS_Edge > edges;
 794   list< int > nbEdgesInWire;
 795   int nbWire = SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
 796   if (nbWire != 1) {
 797     error(COMPERR_BAD_SHAPE, TComm("Wrong number of wires: ") << nbWire);
 798     return 0;
 799   }
 800   FaceQuadStruct* quad = new FaceQuadStruct;
 801   quad->uv_grid = 0;
 802   quad->side.reserve(nbEdgesInWire.front());
 803   quad->face = F;
 804
 805   int nbSides = 0;
 806   list< TopoDS_Edge >::iterator edgeIt = edges.begin();
 807   if (nbEdgesInWire.front() == 3) // exactly 3 edges
 808   {
 809     SMESH_Comment comment;
 810     SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 811     if (myTriaVertexID == -1)
 812     {
 813       comment << "No Base vertex parameter provided for a trilateral geometrical face";
 814     }
 815     else
 816     {
 817       TopoDS_Vertex V = TopoDS::Vertex(meshDS->IndexToShape(myTriaVertexID));
 818       if (!V.IsNull()) {
 819         TopoDS_Edge E1,E2,E3;
 820         for (; edgeIt != edges.end(); ++edgeIt) {
 821           TopoDS_Edge E =  *edgeIt;
 822           TopoDS_Vertex VF, VL;
 823           TopExp::Vertices(E, VF, VL, true);
 824           if (VF.IsSame(V))
 825             E1 = E;
 826           else if (VL.IsSame(V))
 827             E3 = E;
 828           else
 829             E2 = E;
 830         }
 831         if (!E1.IsNull() && !E2.IsNull() && !E3.IsNull())
 832         {
 833           quad->side.push_back(new StdMeshers_FaceSide(F, E1, &aMesh, true, ignoreMediumNodes));
 834           quad->side.push_back(new StdMeshers_FaceSide(F, E2, &aMesh, true, ignoreMediumNodes));
 835           quad->side.push_back(new StdMeshers_FaceSide(F, E3, &aMesh, false,ignoreMediumNodes));
 836           const vector<UVPtStruct>& UVPSleft  = quad->side[0]->GetUVPtStruct(true,0);
 837           /*  vector<UVPtStruct>& UVPStop   = */quad->side[1]->GetUVPtStruct(false,1);
 838           /*  vector<UVPtStruct>& UVPSright = */quad->side[2]->GetUVPtStruct(true,1);
 839           const SMDS_MeshNode* aNode = UVPSleft[0].node;
 840           gp_Pnt2d aPnt2d(UVPSleft[0].u, UVPSleft[0].v);
 841           quad->side.push_back(new StdMeshers_FaceSide(aNode, aPnt2d, quad->side[1]));
 842           return quad;
 843         }
 844       }
 845       comment << "Invalid Base vertex parameter: " << myTriaVertexID << " is not among [";
 846       TopTools_MapOfShape vMap;
 847       for (TopExp_Explorer v(aShape, TopAbs_VERTEX); v.More(); v.Next())
 848         if (vMap.Add(v.Current()))
 849           comment << meshDS->ShapeToIndex(v.Current()) << (vMap.Extent()==3 ? "]" : ", ");
 850     }
 851     error(comment);
 852     delete quad;
 853     return quad = 0;
 854   }
 855   else if (nbEdgesInWire.front() == 4) // exactly 4 edges
 856   {
 857     for (; edgeIt != edges.end(); ++edgeIt, nbSides++)
 858       quad->side.push_back(new StdMeshers_FaceSide(F, *edgeIt, &aMesh,
 859                                                     nbSides<TOP_SIDE, ignoreMediumNodes));
 860   }
 861   else if (nbEdgesInWire.front() > 4) // more than 4 edges - try to unite some
 862   {
 863     list< TopoDS_Edge > sideEdges;
 864     vector< int > degenSides;
 865     while (!edges.empty()) {
 866       sideEdges.clear();
 867       sideEdges.splice(sideEdges.end(), edges, edges.begin()); // edges.front() -> sideEdges.end()
 868       bool sameSide = true;
 869       while (!edges.empty() && sameSide) {
 870         sameSide = SMESH_Algo::IsContinuous(sideEdges.back(), edges.front());
 871         if (sameSide)
 872           sideEdges.splice(sideEdges.end(), edges, edges.begin());
 873       }
 874       if (nbSides == 0) { // go backward from the first edge
 875         sameSide = true;
 876         while (!edges.empty() && sameSide) {
 877           sameSide = SMESH_Algo::IsContinuous(sideEdges.front(), edges.back());
 878           if (sameSide)
 879             sideEdges.splice(sideEdges.begin(), edges, --edges.end());
 880         }
 881       }
 882       if ( sideEdges.size() == 1 && BRep_Tool::Degenerated( sideEdges.front() ))
 883         degenSides.push_back( nbSides );
 884
 885       quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh,
 886                                                     nbSides<TOP_SIDE, ignoreMediumNodes));
 887       ++nbSides;
 888     }
 889     if ( !degenSides.empty() && nbSides - degenSides.size() == 4 )
 890     {
 891       myNeedSmooth = true;
 892       for ( unsigned i = TOP_SIDE; i < quad->side.size(); ++i )
 893         quad->side[i]->Reverse();
 894
 895       for ( int i = degenSides.size()-1; i > -1; --i )
 896       {
 897         StdMeshers_FaceSide * & degenSide = quad->side[ degenSides[ i ]];
 898         delete degenSide;
 899         quad->side.erase( vector<StdMeshers_FaceSide*>::iterator( & degenSide ));
 900       }
 901       for ( unsigned i = TOP_SIDE; i < quad->side.size(); ++i )
 902         quad->side[i]->Reverse();
 903
 904       nbSides -= degenSides.size();
 905     }
 906     // issue 20222. Try to unite only edges shared by two same faces
 907     if (nbSides < 4) {
 908       // delete found sides
 909       { FaceQuadStruct cleaner(*quad); }
 910       quad->side.clear();
 911       quad->side.reserve(nbEdgesInWire.front());
 912       nbSides = 0;
 913
 914       SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
 915       while (!edges.empty()) {
 916         sideEdges.clear();
 917         sideEdges.splice(sideEdges.end(), edges, edges.begin());
 918         bool sameSide = true;
 919         while (!edges.empty() && sameSide) {
 920           sameSide =
 921             SMESH_Algo::IsContinuous(sideEdges.back(), edges.front()) &&
 922             twoEdgesMeatAtVertex(sideEdges.back(), edges.front(), aMesh);
 923           if (sameSide)
 924             sideEdges.splice(sideEdges.end(), edges, edges.begin());
 925         }
 926         if (nbSides == 0) { // go backward from the first edge
 927           sameSide = true;
 928           while (!edges.empty() && sameSide) {
 929             sameSide =
 930               SMESH_Algo::IsContinuous(sideEdges.front(), edges.back()) &&
 931               twoEdgesMeatAtVertex(sideEdges.front(), edges.back(), aMesh);
 932             if (sameSide)
 933               sideEdges.splice(sideEdges.begin(), edges, --edges.end());
 934           }
 935         }
 936         quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh,
 937                                                       nbSides<TOP_SIDE, ignoreMediumNodes));
 938         ++nbSides;
 939       }
 940     }
 941   }
 942   if (nbSides != 4) {
 943 #ifdef _DEBUG_
 944     MESSAGE ("StdMeshers_Quadrangle_2D. Edge IDs of " << nbSides << " sides:\n");
 945     for (int i = 0; i < nbSides; ++i) {
 946       MESSAGE (" (");
 947       for (int e = 0; e < quad->side[i]->NbEdges(); ++e)
 948         MESSAGE (myHelper->GetMeshDS()->ShapeToIndex(quad->side[i]->Edge(e)) << " ");
 949       MESSAGE (")\n");
 950     }
 951     //cout << endl;
 952 #endif
 953     if (!nbSides)
 954       nbSides = nbEdgesInWire.front();
 955     error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides);
 956     delete quad;
 957     quad = 0;
 958   }
 959
 960   return quad;
 961 }
 962
 963
 964 //=============================================================================
 965 /*!
 966  *
 967  */
 968 //=============================================================================
 969
 970 bool StdMeshers_Quadrangle_2D::CheckNbEdgesForEvaluate(SMESH_Mesh& aMesh,
 971                                                        const TopoDS_Shape & aShape,
 972                                                        MapShapeNbElems& aResMap,
 973                                                        std::vector<int>& aNbNodes,
 974                                                        bool& IsQuadratic)
 975
 976 {
 977   const TopoDS_Face & F = TopoDS::Face(aShape);
 978
 979   // verify 1 wire only, with 4 edges
 980   TopoDS_Vertex V;
 981   list< TopoDS_Edge > edges;
 982   list< int > nbEdgesInWire;
 983   int nbWire = SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
 984   if (nbWire != 1) {
 985     return false;
 986   }
 987
 988   aNbNodes.resize(4);
 989
 990   int nbSides = 0;
 991   list< TopoDS_Edge >::iterator edgeIt = edges.begin();
 992   SMESH_subMesh * sm = aMesh.GetSubMesh(*edgeIt);
 993   MapShapeNbElemsItr anIt = aResMap.find(sm);
 994   if (anIt==aResMap.end()) {
 995     return false;
 996   }
 997   std::vector<int> aVec = (*anIt).second;
 998   IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
 999   if (nbEdgesInWire.front() == 3) { // exactly 3 edges
1000     if (myTriaVertexID>0) {
1001       SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
1002       TopoDS_Vertex V = TopoDS::Vertex(meshDS->IndexToShape(myTriaVertexID));
1003       if (!V.IsNull()) {
1004         TopoDS_Edge E1,E2,E3;
1005         for (; edgeIt != edges.end(); ++edgeIt) {
1006           TopoDS_Edge E =  TopoDS::Edge(*edgeIt);
1007           TopoDS_Vertex VF, VL;
1008           TopExp::Vertices(E, VF, VL, true);
1009           if (VF.IsSame(V))
1010             E1 = E;
1011           else if (VL.IsSame(V))
1012             E3 = E;
1013           else
1014             E2 = E;
1015         }
1016         SMESH_subMesh * sm = aMesh.GetSubMesh(E1);
1017         MapShapeNbElemsItr anIt = aResMap.find(sm);
1018         if (anIt==aResMap.end()) return false;
1019         std::vector<int> aVec = (*anIt).second;
1020         if (IsQuadratic)
1021           aNbNodes[0] = (aVec[SMDSEntity_Node]-1)/2 + 2;
1022         else
1023           aNbNodes[0] = aVec[SMDSEntity_Node] + 2;
1024         sm = aMesh.GetSubMesh(E2);
1025         anIt = aResMap.find(sm);
1026         if (anIt==aResMap.end()) return false;
1027         aVec = (*anIt).second;
1028         if (IsQuadratic)
1029           aNbNodes[1] = (aVec[SMDSEntity_Node]-1)/2 + 2;
1030         else
1031           aNbNodes[1] = aVec[SMDSEntity_Node] + 2;
1032         sm = aMesh.GetSubMesh(E3);
1033         anIt = aResMap.find(sm);
1034         if (anIt==aResMap.end()) return false;
1035         aVec = (*anIt).second;
1036         if (IsQuadratic)
1037           aNbNodes[2] = (aVec[SMDSEntity_Node]-1)/2 + 2;
1038         else
1039           aNbNodes[2] = aVec[SMDSEntity_Node] + 2;
1040         aNbNodes[3] = aNbNodes[1];
1041         aNbNodes.resize(5);
1042         nbSides = 4;
1043       }
1044     }
1045   }
1046   if (nbEdgesInWire.front() == 4) { // exactly 4 edges
1047     for (; edgeIt != edges.end(); edgeIt++) {
1048       SMESH_subMesh * sm = aMesh.GetSubMesh(*edgeIt);
1049       MapShapeNbElemsItr anIt = aResMap.find(sm);
1050       if (anIt==aResMap.end()) {
1051         return false;
1052       }
1053       std::vector<int> aVec = (*anIt).second;
1054       if (IsQuadratic)
1055         aNbNodes[nbSides] = (aVec[SMDSEntity_Node]-1)/2 + 2;
1056       else
1057         aNbNodes[nbSides] = aVec[SMDSEntity_Node] + 2;
1058       nbSides++;
1059     }
1060   }
1061   else if (nbEdgesInWire.front() > 4) { // more than 4 edges - try to unite some
1062     list< TopoDS_Edge > sideEdges;
1063     while (!edges.empty()) {
1064       sideEdges.clear();
1065       sideEdges.splice(sideEdges.end(), edges, edges.begin()); // edges.front() -> sideEdges.end()
1066       bool sameSide = true;
1067       while (!edges.empty() && sameSide) {
1068         sameSide = SMESH_Algo::IsContinuous(sideEdges.back(), edges.front());
1069         if (sameSide)
1070           sideEdges.splice(sideEdges.end(), edges, edges.begin());
1071       }
1072       if (nbSides == 0) { // go backward from the first edge
1073         sameSide = true;
1074         while (!edges.empty() && sameSide) {
1075           sameSide = SMESH_Algo::IsContinuous(sideEdges.front(), edges.back());
1076           if (sameSide)
1077             sideEdges.splice(sideEdges.begin(), edges, --edges.end());
1078         }
1079       }
1080       list<TopoDS_Edge>::iterator ite = sideEdges.begin();
1081       aNbNodes[nbSides] = 1;
1082       for (; ite!=sideEdges.end(); ite++) {
1083         SMESH_subMesh * sm = aMesh.GetSubMesh(*ite);
1084         MapShapeNbElemsItr anIt = aResMap.find(sm);
1085         if (anIt==aResMap.end()) {
1086           return false;
1087         }
1088         std::vector<int> aVec = (*anIt).second;
1089         if (IsQuadratic)
1090           aNbNodes[nbSides] += (aVec[SMDSEntity_Node]-1)/2 + 1;
1091         else
1092           aNbNodes[nbSides] += aVec[SMDSEntity_Node] + 1;
1093       }
1094       ++nbSides;
1095     }
1096     // issue 20222. Try to unite only edges shared by two same faces
1097     if (nbSides < 4) {
1098       nbSides = 0;
1099       SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
1100       while (!edges.empty()) {
1101         sideEdges.clear();
1102         sideEdges.splice(sideEdges.end(), edges, edges.begin());
1103         bool sameSide = true;
1104         while (!edges.empty() && sameSide) {
1105           sameSide =
1106             SMESH_Algo::IsContinuous(sideEdges.back(), edges.front()) &&
1107             twoEdgesMeatAtVertex(sideEdges.back(), edges.front(), aMesh);
1108           if (sameSide)
1109             sideEdges.splice(sideEdges.end(), edges, edges.begin());
1110         }
1111         if (nbSides == 0) { // go backward from the first edge
1112           sameSide = true;
1113           while (!edges.empty() && sameSide) {
1114             sameSide =
1115               SMESH_Algo::IsContinuous(sideEdges.front(), edges.back()) &&
1116               twoEdgesMeatAtVertex(sideEdges.front(), edges.back(), aMesh);
1117             if (sameSide)
1118               sideEdges.splice(sideEdges.begin(), edges, --edges.end());
1119           }
1120         }
1121         list<TopoDS_Edge>::iterator ite = sideEdges.begin();
1122         aNbNodes[nbSides] = 1;
1123         for (; ite!=sideEdges.end(); ite++) {
1124           SMESH_subMesh * sm = aMesh.GetSubMesh(*ite);
1125           MapShapeNbElemsItr anIt = aResMap.find(sm);
1126           if (anIt==aResMap.end()) {
1127             return false;
1128           }
1129           std::vector<int> aVec = (*anIt).second;
1130           if (IsQuadratic)
1131             aNbNodes[nbSides] += (aVec[SMDSEntity_Node]-1)/2 + 1;
1132           else
1133             aNbNodes[nbSides] += aVec[SMDSEntity_Node] + 1;
1134         }
1135         ++nbSides;
1136       }
1137     }
1138   }
1139   if (nbSides != 4) {
1140     if (!nbSides)
1141       nbSides = nbEdgesInWire.front();
1142     error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides);
1143     return false;
1144   }
1145
1146   return true;
1147 }
1148
1149
1150 //=============================================================================
1151 /*!
1152  *  CheckAnd2Dcompute
1153  */
1154 //=============================================================================
1155
1156 FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute
1157                            (SMESH_Mesh &         aMesh,
1158                             const TopoDS_Shape & aShape,
1159                             const bool           CreateQuadratic) //throw(SALOME_Exception)
1160 {
1161   _quadraticMesh = CreateQuadratic;
1162
1163   FaceQuadStruct *quad = CheckNbEdges(aMesh, aShape);
1164
1165   if (!quad) return 0;
1166
1167   // set normalized grid on unit square in parametric domain
1168   bool stat = SetNormalizedGrid(aMesh, aShape, quad);
1169   if (!stat) {
1170     if (quad) delete quad;
1171     quad = 0;
1172   }
1173
1174   return quad;
1175 }
1176
1177 //=============================================================================
1178 /*!
1179  *
1180  */
1181 //=============================================================================
1182
1183 faceQuadStruct::~faceQuadStruct()
1184 {
1185   for (int i = 0; i < side.size(); i++) {
1186     if (side[i])     delete side[i];
1187   }
1188   if (uv_grid)       delete [] uv_grid;
1189 }
1190
1191 namespace {
1192   inline const vector<UVPtStruct>& GetUVPtStructIn(FaceQuadStruct* quad, int i, int nbSeg)
1193   {
1194     bool   isXConst   = (i == BOTTOM_SIDE || i == TOP_SIDE);
1195     double constValue = (i == BOTTOM_SIDE || i == LEFT_SIDE) ? 0 : 1;
1196     return
1197       quad->isEdgeOut[i] ?
1198       quad->side[i]->SimulateUVPtStruct(nbSeg,isXConst,constValue) :
1199       quad->side[i]->GetUVPtStruct(isXConst,constValue);
1200   }
1201   inline gp_UV CalcUV(double x, double y,
1202                       const gp_UV& a0,const gp_UV& a1,const gp_UV& a2,const gp_UV& a3,
1203                       const gp_UV& p0,const gp_UV& p1,const gp_UV& p2,const gp_UV& p3)
1204   {
1205     return
1206       ((1 - y) * p0 + x * p1 + y * p2 + (1 - x) * p3 ) -
1207       ((1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3);
1208   }
1209 }
1210
1211 //=============================================================================
1212 /*!
1213  *
1214  */
1215 //=============================================================================
1216
1217 bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh,
1218                                                   const TopoDS_Shape& aShape,
1219                                                   FaceQuadStruct* & quad) //throw (SALOME_Exception)
1220 {
1221   // Algorithme décrit dans "Génération automatique de maillages"
1222   // P.L. GEORGE, MASSON, § 6.4.1 p. 84-85
1223   // traitement dans le domaine paramétrique 2d u,v
1224   // transport - projection sur le carré unité
1225
1226 //  MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid");
1227 //  const TopoDS_Face& F = TopoDS::Face(aShape);
1228
1229   // 1 --- find orientation of the 4 edges, by test on extrema
1230
1231   //      max             min                    0     x1     1
1232   //     |<----north-2-------^                a3 -------------> a2
1233   //     |                   |                   ^1          1^
1234   //    west-3            east-1 =right          |            |
1235   //     |                   |         ==>       |            |
1236   //  y0 |                   | y1                |            |
1237   //     |                   |                   |0          0|
1238   //     v----south-0-------->                a0 -------------> a1
1239   //      min             max                    0     x0     1
1240   //             =down
1241   //
1242
1243   // 3 --- 2D normalized values on unit square [0..1][0..1]
1244
1245   int nbhoriz  = Min(quad->side[0]->NbPoints(), quad->side[2]->NbPoints());
1246   int nbvertic = Min(quad->side[1]->NbPoints(), quad->side[3]->NbPoints());
1247
1248   quad->isEdgeOut[0] = (quad->side[0]->NbPoints() > quad->side[2]->NbPoints());
1249   quad->isEdgeOut[1] = (quad->side[1]->NbPoints() > quad->side[3]->NbPoints());
1250   quad->isEdgeOut[2] = (quad->side[2]->NbPoints() > quad->side[0]->NbPoints());
1251   quad->isEdgeOut[3] = (quad->side[3]->NbPoints() > quad->side[1]->NbPoints());
1252
1253   UVPtStruct *uv_grid = quad->uv_grid = new UVPtStruct[nbvertic * nbhoriz];
1254
1255   const vector<UVPtStruct>& uv_e0 = GetUVPtStructIn(quad, 0, nbhoriz - 1);
1256   const vector<UVPtStruct>& uv_e1 = GetUVPtStructIn(quad, 1, nbvertic - 1);
1257   const vector<UVPtStruct>& uv_e2 = GetUVPtStructIn(quad, 2, nbhoriz - 1);
1258   const vector<UVPtStruct>& uv_e3 = GetUVPtStructIn(quad, 3, nbvertic - 1);
1259
1260   if (uv_e0.empty() || uv_e1.empty() || uv_e2.empty() || uv_e3.empty())
1261     //return error("Can't find nodes on sides");
1262     return error(COMPERR_BAD_INPUT_MESH);
1263
1264   if ( myNeedSmooth )
1265     UpdateDegenUV( quad );
1266
1267   // nodes Id on "in" edges
1268   if (! quad->isEdgeOut[0]) {
1269     int j = 0;
1270     for (int i = 0; i < nbhoriz; i++) { // down
1271       int ij = j * nbhoriz + i;
1272       uv_grid[ij].node = uv_e0[i].node;
1273     }
1274   }
1275   if (! quad->isEdgeOut[1]) {
1276     int i = nbhoriz - 1;
1277     for (int j = 0; j < nbvertic; j++) { // right
1278       int ij = j * nbhoriz + i;
1279       uv_grid[ij].node = uv_e1[j].node;
1280     }
1281   }
1282   if (! quad->isEdgeOut[2]) {
1283     int j = nbvertic - 1;
1284     for (int i = 0; i < nbhoriz; i++) { // up
1285       int ij = j * nbhoriz + i;
1286       uv_grid[ij].node = uv_e2[i].node;
1287     }
1288   }
1289   if (! quad->isEdgeOut[3]) {
1290     int i = 0;
1291     for (int j = 0; j < nbvertic; j++) { // left
1292       int ij = j * nbhoriz + i;
1293       uv_grid[ij].node = uv_e3[j].node;
1294     }
1295   }
1296
1297   // normalized 2d values on grid
1298   for (int i = 0; i < nbhoriz; i++) {
1299     for (int j = 0; j < nbvertic; j++) {
1300       int ij = j * nbhoriz + i;
1301       // --- droite i cste : x = x0 + y(x1-x0)
1302       double x0 = uv_e0[i].normParam;   // bas - sud
1303       double x1 = uv_e2[i].normParam;   // haut - nord
1304       // --- droite j cste : y = y0 + x(y1-y0)
1305       double y0 = uv_e3[j].normParam;   // gauche-ouest
1306       double y1 = uv_e1[j].normParam;   // droite - est
1307       // --- intersection : x=x0+(y0+x(y1-y0))(x1-x0)
1308       double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
1309       double y = y0 + x * (y1 - y0);
1310       uv_grid[ij].x = x;
1311       uv_grid[ij].y = y;
1312       //MESSAGE("-xy-01 "<<x0<<" "<<x1<<" "<<y0<<" "<<y1);
1313       //MESSAGE("-xy-norm "<<i<<" "<<j<<" "<<x<<" "<<y);
1314     }
1315   }
1316
1317   // 4 --- projection on 2d domain (u,v)
1318   gp_UV a0(uv_e0.front().u, uv_e0.front().v);
1319   gp_UV a1(uv_e0.back().u,  uv_e0.back().v);
1320   gp_UV a2(uv_e2.back().u,  uv_e2.back().v);
1321   gp_UV a3(uv_e2.front().u, uv_e2.front().v);
1322
1323   for (int i = 0; i < nbhoriz; i++) {
1324     for (int j = 0; j < nbvertic; j++) {
1325       int ij = j * nbhoriz + i;
1326       double x = uv_grid[ij].x;
1327       double y = uv_grid[ij].y;
1328       double param_0 = uv_e0[0].normParam + x * (uv_e0.back().normParam - uv_e0[0].normParam); // sud
1329       double param_2 = uv_e2[0].normParam + x * (uv_e2.back().normParam - uv_e2[0].normParam); // nord
1330       double param_1 = uv_e1[0].normParam + y * (uv_e1.back().normParam - uv_e1[0].normParam); // est
1331       double param_3 = uv_e3[0].normParam + y * (uv_e3.back().normParam - uv_e3[0].normParam); // ouest
1332
1333       //MESSAGE("params "<<param_0<<" "<<param_1<<" "<<param_2<<" "<<param_3);
1334       gp_UV p0 = quad->side[0]->Value2d(param_0).XY();
1335       gp_UV p1 = quad->side[1]->Value2d(param_1).XY();
1336       gp_UV p2 = quad->side[2]->Value2d(param_2).XY();
1337       gp_UV p3 = quad->side[3]->Value2d(param_3).XY();
1338
1339       gp_UV uv = CalcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3);
1340
1341       uv_grid[ij].u = uv.X();
1342       uv_grid[ij].v = uv.Y();
1343     }
1344   }
1345   return true;
1346 }
1347
1348 //=======================================================================
1349 //function : ShiftQuad
1350 //purpose  : auxilary function for ComputeQuadPref
1351 //=======================================================================
1352
1353 static void ShiftQuad(FaceQuadStruct* quad, const int num, bool)
1354 {
1355   StdMeshers_FaceSide* side[4] = { quad->side[0], quad->side[1], quad->side[2], quad->side[3] };
1356   for (int i = BOTTOM_SIDE; i < NB_SIDES; ++i) {
1357     int id = (i + num) % NB_SIDES;
1358     bool wasForward = (i < TOP_SIDE);
1359     bool newForward = (id < TOP_SIDE);
1360     if (wasForward != newForward)
1361       side[ i ]->Reverse();
1362     quad->side[ id ] = side[ i ];
1363   }
1364 }
1365
1366 //=======================================================================
1367 //function : CalcUV
1368 //purpose  : auxilary function for ComputeQuadPref
1369 //=======================================================================
1370
1371 static gp_UV CalcUV(double x0, double x1, double y0, double y1,
1372                     FaceQuadStruct* quad,
1373                     const gp_UV& a0, const gp_UV& a1,
1374                     const gp_UV& a2, const gp_UV& a3)
1375 {
1376   const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
1377   const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
1378   const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1);
1379   const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
1380
1381   double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
1382   double y = y0 + x * (y1 - y0);
1383
1384   double param_b = uv_eb[0].normParam + x * (uv_eb.back().normParam - uv_eb[0].normParam);
1385   double param_t = uv_et[0].normParam + x * (uv_et.back().normParam - uv_et[0].normParam);
1386   double param_r = uv_er[0].normParam + y * (uv_er.back().normParam - uv_er[0].normParam);
1387   double param_l = uv_el[0].normParam + y * (uv_el.back().normParam - uv_el[0].normParam);
1388
1389   gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(param_b).XY();
1390   gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(param_r).XY();
1391   gp_UV p2 = quad->side[TOP_SIDE   ]->Value2d(param_t).XY();
1392   gp_UV p3 = quad->side[LEFT_SIDE  ]->Value2d(param_l).XY();
1393
1394   gp_UV uv = CalcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3);
1395
1396   return uv;
1397 }
1398
1399 //=======================================================================
1400 //function : CalcUV2
1401 //purpose  : auxilary function for ComputeQuadPref
1402 //=======================================================================
1403
1404 static gp_UV CalcUV2(double x, double y,
1405                      FaceQuadStruct* quad,
1406                      const gp_UV& a0, const gp_UV& a1,
1407                      const gp_UV& a2, const gp_UV& a3)
1408 {
1409   gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(x).XY();
1410   gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(y).XY();
1411   gp_UV p2 = quad->side[TOP_SIDE   ]->Value2d(x).XY();
1412   gp_UV p3 = quad->side[LEFT_SIDE  ]->Value2d(y).XY();
1413
1414   gp_UV uv = CalcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3);
1415
1416   return uv;
1417 }
1418
1419
1420 //=======================================================================
1421 /*!
1422  * Create only quandrangle faces
1423  */
1424 //=======================================================================
1425
1426 bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh &        aMesh,
1427                                                 const TopoDS_Shape& aShape,
1428                                                 FaceQuadStruct*     quad)
1429 {
1430   // Auxilary key in order to keep old variant
1431   // of meshing after implementation new variant
1432   // for bug 0016220 from Mantis.
1433   bool OldVersion = false;
1434   if (myQuadType == QUAD_QUADRANGLE_PREF_REVERSED)
1435     OldVersion = true;
1436
1437   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
1438   const TopoDS_Face& F = TopoDS::Face(aShape);
1439   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
1440   bool WisF = true;
1441   int i,j,geomFaceID = meshDS->ShapeToIndex(F);
1442
1443   int nb = quad->side[0]->NbPoints();
1444   int nr = quad->side[1]->NbPoints();
1445   int nt = quad->side[2]->NbPoints();
1446   int nl = quad->side[3]->NbPoints();
1447   int dh = abs(nb-nt);
1448   int dv = abs(nr-nl);
1449
1450   if (dh>=dv) {
1451     if (nt>nb) {
1452       // it is a base case => not shift quad but me be replacement is need
1453       ShiftQuad(quad,0,WisF);
1454     }
1455     else {
1456       // we have to shift quad on 2
1457       ShiftQuad(quad,2,WisF);
1458     }
1459   }
1460   else {
1461     if (nr>nl) {
1462       // we have to shift quad on 1
1463       ShiftQuad(quad,1,WisF);
1464     }
1465     else {
1466       // we have to shift quad on 3
1467       ShiftQuad(quad,3,WisF);
1468     }
1469   }
1470
1471   nb = quad->side[0]->NbPoints();
1472   nr = quad->side[1]->NbPoints();
1473   nt = quad->side[2]->NbPoints();
1474   nl = quad->side[3]->NbPoints();
1475   dh = abs(nb-nt);
1476   dv = abs(nr-nl);
1477   int nbh  = Max(nb,nt);
1478   int nbv = Max(nr,nl);
1479   int addh = 0;
1480   int addv = 0;
1481
1482   // ----------- Old version ---------------
1483   // orientation of face and 3 main domain for future faces
1484   //       0   top    1
1485   //      1------------1
1486   //       |   |  |   |
1487   //       |   |  |   |
1488   //       | L |  | R |
1489   //  left |   |  |   | rigth
1490   //       |  /    \  |
1491   //       | /  C   \ |
1492   //       |/        \|
1493   //      0------------0
1494   //       0  bottom  1
1495
1496   // ----------- New version ---------------
1497   // orientation of face and 3 main domain for future faces
1498   //       0   top    1
1499   //      1------------1
1500   //       |  |____|  |
1501   //       |  /    \  |
1502   //       | /  C   \ |
1503   //  left |/________\| rigth
1504   //       |          |
1505   //       |          |
1506   //       |          |
1507   //      0------------0
1508   //       0  bottom  1
1509
1510   if (dh>dv) {
1511     addv = (dh-dv)/2;
1512     nbv = nbv + addv;
1513   }
1514   else { // dv>=dh
1515     addh = (dv-dh)/2;
1516     nbh = nbh + addh;
1517   }
1518
1519   const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
1520   const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
1521   const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1);
1522   const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
1523
1524   if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl)
1525     return error(COMPERR_BAD_INPUT_MESH);
1526
1527   if ( myNeedSmooth )
1528     UpdateDegenUV( quad );
1529
1530   // arrays for normalized params
1531   //cout<<"Dump B:"<<endl;
1532   TColStd_SequenceOfReal npb, npr, npt, npl;
1533   for (i=0; i<nb; i++) {
1534     npb.Append(uv_eb[i].normParam);
1535     //cout<<"i="<<i<<" par="<<uv_eb[i].normParam<<" npar="<<uv_eb[i].normParam;
1536     //const SMDS_MeshNode* N = uv_eb[i].node;
1537     //cout<<" node("<<N->X()<<","<<N->Y()<<","<<N->Z()<<")"<<endl;
1538   }
1539   for (i=0; i<nr; i++) {
1540     npr.Append(uv_er[i].normParam);
1541   }
1542   for (i=0; i<nt; i++) {
1543     npt.Append(uv_et[i].normParam);
1544   }
1545   for (i=0; i<nl; i++) {
1546     npl.Append(uv_el[i].normParam);
1547   }
1548
1549   int dl,dr;
1550   if (OldVersion) {
1551     // add some params to right and left after the first param
1552     // insert to right
1553     dr = nbv - nr;
1554     double dpr = (npr.Value(2) - npr.Value(1))/(dr+1);
1555     for (i=1; i<=dr; i++) {
1556       npr.InsertAfter(1,npr.Value(2)-dpr);
1557     }
1558     // insert to left
1559     dl = nbv - nl;
1560     dpr = (npl.Value(2) - npl.Value(1))/(dl+1);
1561     for (i=1; i<=dl; i++) {
1562       npl.InsertAfter(1,npl.Value(2)-dpr);
1563     }
1564   }
1565   //cout<<"npb:";
1566   //for (i=1; i<=npb.Length(); i++) {
1567   //  cout<<" "<<npb.Value(i);
1568   //}
1569   //cout<<endl;
1570
1571   gp_XY a0(uv_eb.front().u, uv_eb.front().v);
1572   gp_XY a1(uv_eb.back().u,  uv_eb.back().v);
1573   gp_XY a2(uv_et.back().u,  uv_et.back().v);
1574   gp_XY a3(uv_et.front().u, uv_et.front().v);
1575   //cout<<" a0("<<a0.X()<<","<<a0.Y()<<")"<<" a1("<<a1.X()<<","<<a1.Y()<<")"
1576   //    <<" a2("<<a2.X()<<","<<a2.Y()<<")"<<" a3("<<a3.X()<<","<<a3.Y()<<")"<<endl;
1577
1578   int nnn = Min(nr,nl);
1579   // auxilary sequence of XY for creation nodes
1580   // in the bottom part of central domain
1581   // Length of UVL and UVR must be == nbv-nnn
1582   TColgp_SequenceOfXY UVL, UVR, UVT;
1583
1584   if (OldVersion) {
1585     // step1: create faces for left domain
1586     StdMeshers_Array2OfNode NodesL(1,dl+1,1,nl);
1587     // add left nodes
1588     for (j=1; j<=nl; j++)
1589       NodesL.SetValue(1,j,uv_el[j-1].node);
1590     if (dl>0) {
1591       // add top nodes
1592       for (i=1; i<=dl; i++)
1593         NodesL.SetValue(i+1,nl,uv_et[i].node);
1594       // create and add needed nodes
1595       TColgp_SequenceOfXY UVtmp;
1596       for (i=1; i<=dl; i++) {
1597         double x0 = npt.Value(i+1);
1598         double x1 = x0;
1599         // diagonal node
1600         double y0 = npl.Value(i+1);
1601         double y1 = npr.Value(i+1);
1602         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1603         gp_Pnt P = S->Value(UV.X(),UV.Y());
1604         SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1605         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1606         NodesL.SetValue(i+1,1,N);
1607         if (UVL.Length()<nbv-nnn) UVL.Append(UV);
1608         // internal nodes
1609         for (j=2; j<nl; j++) {
1610           double y0 = npl.Value(dl+j);
1611           double y1 = npr.Value(dl+j);
1612           gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1613           gp_Pnt P = S->Value(UV.X(),UV.Y());
1614           SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1615           meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1616           NodesL.SetValue(i+1,j,N);
1617           if (i==dl) UVtmp.Append(UV);
1618         }
1619       }
1620       for (i=1; i<=UVtmp.Length() && UVL.Length()<nbv-nnn; i++) {
1621         UVL.Append(UVtmp.Value(i));
1622       }
1623       //cout<<"Dump NodesL:"<<endl;
1624       //for (i=1; i<=dl+1; i++) {
1625       //  cout<<"i="<<i;
1626       //  for (j=1; j<=nl; j++) {
1627       //    cout<<" ("<<NodesL.Value(i,j)->X()<<","<<NodesL.Value(i,j)->Y()<<","<<NodesL.Value(i,j)->Z()<<")";
1628       //  }
1629       //  cout<<endl;
1630       //}
1631       // create faces
1632       for (i=1; i<=dl; i++) {
1633         for (j=1; j<nl; j++) {
1634           if (WisF) {
1635             SMDS_MeshFace* F =
1636               myHelper->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j),
1637                               NodesL.Value(i+1,j+1), NodesL.Value(i,j+1));
1638             if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1639           }
1640           else {
1641             SMDS_MeshFace* F =
1642               myHelper->AddFace(NodesL.Value(i,j), NodesL.Value(i,j+1),
1643                               NodesL.Value(i+1,j+1), NodesL.Value(i+1,j));
1644             if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1645           }
1646         }
1647       }
1648     }
1649     else {
1650       // fill UVL using c2d
1651       for (i=1; i<npl.Length() && UVL.Length()<nbv-nnn; i++) {
1652         UVL.Append(gp_UV (uv_el[i].u, uv_el[i].v));
1653       }
1654     }
1655
1656     // step2: create faces for right domain
1657     StdMeshers_Array2OfNode NodesR(1,dr+1,1,nr);
1658     // add right nodes
1659     for (j=1; j<=nr; j++)
1660       NodesR.SetValue(1,j,uv_er[nr-j].node);
1661     if (dr>0) {
1662       // add top nodes
1663       for (i=1; i<=dr; i++)
1664         NodesR.SetValue(i+1,1,uv_et[nt-1-i].node);
1665       // create and add needed nodes
1666       TColgp_SequenceOfXY UVtmp;
1667       for (i=1; i<=dr; i++) {
1668         double x0 = npt.Value(nt-i);
1669         double x1 = x0;
1670         // diagonal node
1671         double y0 = npl.Value(i+1);
1672         double y1 = npr.Value(i+1);
1673         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1674         gp_Pnt P = S->Value(UV.X(),UV.Y());
1675         SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1676         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1677         NodesR.SetValue(i+1,nr,N);
1678         if (UVR.Length()<nbv-nnn) UVR.Append(UV);
1679         // internal nodes
1680         for (j=2; j<nr; j++) {
1681           double y0 = npl.Value(nbv-j+1);
1682           double y1 = npr.Value(nbv-j+1);
1683           gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1684           gp_Pnt P = S->Value(UV.X(),UV.Y());
1685           SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1686           meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1687           NodesR.SetValue(i+1,j,N);
1688           if (i==dr) UVtmp.Prepend(UV);
1689         }
1690       }
1691       for (i=1; i<=UVtmp.Length() && UVR.Length()<nbv-nnn; i++) {
1692         UVR.Append(UVtmp.Value(i));
1693       }
1694       // create faces
1695       for (i=1; i<=dr; i++) {
1696         for (j=1; j<nr; j++) {
1697           if (WisF) {
1698             SMDS_MeshFace* F =
1699               myHelper->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j),
1700                               NodesR.Value(i+1,j+1), NodesR.Value(i,j+1));
1701             if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1702           }
1703           else {
1704             SMDS_MeshFace* F =
1705               myHelper->AddFace(NodesR.Value(i,j), NodesR.Value(i,j+1),
1706                               NodesR.Value(i+1,j+1), NodesR.Value(i+1,j));
1707             if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1708           }
1709         }
1710       }
1711     }
1712     else {
1713       // fill UVR using c2d
1714       for (i=1; i<npr.Length() && UVR.Length()<nbv-nnn; i++) {
1715         UVR.Append(gp_UV(uv_er[i].u, uv_er[i].v));
1716       }
1717     }
1718
1719     // step3: create faces for central domain
1720     StdMeshers_Array2OfNode NodesC(1,nb,1,nbv);
1721     // add first line using NodesL
1722     for (i=1; i<=dl+1; i++)
1723       NodesC.SetValue(1,i,NodesL(i,1));
1724     for (i=2; i<=nl; i++)
1725       NodesC.SetValue(1,dl+i,NodesL(dl+1,i));
1726     // add last line using NodesR
1727     for (i=1; i<=dr+1; i++)
1728       NodesC.SetValue(nb,i,NodesR(i,nr));
1729     for (i=1; i<nr; i++)
1730       NodesC.SetValue(nb,dr+i+1,NodesR(dr+1,nr-i));
1731     // add top nodes (last columns)
1732     for (i=dl+2; i<nbh-dr; i++)
1733       NodesC.SetValue(i-dl,nbv,uv_et[i-1].node);
1734     // add bottom nodes (first columns)
1735     for (i=2; i<nb; i++)
1736       NodesC.SetValue(i,1,uv_eb[i-1].node);
1737
1738     // create and add needed nodes
1739     // add linear layers
1740     for (i=2; i<nb; i++) {
1741       double x0 = npt.Value(dl+i);
1742       double x1 = x0;
1743       for (j=1; j<nnn; j++) {
1744         double y0 = npl.Value(nbv-nnn+j);
1745         double y1 = npr.Value(nbv-nnn+j);
1746         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1747         gp_Pnt P = S->Value(UV.X(),UV.Y());
1748         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1749         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1750         NodesC.SetValue(i,nbv-nnn+j,N);
1751         if ( j==1 )
1752           UVT.Append( UV );
1753       }
1754     }
1755     // add diagonal layers
1756     //cout<<"UVL.Length()="<<UVL.Length()<<" UVR.Length()="<<UVR.Length()<<endl;
1757     //cout<<"Dump UVL:"<<endl;
1758     //for (i=1; i<=UVL.Length(); i++) {
1759     //  cout<<" ("<<UVL.Value(i).X()<<","<<UVL.Value(i).Y()<<")";
1760     //}
1761     //cout<<endl;
1762     gp_UV A2 = UVR.Value(nbv-nnn);
1763     gp_UV A3 = UVL.Value(nbv-nnn);
1764     for (i=1; i<nbv-nnn; i++) {
1765       gp_UV p1 = UVR.Value(i);
1766       gp_UV p3 = UVL.Value(i);
1767       double y = i / double(nbv-nnn);
1768       for (j=2; j<nb; j++) {
1769         double x = npb.Value(j);
1770         gp_UV p0( uv_eb[j-1].u, uv_eb[j-1].v );
1771         gp_UV p2 = UVT.Value( j-1 );
1772         gp_UV UV = CalcUV(x, y, a0, a1, A2, A3, p0,p1,p2,p3 );
1773         gp_Pnt P = S->Value(UV.X(),UV.Y());
1774         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1775         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(),UV.Y());
1776         NodesC.SetValue(j,i+1,N);
1777       }
1778     }
1779     // create faces
1780     for (i=1; i<nb; i++) {
1781       for (j=1; j<nbv; j++) {
1782         if (WisF) {
1783           SMDS_MeshFace* F =
1784             myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
1785                             NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
1786           if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1787         }
1788         else {
1789           SMDS_MeshFace* F =
1790             myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1),
1791                             NodesC.Value(i+1,j+1), NodesC.Value(i+1,j));
1792           if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1793         }
1794       }
1795     }
1796   }
1797
1798   else { // New version (!OldVersion)
1799     // step1: create faces for bottom rectangle domain
1800     StdMeshers_Array2OfNode NodesBRD(1,nb,1,nnn-1);
1801     // fill UVL and UVR using c2d
1802     for (j=0; j<nb; j++) {
1803       NodesBRD.SetValue(j+1,1,uv_eb[j].node);
1804     }
1805     for (i=1; i<nnn-1; i++) {
1806       NodesBRD.SetValue(1,i+1,uv_el[i].node);
1807       NodesBRD.SetValue(nb,i+1,uv_er[i].node);
1808       for (j=2; j<nb; j++) {
1809         double x = npb.Value(j);
1810         double y = (1-x) * npl.Value(i+1) + x * npr.Value(i+1);
1811         gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1812         gp_Pnt P = S->Value(UV.X(),UV.Y());
1813         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1814         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(),UV.Y());
1815         NodesBRD.SetValue(j,i+1,N);
1816       }
1817     }
1818     for (j=1; j<nnn-1; j++) {
1819       for (i=1; i<nb; i++) {
1820         if (WisF) {
1821           SMDS_MeshFace* F =
1822             myHelper->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i+1,j),
1823                             NodesBRD.Value(i+1,j+1), NodesBRD.Value(i,j+1));
1824           if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1825         }
1826         else {
1827           SMDS_MeshFace* F =
1828             myHelper->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i,j+1),
1829                             NodesBRD.Value(i+1,j+1), NodesBRD.Value(i+1,j));
1830           if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1831         }
1832       }
1833     }
1834     int drl = abs(nr-nl);
1835     // create faces for region C
1836     StdMeshers_Array2OfNode NodesC(1,nb,1,drl+1+addv);
1837     // add nodes from previous region
1838     for (j=1; j<=nb; j++) {
1839       NodesC.SetValue(j,1,NodesBRD.Value(j,nnn-1));
1840     }
1841     if ((drl+addv) > 0) {
1842       int n1,n2;
1843       if (nr>nl) {
1844         n1 = 1;
1845         n2 = drl + 1;
1846         TColgp_SequenceOfXY UVtmp;
1847         double drparam = npr.Value(nr) - npr.Value(nnn-1);
1848         double dlparam = npl.Value(nnn) - npl.Value(nnn-1);
1849         double y0,y1;
1850         for (i=1; i<=drl; i++) {
1851           // add existed nodes from right edge
1852           NodesC.SetValue(nb,i+1,uv_er[nnn+i-2].node);
1853           //double dtparam = npt.Value(i+1);
1854           y1 = npr.Value(nnn+i-1); // param on right edge
1855           double dpar = (y1 - npr.Value(nnn-1))/drparam;
1856           y0 = npl.Value(nnn-1) + dpar*dlparam; // param on left edge
1857           double dy = y1 - y0;
1858           for (j=1; j<nb; j++) {
1859             double x = npt.Value(i+1) + npb.Value(j)*(1-npt.Value(i+1));
1860             double y = y0 + dy*x;
1861             gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1862             gp_Pnt P = S->Value(UV.X(),UV.Y());
1863             SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1864             meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1865             NodesC.SetValue(j,i+1,N);
1866           }
1867         }
1868         double dy0 = (1-y0)/(addv+1);
1869         double dy1 = (1-y1)/(addv+1);
1870         for (i=1; i<=addv; i++) {
1871           double yy0 = y0 + dy0*i;
1872           double yy1 = y1 + dy1*i;
1873           double dyy = yy1 - yy0;
1874           for (j=1; j<=nb; j++) {
1875             double x = npt.Value(i+1+drl) +
1876               npb.Value(j) * (npt.Value(nt-i) - npt.Value(i+1+drl));
1877             double y = yy0 + dyy*x;
1878             gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1879             gp_Pnt P = S->Value(UV.X(),UV.Y());
1880             SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1881             meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1882             NodesC.SetValue(j,i+drl+1,N);
1883           }
1884         }
1885       }
1886       else { // nr<nl
1887         n2 = 1;
1888         n1 = drl + 1;
1889         TColgp_SequenceOfXY UVtmp;
1890         double dlparam = npl.Value(nl) - npl.Value(nnn-1);
1891         double drparam = npr.Value(nnn) - npr.Value(nnn-1);
1892         double y0 = npl.Value(nnn-1);
1893         double y1 = npr.Value(nnn-1);
1894         for (i=1; i<=drl; i++) {
1895           // add existed nodes from right edge
1896           NodesC.SetValue(1,i+1,uv_el[nnn+i-2].node);
1897           y0 = npl.Value(nnn+i-1); // param on left edge
1898           double dpar = (y0 - npl.Value(nnn-1))/dlparam;
1899           y1 = npr.Value(nnn-1) + dpar*drparam; // param on right edge
1900           double dy = y1 - y0;
1901           for (j=2; j<=nb; j++) {
1902             double x = npb.Value(j)*npt.Value(nt-i);
1903             double y = y0 + dy*x;
1904             gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1905             gp_Pnt P = S->Value(UV.X(),UV.Y());
1906             SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1907             meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1908             NodesC.SetValue(j,i+1,N);
1909           }
1910         }
1911         double dy0 = (1-y0)/(addv+1);
1912         double dy1 = (1-y1)/(addv+1);
1913         for (i=1; i<=addv; i++) {
1914           double yy0 = y0 + dy0*i;
1915           double yy1 = y1 + dy1*i;
1916           double dyy = yy1 - yy0;
1917           for (j=1; j<=nb; j++) {
1918             double x = npt.Value(i+1) +
1919               npb.Value(j) * (npt.Value(nt-i-drl) - npt.Value(i+1));
1920             double y = yy0 + dyy*x;
1921             gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1922             gp_Pnt P = S->Value(UV.X(),UV.Y());
1923             SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1924             meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1925             NodesC.SetValue(j,i+drl+1,N);
1926           }
1927         }
1928       }
1929       // create faces
1930       for (j=1; j<=drl+addv; j++) {
1931         for (i=1; i<nb; i++) {
1932           if (WisF) {
1933             SMDS_MeshFace* F =
1934               myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
1935                               NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
1936             if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1937           }
1938           else {
1939             SMDS_MeshFace* F =
1940               myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1),
1941                               NodesC.Value(i+1,j+1), NodesC.Value(i+1,j));
1942             if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1943           }
1944         }
1945       } // end nr<nl
1946
1947       StdMeshers_Array2OfNode NodesLast(1,nt,1,2);
1948       for (i=1; i<=nt; i++) {
1949         NodesLast.SetValue(i,2,uv_et[i-1].node);
1950       }
1951       int nnn=0;
1952       for (i=n1; i<drl+addv+1; i++) {
1953         nnn++;
1954         NodesLast.SetValue(nnn,1,NodesC.Value(1,i));
1955       }
1956       for (i=1; i<=nb; i++) {
1957         nnn++;
1958         NodesLast.SetValue(nnn,1,NodesC.Value(i,drl+addv+1));
1959       }
1960       for (i=drl+addv; i>=n2; i--) {
1961         nnn++;
1962         NodesLast.SetValue(nnn,1,NodesC.Value(nb,i));
1963       }
1964       for (i=1; i<nt; i++) {
1965         if (WisF) {
1966           SMDS_MeshFace* F =
1967             myHelper->AddFace(NodesLast.Value(i,1), NodesLast.Value(i+1,1),
1968                             NodesLast.Value(i+1,2), NodesLast.Value(i,2));
1969           if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1970         }
1971         else {
1972           SMDS_MeshFace* F =
1973             myHelper->AddFace(NodesLast.Value(i,1), NodesLast.Value(i,2),
1974                             NodesLast.Value(i+1,2), NodesLast.Value(i+1,2));
1975           if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1976         }
1977       }
1978     } // if ((drl+addv) > 0)
1979
1980   } // end new version implementation
1981
1982   bool isOk = true;
1983   return isOk;
1984 }
1985
1986
1987 //=======================================================================
1988 /*!
1989  * Evaluate only quandrangle faces
1990  */
1991 //=======================================================================
1992
1993 bool StdMeshers_Quadrangle_2D::EvaluateQuadPref(SMESH_Mesh &        aMesh,
1994                                                 const TopoDS_Shape& aShape,
1995                                                 std::vector<int>& aNbNodes,
1996                                                 MapShapeNbElems& aResMap,
1997                                                 bool IsQuadratic)
1998 {
1999   // Auxilary key in order to keep old variant
2000   // of meshing after implementation new variant
2001   // for bug 0016220 from Mantis.
2002   bool OldVersion = false;
2003   if (myQuadType == QUAD_QUADRANGLE_PREF_REVERSED)
2004     OldVersion = true;
2005
2006   const TopoDS_Face& F = TopoDS::Face(aShape);
2007   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
2008
2009   int nb = aNbNodes[0];
2010   int nr = aNbNodes[1];
2011   int nt = aNbNodes[2];
2012   int nl = aNbNodes[3];
2013   int dh = abs(nb-nt);
2014   int dv = abs(nr-nl);
2015
2016   if (dh>=dv) {
2017     if (nt>nb) {
2018       // it is a base case => not shift
2019     }
2020     else {
2021       // we have to shift on 2
2022       nb = aNbNodes[2];
2023       nr = aNbNodes[3];
2024       nt = aNbNodes[0];
2025       nl = aNbNodes[1];
2026     }
2027   }
2028   else {
2029     if (nr>nl) {
2030       // we have to shift quad on 1
2031       nb = aNbNodes[3];
2032       nr = aNbNodes[0];
2033       nt = aNbNodes[1];
2034       nl = aNbNodes[2];
2035     }
2036     else {
2037       // we have to shift quad on 3
2038       nb = aNbNodes[1];
2039       nr = aNbNodes[2];
2040       nt = aNbNodes[3];
2041       nl = aNbNodes[0];
2042     }
2043   }
2044
2045   dh = abs(nb-nt);
2046   dv = abs(nr-nl);
2047   int nbh  = Max(nb,nt);
2048   int nbv = Max(nr,nl);
2049   int addh = 0;
2050   int addv = 0;
2051
2052   if (dh>dv) {
2053     addv = (dh-dv)/2;
2054     nbv = nbv + addv;
2055   }
2056   else { // dv>=dh
2057     addh = (dv-dh)/2;
2058     nbh = nbh + addh;
2059   }
2060
2061   int dl,dr;
2062   if (OldVersion) {
2063     // add some params to right and left after the first param
2064     // insert to right
2065     dr = nbv - nr;
2066     // insert to left
2067     dl = nbv - nl;
2068   }
2069
2070   int nnn = Min(nr,nl);
2071
2072   int nbNodes = 0;
2073   int nbFaces = 0;
2074   if (OldVersion) {
2075     // step1: create faces for left domain
2076     if (dl>0) {
2077       nbNodes += dl*(nl-1);
2078       nbFaces += dl*(nl-1);
2079     }
2080     // step2: create faces for right domain
2081     if (dr>0) {
2082       nbNodes += dr*(nr-1);
2083       nbFaces += dr*(nr-1);
2084     }
2085     // step3: create faces for central domain
2086     nbNodes += (nb-2)*(nnn-1) + (nbv-nnn-1)*(nb-2);
2087     nbFaces += (nb-1)*(nbv-1);
2088   }
2089   else { // New version (!OldVersion)
2090     nbNodes += (nnn-2)*(nb-2);
2091     nbFaces += (nnn-2)*(nb-1);
2092     int drl = abs(nr-nl);
2093     nbNodes += drl*(nb-1) + addv*nb;
2094     nbFaces += (drl+addv)*(nb-1) + (nt-1);
2095   } // end new version implementation
2096
2097   std::vector<int> aVec(SMDSEntity_Last);
2098   for (int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i] = 0;
2099   if (IsQuadratic) {
2100     aVec[SMDSEntity_Quad_Quadrangle] = nbFaces;
2101     aVec[SMDSEntity_Node] = nbNodes + nbFaces*4;
2102     if (aNbNodes.size()==5) {
2103       aVec[SMDSEntity_Quad_Triangle] = aNbNodes[3] - 1;
2104       aVec[SMDSEntity_Quad_Quadrangle] = nbFaces - aNbNodes[3] + 1;
2105     }
2106   }
2107   else {
2108     aVec[SMDSEntity_Node] = nbNodes;
2109     aVec[SMDSEntity_Quadrangle] = nbFaces;
2110     if (aNbNodes.size()==5) {
2111       aVec[SMDSEntity_Triangle] = aNbNodes[3] - 1;
2112       aVec[SMDSEntity_Quadrangle] = nbFaces - aNbNodes[3] + 1;
2113     }
2114   }
2115   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
2116   aResMap.insert(std::make_pair(sm,aVec));
2117
2118   return true;
2119 }
2120
2121
2122 //=============================================================================
2123 /*! Split quadrangle in to 2 triangles by smallest diagonal
2124  *
2125  */
2126 //=============================================================================
2127 void StdMeshers_Quadrangle_2D::SplitQuad(SMESHDS_Mesh *theMeshDS,
2128                                          int theFaceID,
2129                                          const SMDS_MeshNode* theNode1,
2130                                          const SMDS_MeshNode* theNode2,
2131                                          const SMDS_MeshNode* theNode3,
2132                                          const SMDS_MeshNode* theNode4)
2133 {
2134   gp_Pnt a(theNode1->X(),theNode1->Y(),theNode1->Z());
2135   gp_Pnt b(theNode2->X(),theNode2->Y(),theNode2->Z());
2136   gp_Pnt c(theNode3->X(),theNode3->Y(),theNode3->Z());
2137   gp_Pnt d(theNode4->X(),theNode4->Y(),theNode4->Z());
2138   SMDS_MeshFace* face;
2139   if (a.Distance(c) > b.Distance(d)){
2140     face = myHelper->AddFace(theNode2, theNode4 , theNode1);
2141     if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
2142     face = myHelper->AddFace(theNode2, theNode3, theNode4);
2143     if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
2144
2145   }
2146   else{
2147     face = myHelper->AddFace(theNode1, theNode2 ,theNode3);
2148     if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
2149     face = myHelper->AddFace(theNode1, theNode3, theNode4);
2150     if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
2151   }
2152 }
2153
2154 //=======================================================================
2155 /*!
2156  *  Implementation of Reduced algorithm (meshing with quadrangles only)
2157  */
2158 //=======================================================================
2159 bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh &        aMesh,
2160                                                const TopoDS_Shape& aShape,
2161                                                FaceQuadStruct*     quad)
2162 {
2163   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
2164   const TopoDS_Face& F = TopoDS::Face(aShape);
2165   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
2166   int i,j,geomFaceID = meshDS->ShapeToIndex(F);
2167
2168   int nb = quad->side[0]->NbPoints();
2169   int nr = quad->side[1]->NbPoints();
2170   int nt = quad->side[2]->NbPoints();
2171   int nl = quad->side[3]->NbPoints();
2172
2173   //  Simple Reduce 8->6->4->2 (3 steps)      Multiple Reduce 8->2 (1 step)
2174   //
2175   //  .-----.-----.-----.-----.               .-----.-----.-----.-----.
2176   //  |    / \    |    / \    |               |    / \    |    / \    |
2177   //  |   /    .--.--.    \   |               |    / \    |    / \    |
2178   //  |   /   /   |   \   \   |               |   /  .----.----.  \   |
2179   //  .---.---.---.---.---.---.               |   / / \   |   / \ \   |
2180   //  |   /  / \  |  / \  \   |               |   / / \   |   / \ \   |
2181   //  |  /  /   .-.-.   \  \  |               |  / /  .---.---.  \ \  |
2182   //  |  /  /  /  |  \  \  \  |               |  / / / \  |  / \ \ \  |
2183   //  .--.--.--.--.--.--.--.--.               |  / / /  \ | /  \ \ \  |
2184   //  |  / /  / \ | / \  \ \  |               | / / /   .-.-.   \ \ \ |
2185   //  | / /  /  .-.-.  \  \ \ |               | / / /  /  |  \  \ \ \ |
2186   //  | / / /  /  |  \  \ \ \ |               | / / /  /  |  \  \ \ \ |
2187   //  .-.-.-.--.--.--.--.-.-.-.               .-.-.-.--.--.--.--.-.-.-.
2188
2189   bool MultipleReduce = false;
2190   {
2191     int nb1 = nb;
2192     int nr1 = nr;
2193     int nt1 = nt;
2194
2195     if (nr == nl) {
2196       if (nb < nt) {
2197         nt1 = nb;
2198         nb1 = nt;
2199       }
2200     }
2201     else if (nb == nt) {
2202       nr1 = nb; // and == nt
2203       if (nl < nr) {
2204         nt1 = nl;
2205         nb1 = nr;
2206       }
2207       else {
2208         nt1 = nr;
2209         nb1 = nl;
2210       }
2211     }
2212     else {
2213       return false;
2214     }
2215
2216     // number of rows and columns
2217     int nrows = nr1 - 1;
2218     int ncol_top = nt1 - 1;
2219     int ncol_bot = nb1 - 1;
2220     // maximum number of bottom elements for "tree" simple reduce 3->1
2221     int max_tree31 = ncol_top * pow(3.0, nrows);
2222     if (ncol_bot > max_tree31)
2223       MultipleReduce = true;
2224   }
2225
2226   if (MultipleReduce) { // == ComputeQuadPref QUAD_QUADRANGLE_PREF_REVERSED
2227     //==================================================
2228     int dh = abs(nb-nt);
2229     int dv = abs(nr-nl);
2230
2231     if (dh >= dv) {
2232       if (nt > nb) {
2233         // it is a base case => not shift quad but may be replacement is need
2234         ShiftQuad(quad,0,true);
2235       }
2236       else {
2237         // we have to shift quad on 2
2238         ShiftQuad(quad,2,true);
2239       }
2240     }
2241     else {
2242       if (nr > nl) {
2243         // we have to shift quad on 1
2244         ShiftQuad(quad,1,true);
2245       }
2246       else {
2247         // we have to shift quad on 3
2248         ShiftQuad(quad,3,true);
2249       }
2250     }
2251
2252     nb = quad->side[0]->NbPoints();
2253     nr = quad->side[1]->NbPoints();
2254     nt = quad->side[2]->NbPoints();
2255     nl = quad->side[3]->NbPoints();
2256     dh = abs(nb-nt);
2257     dv = abs(nr-nl);
2258     int nbh  = Max(nb,nt);
2259     int nbv = Max(nr,nl);
2260     int addh = 0;
2261     int addv = 0;
2262
2263     if (dh>dv) {
2264       addv = (dh-dv)/2;
2265       nbv = nbv + addv;
2266     }
2267     else { // dv>=dh
2268       addh = (dv-dh)/2;
2269       nbh = nbh + addh;
2270     }
2271
2272     const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
2273     const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
2274     const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1);
2275     const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
2276
2277     if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl)
2278       return error(COMPERR_BAD_INPUT_MESH);
2279
2280     if ( myNeedSmooth )
2281       UpdateDegenUV( quad );
2282
2283     // arrays for normalized params
2284     TColStd_SequenceOfReal npb, npr, npt, npl;
2285     for (j = 0; j < nb; j++) {
2286       npb.Append(uv_eb[j].normParam);
2287     }
2288     for (i = 0; i < nr; i++) {
2289       npr.Append(uv_er[i].normParam);
2290     }
2291     for (j = 0; j < nt; j++) {
2292       npt.Append(uv_et[j].normParam);
2293     }
2294     for (i = 0; i < nl; i++) {
2295       npl.Append(uv_el[i].normParam);
2296     }
2297
2298     int dl,dr;
2299     // orientation of face and 3 main domain for future faces
2300     //       0   top    1
2301     //      1------------1
2302     //       |   |  |   |
2303     //       |   |  |   |
2304     //       | L |  | R |
2305     //  left |   |  |   | rigth
2306     //       |  /    \  |
2307     //       | /  C   \ |
2308     //       |/        \|
2309     //      0------------0
2310     //       0  bottom  1
2311
2312     // add some params to right and left after the first param
2313     // insert to right
2314     dr = nbv - nr;
2315     double dpr = (npr.Value(2) - npr.Value(1))/(dr+1);
2316     for (i=1; i<=dr; i++) {
2317       npr.InsertAfter(1,npr.Value(2)-dpr);
2318     }
2319     // insert to left
2320     dl = nbv - nl;
2321     dpr = (npl.Value(2) - npl.Value(1))/(dl+1);
2322     for (i=1; i<=dl; i++) {
2323       npl.InsertAfter(1,npl.Value(2)-dpr);
2324     }
2325
2326     gp_XY a0 (uv_eb.front().u, uv_eb.front().v);
2327     gp_XY a1 (uv_eb.back().u,  uv_eb.back().v);
2328     gp_XY a2 (uv_et.back().u,  uv_et.back().v);
2329     gp_XY a3 (uv_et.front().u, uv_et.front().v);
2330
2331     int nnn = Min(nr,nl);
2332     // auxilary sequence of XY for creation nodes
2333     // in the bottom part of central domain
2334     // it's length must be == nbv-nnn-1
2335     TColgp_SequenceOfXY UVL;
2336     TColgp_SequenceOfXY UVR;
2337     //==================================================
2338
2339     // step1: create faces for left domain
2340     StdMeshers_Array2OfNode NodesL(1,dl+1,1,nl);
2341     // add left nodes
2342     for (j=1; j<=nl; j++)
2343       NodesL.SetValue(1,j,uv_el[j-1].node);
2344     if (dl>0) {
2345       // add top nodes
2346       for (i=1; i<=dl; i++)
2347         NodesL.SetValue(i+1,nl,uv_et[i].node);
2348       // create and add needed nodes
2349       TColgp_SequenceOfXY UVtmp;
2350       for (i=1; i<=dl; i++) {
2351         double x0 = npt.Value(i+1);
2352         double x1 = x0;
2353         // diagonal node
2354         double y0 = npl.Value(i+1);
2355         double y1 = npr.Value(i+1);
2356         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
2357         gp_Pnt P = S->Value(UV.X(),UV.Y());
2358         SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2359         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
2360         NodesL.SetValue(i+1,1,N);
2361         if (UVL.Length()<nbv-nnn-1) UVL.Append(UV);
2362         // internal nodes
2363         for (j=2; j<nl; j++) {
2364           double y0 = npl.Value(dl+j);
2365           double y1 = npr.Value(dl+j);
2366           gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
2367           gp_Pnt P = S->Value(UV.X(),UV.Y());
2368           SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2369           meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
2370           NodesL.SetValue(i+1,j,N);
2371           if (i==dl) UVtmp.Append(UV);
2372         }
2373       }
2374       for (i=1; i<=UVtmp.Length() && UVL.Length()<nbv-nnn-1; i++) {
2375         UVL.Append(UVtmp.Value(i));
2376       }
2377       // create faces
2378       for (i=1; i<=dl; i++) {
2379         for (j=1; j<nl; j++) {
2380             SMDS_MeshFace* F =
2381               myHelper->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j),
2382                               NodesL.Value(i+1,j+1), NodesL.Value(i,j+1));
2383             if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
2384         }
2385       }
2386     }
2387     else {
2388       // fill UVL using c2d
2389       for (i=1; i<npl.Length() && UVL.Length()<nbv-nnn-1; i++) {
2390         UVL.Append(gp_UV (uv_el[i].u, uv_el[i].v));
2391       }
2392     }
2393
2394     // step2: create faces for right domain
2395     StdMeshers_Array2OfNode NodesR(1,dr+1,1,nr);
2396     // add right nodes
2397     for (j=1; j<=nr; j++)
2398       NodesR.SetValue(1,j,uv_er[nr-j].node);
2399     if (dr>0) {
2400       // add top nodes
2401       for (i=1; i<=dr; i++)
2402         NodesR.SetValue(i+1,1,uv_et[nt-1-i].node);
2403       // create and add needed nodes
2404       TColgp_SequenceOfXY UVtmp;
2405       for (i=1; i<=dr; i++) {
2406         double x0 = npt.Value(nt-i);
2407         double x1 = x0;
2408         // diagonal node
2409         double y0 = npl.Value(i+1);
2410         double y1 = npr.Value(i+1);
2411         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
2412         gp_Pnt P = S->Value(UV.X(),UV.Y());
2413         SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2414         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
2415         NodesR.SetValue(i+1,nr,N);
2416         if (UVR.Length()<nbv-nnn-1) UVR.Append(UV);
2417         // internal nodes
2418         for (j=2; j<nr; j++) {
2419           double y0 = npl.Value(nbv-j+1);
2420           double y1 = npr.Value(nbv-j+1);
2421           gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
2422           gp_Pnt P = S->Value(UV.X(),UV.Y());
2423           SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2424           meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
2425           NodesR.SetValue(i+1,j,N);
2426           if (i==dr) UVtmp.Prepend(UV);
2427         }
2428       }
2429       for (i=1; i<=UVtmp.Length() && UVR.Length()<nbv-nnn-1; i++) {
2430         UVR.Append(UVtmp.Value(i));
2431       }
2432       // create faces
2433       for (i=1; i<=dr; i++) {
2434         for (j=1; j<nr; j++) {
2435             SMDS_MeshFace* F =
2436               myHelper->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j),
2437                               NodesR.Value(i+1,j+1), NodesR.Value(i,j+1));
2438             if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
2439         }
2440       }
2441     }
2442     else {
2443       // fill UVR using c2d
2444       for (i=1; i<npr.Length() && UVR.Length()<nbv-nnn-1; i++) {
2445         UVR.Append(gp_UV(uv_er[i].u, uv_er[i].v));
2446       }
2447     }
2448
2449     // step3: create faces for central domain
2450     StdMeshers_Array2OfNode NodesC(1,nb,1,nbv);
2451     // add first line using NodesL
2452     for (i=1; i<=dl+1; i++)
2453       NodesC.SetValue(1,i,NodesL(i,1));
2454     for (i=2; i<=nl; i++)
2455       NodesC.SetValue(1,dl+i,NodesL(dl+1,i));
2456     // add last line using NodesR
2457     for (i=1; i<=dr+1; i++)
2458       NodesC.SetValue(nb,i,NodesR(i,nr));
2459     for (i=1; i<nr; i++)
2460       NodesC.SetValue(nb,dr+i+1,NodesR(dr+1,nr-i));
2461     // add top nodes (last columns)
2462     for (i=dl+2; i<nbh-dr; i++)
2463       NodesC.SetValue(i-dl,nbv,uv_et[i-1].node);
2464     // add bottom nodes (first columns)
2465     for (i=2; i<nb; i++)
2466       NodesC.SetValue(i,1,uv_eb[i-1].node);
2467
2468     // create and add needed nodes
2469     // add linear layers
2470     for (i=2; i<nb; i++) {
2471       double x0 = npt.Value(dl+i);
2472       double x1 = x0;
2473       for (j=1; j<nnn; j++) {
2474         double y0 = npl.Value(nbv-nnn+j);
2475         double y1 = npr.Value(nbv-nnn+j);
2476         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
2477         gp_Pnt P = S->Value(UV.X(),UV.Y());
2478         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2479         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
2480         NodesC.SetValue(i,nbv-nnn+j,N);
2481       }
2482     }
2483     // add diagonal layers
2484     for (i=1; i<nbv-nnn; i++) {
2485       double du = UVR.Value(i).X() - UVL.Value(i).X();
2486       double dv = UVR.Value(i).Y() - UVL.Value(i).Y();
2487       for (j=2; j<nb; j++) {
2488         double u = UVL.Value(i).X() + du*npb.Value(j);
2489         double v = UVL.Value(i).Y() + dv*npb.Value(j);
2490         gp_Pnt P = S->Value(u,v);
2491         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2492         meshDS->SetNodeOnFace(N, geomFaceID, u, v);
2493         NodesC.SetValue(j,i+1,N);
2494       }
2495     }
2496     // create faces
2497     for (i=1; i<nb; i++) {
2498       for (j=1; j<nbv; j++) {
2499           SMDS_MeshFace* F =
2500             myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
2501                             NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
2502           if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
2503       }
2504     }
2505     // TODO ???
2506   } // end Multiple Reduce implementation
2507   else { // Simple Reduce (!MultipleReduce)
2508     //=========================================================
2509     if (nr == nl) {
2510       if (nt < nb) {
2511         // it is a base case => not shift quad
2512         //ShiftQuad(quad,0,true);
2513       }
2514       else {
2515         // we have to shift quad on 2
2516         ShiftQuad(quad,2,true);
2517       }
2518     }
2519     else {
2520       if (nl > nr) {
2521         // we have to shift quad on 1
2522         ShiftQuad(quad,1,true);
2523       }
2524       else {
2525         // we have to shift quad on 3
2526         ShiftQuad(quad,3,true);
2527       }
2528     }
2529
2530     nb = quad->side[0]->NbPoints();
2531     nr = quad->side[1]->NbPoints();
2532     nt = quad->side[2]->NbPoints();
2533     nl = quad->side[3]->NbPoints();
2534
2535     // number of rows and columns
2536     int nrows = nr - 1; // and also == nl - 1
2537     int ncol_top = nt - 1;
2538     int ncol_bot = nb - 1;
2539     int npair_top = ncol_top / 2;
2540     // maximum number of bottom elements for "linear" simple reduce 4->2
2541     int max_lin = ncol_top + npair_top * 2 * nrows;
2542     // maximum number of bottom elements for "linear" simple reduce 4->2
2543     int max_lin31 = ncol_top + ncol_top * 2 * nrows;
2544     // maximum number of bottom elements for "tree" simple reduce 4->2
2545     int max_tree42 = npair_top * pow(2.0, nrows + 1);
2546     if (ncol_top > npair_top * 2) {
2547       int delta = ncol_bot - max_tree42;
2548       for (int irow = 1; irow < nrows; irow++) {
2549         int nfour = delta / 4;
2550         delta -= nfour * 2;
2551       }
2552       if (delta <= (ncol_top - npair_top * 2))
2553         max_tree42 = ncol_bot;
2554     }
2555     // maximum number of bottom elements for "tree" simple reduce 3->1
2556     //int max_tree31 = ncol_top * pow(3.0, nrows);
2557     bool is_lin_31 = false;
2558     bool is_lin_42 = false;
2559     bool is_tree_31 = false;
2560     bool is_tree_42 = false;
2561     if (ncol_bot > max_lin) {
2562       if (ncol_bot <= max_lin31) {
2563         is_lin_31 = true;
2564         max_lin = max_lin31;
2565       }
2566     }
2567     else {
2568       // if ncol_bot is a 3*n or not 2*n
2569       if ((ncol_bot/3)*3 == ncol_bot || (ncol_bot/2)*2 != ncol_bot) {
2570         is_lin_31 = true;
2571         max_lin = max_lin31;
2572       }
2573       else {
2574         is_lin_42 = true;
2575       }
2576     }
2577     if (ncol_bot > max_lin) { // not "linear"
2578       is_tree_31 = (ncol_bot > max_tree42);
2579       if (ncol_bot <= max_tree42) {
2580         if ((ncol_bot/3)*3 == ncol_bot || (ncol_bot/2)*2 != ncol_bot) {
2581           is_tree_31 = true;
2582         }
2583         else {
2584           is_tree_42 = true;
2585         }
2586       }
2587     }
2588
2589     const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
2590     const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
2591     const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1);
2592     const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
2593
2594     if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl)
2595       return error(COMPERR_BAD_INPUT_MESH);
2596
2597     // arrays for normalized params
2598     TColStd_SequenceOfReal npb, npr, npt, npl;
2599     for (j = 0; j < nb; j++) {
2600       npb.Append(uv_eb[j].normParam);
2601     }
2602     for (i = 0; i < nr; i++) {
2603       npr.Append(uv_er[i].normParam);
2604     }
2605     for (j = 0; j < nt; j++) {
2606       npt.Append(uv_et[j].normParam);
2607     }
2608     for (i = 0; i < nl; i++) {
2609       npl.Append(uv_el[i].normParam);
2610     }
2611
2612     // We will ajust new points to this grid
2613     if (!SetNormalizedGrid(aMesh, aShape, quad))
2614       return false;
2615
2616     // TODO ???
2617     gp_XY a0 (uv_eb.front().u, uv_eb.front().v);
2618     gp_XY a1 (uv_eb.back().u,  uv_eb.back().v);
2619     gp_XY a2 (uv_et.back().u,  uv_et.back().v);
2620     gp_XY a3 (uv_et.front().u, uv_et.front().v);
2621     //=========================================================
2622
2623     TColStd_SequenceOfInteger curr_base, next_base;
2624     TColStd_SequenceOfReal curr_par_u, curr_par_v;
2625     TColStd_SequenceOfReal next_par_u, next_par_v;
2626     StdMeshers_Array2OfNode NodesBRD (1,nb, 1,nr);
2627     for (j = 1; j <= nb; j++) {
2628       NodesBRD.SetValue(j, 1, uv_eb[j - 1].node); // bottom
2629       curr_base.Append(j);
2630       next_base.Append(-1);
2631       curr_par_u.Append(uv_eb[j-1].u);
2632       curr_par_v.Append(uv_eb[j-1].v);
2633       next_par_u.Append(0.);
2634       next_par_v.Append(0.);
2635     }
2636     for (j = 1; j <= nt; j++) {
2637       NodesBRD.SetValue(j, nr, uv_et[j - 1].node); // top
2638     }
2639
2640     int curr_base_len = nb;
2641     int next_base_len = 0;
2642
2643     if (is_tree_42) {
2644       // "tree" simple reduce "42": 2->4->8->16->32->...
2645       //
2646       //  .-------------------------------.-------------------------------. nr
2647       //  |    \                          |                          /    |
2648       //  |         \     .---------------.---------------.     /         |
2649       //  |               |               |               |               |
2650       //  .---------------.---------------.---------------.---------------.
2651       //  | \             |             / | \             |             / |
2652       //  |     \ .-------.-------. /     |     \ .-------.-------. /     |
2653       //  |       |       |       |       |       |       |       |       |
2654       //  .-------.-------.-------.-------.-------.-------.-------.-------. i
2655       //  |\      |      /|\      |      /|\      |      /|\      |      /|
2656       //  |  \.---.---./  |  \.---.---./  |  \.---.---./  |  \.---.---./  |
2657       //  |   |   |   |   |   |   |   |   |   |   |   |   |   |   |   |   |
2658       //  .---.---.---.---.---.---.---.---.---.---.---.---.---.---.---.---.
2659       //  |\  |  /|\  |  /|\  |  /|\  |  /|\  |  /|\  |  /|\  |  /|\  |  /|
2660       //  | .-.-. | .-.-. | .-.-. | .-.-. | .-.-. | .-.-. | .-.-. | .-.-. |
2661       //  | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | |
2662       //  .-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. 1
2663       //  1                               j                               nb
2664
2665       for (i = 1; i < nr; i++) { // layer by layer
2666         // left
2667         NodesBRD.SetValue(1, i+1, uv_el[i].node);
2668         next_base.SetValue(++next_base_len, 1);
2669         // right
2670         NodesBRD.SetValue(nb, i+1, uv_er[i].node);
2671
2672         next_par_u.SetValue(next_base_len, uv_el[i].u);
2673         next_par_v.SetValue(next_base_len, uv_el[i].v);
2674
2675         // to stop reducing, if number of nodes reaches nt
2676         int delta = curr_base_len - nt;
2677
2678         //double du = uv_er[i].u - uv_el[i].u;
2679         //double dv = uv_er[i].v - uv_el[i].v;
2680
2681         // to calculate normalized parameter, we must know number of points in next layer
2682         int nb_four = (curr_base_len - 1) / 4;
2683         int nb_next = nb_four*2 + (curr_base_len - nb_four*4);
2684         if (nb_next < nt) nb_next = nt;
2685
2686         for (j = 1; j + 4 <= curr_base_len && delta > 0; j += 4, delta -= 2) {
2687           // add one "HH": nodes a,b,c,d,e and faces 1,2,3,4,5,6
2688           //
2689           //  .-----a-----b i + 1
2690           //  |\ 5  | 6  /|
2691           //  | \   |   / |
2692           //  |  c--d--e  |
2693           //  |1 |2 |3 |4 |
2694           //  |  |  |  |  |
2695           //  .--.--.--.--. i
2696           //
2697           //  j     j+2   j+4
2698
2699           double u,v;
2700
2701           // a (i + 1, j + 2)
2702           const SMDS_MeshNode* Na;
2703           next_base_len++;
2704           next_base.SetValue(next_base_len, curr_base.Value(j + 2));
2705           if (i + 1 == nr) { // top
2706             Na = uv_et[next_base_len - 1].node;
2707             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na);
2708             u = uv_et[next_base_len - 1].u;
2709             v = uv_et[next_base_len - 1].v;
2710           }
2711           else {
2712             //double norm_par = double(next_base_len - 1)/double(nb_next - 1);
2713             //u = uv_el[i].u + du * norm_par;
2714             //v = uv_el[i].v + dv * norm_par;
2715             {
2716               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
2717               int nearest_node_j = (int)rel;
2718               rel -= nearest_node_j;
2719               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
2720               double u1 = quad->uv_grid[ij].u;
2721               double v1 = quad->uv_grid[ij].v;
2722               double u2 = quad->uv_grid[ij + 1].u;
2723               double v2 = quad->uv_grid[ij + 1].v;
2724               double duj = (u2 - u1) * rel;
2725               double dvj = (v2 - v1) * rel;
2726               u = u1 + duj;
2727               v = v1 + dvj;
2728             }
2729             //u = uv_el[i].u + du*npb.Value(curr_base.Value(j + 2));
2730             //v = uv_el[i].v + dv*npb.Value(curr_base.Value(j + 2));
2731             gp_Pnt P = S->Value(u,v);
2732             SMDS_MeshNode* Na1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
2733             meshDS->SetNodeOnFace(Na1, geomFaceID, u, v);
2734             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na1);
2735             Na = Na1;
2736           }
2737           next_par_u.SetValue(next_base_len, u);
2738           next_par_v.SetValue(next_base_len, v);
2739
2740           // b (i + 1, j + 4)
2741           const SMDS_MeshNode* Nb;
2742           next_base_len++;
2743           next_base.SetValue(next_base_len, curr_base.Value(j + 4));
2744           if (i + 1 == nr) { // top
2745             Nb = uv_et[next_base_len - 1].node;
2746             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb);
2747             u = uv_et[next_base_len - 1].u;
2748             v = uv_et[next_base_len - 1].v;
2749           }
2750           else if (j + 4 == curr_base_len) { // right
2751             Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
2752             u = uv_er[i].u;
2753             v = uv_er[i].v;
2754           }
2755           else {
2756             //double norm_par = double(next_base_len - 1)/double(nb_next - 1);
2757             //u = uv_el[i].u + du * norm_par;
2758             //v = uv_el[i].v + dv * norm_par;
2759             {
2760               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
2761               int nearest_node_j = (int)rel;
2762               rel -= nearest_node_j;
2763               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
2764               double u1 = quad->uv_grid[ij].u;
2765               double v1 = quad->uv_grid[ij].v;
2766               double u2 = quad->uv_grid[ij + 1].u;
2767               double v2 = quad->uv_grid[ij + 1].v;
2768               double duj = (u2 - u1) * rel;
2769               double dvj = (v2 - v1) * rel;
2770               u = u1 + duj;
2771               v = v1 + dvj;
2772             }
2773             //u = uv_el[i].u + du*npb.Value(curr_base.Value(j + 4));
2774             //v = uv_el[i].v + dv*npb.Value(curr_base.Value(j + 4));
2775             gp_Pnt P = S->Value(u,v);
2776             SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
2777             meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v);
2778             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1);
2779             Nb = Nb1;
2780           }
2781           next_par_u.SetValue(next_base_len, u);
2782           next_par_v.SetValue(next_base_len, v);
2783
2784           // c
2785           u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 2)) / 2.0;
2786           v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 2)) / 2.0;
2787           gp_Pnt P = S->Value(u,v);
2788           SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z());
2789           meshDS->SetNodeOnFace(Nc, geomFaceID, u, v);
2790
2791           // d
2792           u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 1)) / 2.0;
2793           v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 1)) / 2.0;
2794           P = S->Value(u,v);
2795           SMDS_MeshNode* Nd = meshDS->AddNode(P.X(), P.Y(), P.Z());
2796           meshDS->SetNodeOnFace(Nd, geomFaceID, u, v);
2797
2798           // e
2799           u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len)) / 2.0;
2800           v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len)) / 2.0;
2801           P = S->Value(u,v);
2802           SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z());
2803           meshDS->SetNodeOnFace(Ne, geomFaceID, u, v);
2804
2805           // Faces
2806           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i),
2807                                               NodesBRD.Value(curr_base.Value(j + 1), i),
2808                                               Nc,
2809                                               NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1));
2810           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
2811
2812           SMDS_MeshFace* F2 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i),
2813                                               NodesBRD.Value(curr_base.Value(j + 2), i),
2814                                               Nd, Nc);
2815           if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID);
2816
2817           SMDS_MeshFace* F3 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i),
2818                                               NodesBRD.Value(curr_base.Value(j + 3), i),
2819                                               Ne, Nd);
2820           if (F3) meshDS->SetMeshElementOnShape(F3, geomFaceID);
2821
2822           SMDS_MeshFace* F4 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 3), i),
2823                                               NodesBRD.Value(curr_base.Value(j + 4), i),
2824                                               Nb, Ne);
2825           if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID);
2826
2827           SMDS_MeshFace* F5 = myHelper->AddFace(Nc, Nd, Na,
2828                                               NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1));
2829           if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID);
2830
2831           SMDS_MeshFace* F6 = myHelper->AddFace(Nd, Ne, Nb, Na);
2832           if (F6) meshDS->SetMeshElementOnShape(F6, geomFaceID);
2833         }
2834
2835         // not reduced side elements (if any)
2836         for (; j < curr_base_len; j++) {
2837           // f (i + 1, j + 1)
2838           const SMDS_MeshNode* Nf;
2839           double u,v;
2840           next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
2841           if (i + 1 == nr) { // top
2842             Nf = uv_et[next_base_len - 1].node;
2843             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
2844             u = uv_et[next_base_len - 1].u;
2845             v = uv_et[next_base_len - 1].v;
2846           }
2847           else if (j + 1 == curr_base_len) { // right
2848             Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
2849             u = uv_er[i].u;
2850             v = uv_er[i].v;
2851           }
2852           else {
2853             //double norm_par = double(next_base_len - 1)/double(nb_next - 1);
2854             //u = uv_el[i].u + du * norm_par;
2855             //v = uv_el[i].v + dv * norm_par;
2856             {
2857               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
2858               int nearest_node_j = (int)rel;
2859               rel -= nearest_node_j;
2860               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
2861               double u1 = quad->uv_grid[ij].u;
2862               double v1 = quad->uv_grid[ij].v;
2863               double u2 = quad->uv_grid[ij + 1].u;
2864               double v2 = quad->uv_grid[ij + 1].v;
2865               double duj = (u2 - u1) * rel;
2866               double dvj = (v2 - v1) * rel;
2867               u = u1 + duj;
2868               v = v1 + dvj;
2869             }
2870             //u = uv_el[i].u + du*npb.Value(curr_base.Value(j + 1));
2871             //v = uv_el[i].v + dv*npb.Value(curr_base.Value(j + 1));
2872             gp_Pnt P = S->Value(u,v);
2873             SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
2874             meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
2875             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
2876             Nf = Nf1;
2877           }
2878           next_par_u.SetValue(next_base_len, u);
2879           next_par_v.SetValue(next_base_len, v);
2880           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
2881                                               NodesBRD.Value(curr_base.Value(j + 1), i),
2882                                               NodesBRD.Value(next_base.Value(next_base_len), i + 1),
2883                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
2884           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
2885         }
2886
2887         curr_base_len = next_base_len;
2888         curr_base = next_base;
2889         curr_par_u = next_par_u;
2890         curr_par_v = next_par_v;
2891         next_base_len = 0;
2892       }
2893     } // end "tree" simple reduce "42"
2894     else if (is_tree_31) {
2895       // "tree" simple reduce "31": 1->3->9->27->...
2896       //
2897       //  .-----------------------------------------------------. nr
2898       //  |        \                                   /        |
2899       //  |                 .-----------------.                 |
2900       //  |                 |                 |                 |
2901       //  .-----------------.-----------------.-----------------.
2902       //  |   \         /   |   \         /   |   \         /   |
2903       //  |     .-----.     |     .-----.     |     .-----.     | i
2904       //  |     |     |     |     |     |     |     |     |     |
2905       //  .-----.-----.-----.-----.-----.-----.-----.-----.-----.
2906       //  |\   /|\   /|\   /|\   /|\   /|\   /|\   /|\   /|\   /|
2907       //  | .-. | .-. | .-. | .-. | .-. | .-. | .-. | .-. | .-. |
2908       //  | | | | | | | | | | | | | | | | | | | | | | | | | | | |
2909       //  .-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. 1
2910       //  1                          j                          nb
2911
2912       for (i = 1; i < nr; i++) { // layer by layer
2913         // left
2914         NodesBRD.SetValue(1, i+1, uv_el[i].node);
2915         next_base.SetValue(++next_base_len, 1);
2916         // right
2917         NodesBRD.SetValue(nb, i+1, uv_er[i].node);
2918
2919         next_par_u.SetValue(next_base_len, uv_el[i].u);
2920         next_par_v.SetValue(next_base_len, uv_el[i].v);
2921
2922         // to stop reducing, if number of nodes reaches nt
2923         int delta = curr_base_len - nt;
2924
2925         // to calculate normalized parameter, we must know number of points in next layer
2926         int nb_three = (curr_base_len - 1) / 3;
2927         int nb_next = nb_three + (curr_base_len - nb_three*3);
2928         if (nb_next < nt) nb_next = nt;
2929
2930         for (j = 1; j + 3 <= curr_base_len && delta > 0; j += 3, delta -= 2) {
2931           // add one "H": nodes b,c,e and faces 1,2,4,5
2932           //
2933           //  .---------b i + 1
2934           //  |\   5   /|
2935           //  | \     / |
2936           //  |  c---e  |
2937           //  |1 |2  |4 |
2938           //  |  |   |  |
2939           //  .--.---.--. i
2940           //
2941           //  j j+1 j+2 j+3
2942
2943           double u,v;
2944
2945           // b (i + 1, j + 3)
2946           const SMDS_MeshNode* Nb;
2947           next_base_len++;
2948           next_base.SetValue(next_base_len, curr_base.Value(j + 3));
2949           if (i + 1 == nr) { // top
2950             Nb = uv_et[next_base_len - 1].node;
2951             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb);
2952             u = uv_et[next_base_len - 1].u;
2953             v = uv_et[next_base_len - 1].v;
2954           }
2955           else if (j + 3 == curr_base_len) { // right
2956             Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
2957             u = uv_er[i].u;
2958             v = uv_er[i].v;
2959           }
2960           else {
2961             {
2962               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
2963               int nearest_node_j = (int)rel;
2964               rel -= nearest_node_j;
2965               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
2966               double u1 = quad->uv_grid[ij].u;
2967               double v1 = quad->uv_grid[ij].v;
2968               double u2 = quad->uv_grid[ij + 1].u;
2969               double v2 = quad->uv_grid[ij + 1].v;
2970               double duj = (u2 - u1) * rel;
2971               double dvj = (v2 - v1) * rel;
2972               u = u1 + duj;
2973               v = v1 + dvj;
2974             }
2975             gp_Pnt P = S->Value(u,v);
2976             SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
2977             meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v);
2978             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1);
2979             Nb = Nb1;
2980           }
2981           next_par_u.SetValue(next_base_len, u);
2982           next_par_v.SetValue(next_base_len, v);
2983
2984           // c and e
2985           double u1 = (curr_par_u.Value(j) + next_par_u.Value(next_base_len - 1)) / 2.0;
2986           double u2 = (curr_par_u.Value(j + 3) + next_par_u.Value(next_base_len)) / 2.0;
2987           double u3 = (u2 - u1) / 3.0;
2988
2989           double v1 = (curr_par_v.Value(j) + next_par_v.Value(next_base_len - 1)) / 2.0;
2990           double v2 = (curr_par_v.Value(j + 3) + next_par_v.Value(next_base_len)) / 2.0;
2991           double v3 = (v2 - v1) / 3.0;
2992
2993           // c
2994           u = u1 + u3;
2995           v = v1 + v3;
2996           gp_Pnt P = S->Value(u,v);
2997           SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z());
2998           meshDS->SetNodeOnFace(Nc, geomFaceID, u, v);
2999
3000           // e
3001           u = u1 + u3 + u3;
3002           v = v1 + v3 + v3;
3003           P = S->Value(u,v);
3004           SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z());
3005           meshDS->SetNodeOnFace(Ne, geomFaceID, u, v);
3006
3007           // Faces
3008           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i),
3009                                               NodesBRD.Value(curr_base.Value(j + 1), i),
3010                                               Nc,
3011                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3012           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3013
3014           SMDS_MeshFace* F2 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i),
3015                                               NodesBRD.Value(curr_base.Value(j + 2), i),
3016                                               Ne, Nc);
3017           if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID);
3018
3019           SMDS_MeshFace* F4 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i),
3020                                               NodesBRD.Value(curr_base.Value(j + 3), i),
3021                                               Nb, Ne);
3022           if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID);
3023
3024           SMDS_MeshFace* F5 = myHelper->AddFace(Nc, Ne, Nb,
3025                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3026           if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID);
3027         }
3028
3029         // not reduced side elements (if any)
3030         for (; j < curr_base_len; j++) {
3031           // f (i + 1, j + 1)
3032           const SMDS_MeshNode* Nf;
3033           double u,v;
3034           next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
3035           if (i + 1 == nr) { // top
3036             Nf = uv_et[next_base_len - 1].node;
3037             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
3038             u = uv_et[next_base_len - 1].u;
3039             v = uv_et[next_base_len - 1].v;
3040           }
3041           else if (j + 1 == curr_base_len) { // right
3042             Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
3043             u = uv_er[i].u;
3044             v = uv_er[i].v;
3045           }
3046           else {
3047             {
3048               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3049               int nearest_node_j = (int)rel;
3050               rel -= nearest_node_j;
3051               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3052               double u1 = quad->uv_grid[ij].u;
3053               double v1 = quad->uv_grid[ij].v;
3054               double u2 = quad->uv_grid[ij + 1].u;
3055               double v2 = quad->uv_grid[ij + 1].v;
3056               double duj = (u2 - u1) * rel;
3057               double dvj = (v2 - v1) * rel;
3058               u = u1 + duj;
3059               v = v1 + dvj;
3060             }
3061             gp_Pnt P = S->Value(u,v);
3062             SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3063             meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
3064             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
3065             Nf = Nf1;
3066           }
3067           next_par_u.SetValue(next_base_len, u);
3068           next_par_v.SetValue(next_base_len, v);
3069           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
3070                                               NodesBRD.Value(curr_base.Value(j + 1), i),
3071                                               NodesBRD.Value(next_base.Value(next_base_len), i + 1),
3072                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3073           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3074         }
3075
3076         curr_base_len = next_base_len;
3077         curr_base = next_base;
3078         curr_par_u = next_par_u;
3079         curr_par_v = next_par_v;
3080         next_base_len = 0;
3081       }
3082     } // end "tree" simple reduce "31"
3083     else if (is_lin_42) {
3084       // "linear" simple reduce "42": 4->8->12->16
3085       //
3086       //  .---------------.---------------.---------------.---------------. nr
3087       //  | \             |             / | \             |             / |
3088       //  |     \ .-------.-------. /     |     \ .-------.-------. /     |
3089       //  |       |       |       |       |       |       |       |       |
3090       //  .-------.-------.-------.-------.-------.-------.-------.-------.
3091       //  |      / \      |      / \      |      / \      |      / \      |
3092       //  |     /   \.----.----./   \     |     /   \.----.----./   \     | i
3093       //  |     /    |    |    |    \     |     /    |    |    |    \     |
3094       //  .-----.----.----.----.----.-----.-----.----.----.----.----.-----.
3095       //  |     /   / \   |  /  \   \     |     /   / \   |  /  \   \     |
3096       //  |    /   /    .-.-.    \   \    |    /   /    .-.-.    \   \    |
3097       //  |   /   /    /  |  \    \   \   |   /   /    /  |  \    \   \   |
3098       //  .---.---.---.---.---.---.---.---.---.---.---.---.---.---.---.---. 1
3099       //  1                               j                               nb
3100
3101       // nt = 5, nb = 7, nr = 4
3102       //int delta_all = 2;
3103       //int delta_one_col = 6;
3104       //int nb_col = 0;
3105       //int remainder = 2;
3106       //if (remainder > 0) nb_col++;
3107       //nb_col = 1;
3108       //int free_left = 1;
3109       //free_left += 2;
3110       //int free_middle = 4;
3111
3112       int delta_all = nb - nt;
3113       int delta_one_col = (nr - 1) * 2;
3114       int nb_col = delta_all / delta_one_col;
3115       int remainder = delta_all - nb_col * delta_one_col;
3116       if (remainder > 0) {
3117         nb_col++;
3118       }
3119       int free_left = ((nt - 1) - nb_col * 2) / 2;
3120       free_left += nr - 2;
3121       int free_middle = (nr - 2) * 2;
3122       if (remainder > 0 && nb_col == 1) {
3123         int nb_rows_short_col = remainder / 2;
3124         int nb_rows_thrown = (nr - 1) - nb_rows_short_col;
3125         free_left -= nb_rows_thrown;
3126       }
3127
3128       // nt = 5, nb = 17, nr = 4
3129       //int delta_all = 12;
3130       //int delta_one_col = 6;
3131       //int nb_col = 2;
3132       //int remainder = 0;
3133       //int free_left = 2;
3134       //int free_middle = 4;
3135
3136       for (i = 1; i < nr; i++, free_middle -= 2, free_left -= 1) { // layer by layer
3137         // left
3138         NodesBRD.SetValue(1, i+1, uv_el[i].node);
3139         next_base.SetValue(++next_base_len, 1);
3140         // right
3141         NodesBRD.SetValue(nb, i+1, uv_er[i].node);
3142
3143         // left
3144         next_par_u.SetValue(next_base_len, uv_el[i].u);
3145         next_par_v.SetValue(next_base_len, uv_el[i].v);
3146
3147         // to calculate normalized parameter, we must know number of points in next layer
3148         int nb_next = curr_base_len - nb_col * 2;
3149         if (remainder > 0 && i > remainder / 2)
3150           // take into account short "column"
3151           nb_next += 2;
3152         if (nb_next < nt) nb_next = nt;
3153
3154         // not reduced left elements
3155         for (j = 1; j <= free_left; j++) {
3156           // f (i + 1, j + 1)
3157           const SMDS_MeshNode* Nf;
3158           double u,v;
3159           next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
3160           if (i + 1 == nr) { // top
3161             Nf = uv_et[next_base_len - 1].node;
3162             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
3163             u = uv_et[next_base_len - 1].u;
3164             v = uv_et[next_base_len - 1].v;
3165           }
3166           else {
3167             {
3168               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3169               int nearest_node_j = (int)rel;
3170               rel -= nearest_node_j;
3171               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3172               double u1 = quad->uv_grid[ij].u;
3173               double v1 = quad->uv_grid[ij].v;
3174               double u2 = quad->uv_grid[ij + 1].u;
3175               double v2 = quad->uv_grid[ij + 1].v;
3176               double duj = (u2 - u1) * rel;
3177               double dvj = (v2 - v1) * rel;
3178               u = u1 + duj;
3179               v = v1 + dvj;
3180             }
3181             gp_Pnt P = S->Value(u,v);
3182             SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3183             meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
3184             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
3185             Nf = Nf1;
3186           }
3187           next_par_u.SetValue(next_base_len, u);
3188           next_par_v.SetValue(next_base_len, v);
3189           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
3190                                               NodesBRD.Value(curr_base.Value(j + 1), i),
3191                                               NodesBRD.Value(next_base.Value(next_base_len), i + 1),
3192                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3193           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3194         }
3195
3196         for (int icol = 1; icol <= nb_col; icol++) {
3197
3198           if (remainder > 0 && icol == nb_col && i > remainder / 2)
3199             // stop short "column"
3200             break;
3201
3202           // add one "HH": nodes a,b,c,d,e and faces 1,2,3,4,5,6
3203           //
3204           //  .-----a-----b i + 1
3205           //  |\ 5  | 6  /|
3206           //  | \   |   / |
3207           //  |  c--d--e  |
3208           //  |1 |2 |3 |4 |
3209           //  |  |  |  |  |
3210           //  .--.--.--.--. i
3211           //
3212           //  j     j+2   j+4
3213
3214           double u,v;
3215
3216           // a (i + 1, j + 2)
3217           const SMDS_MeshNode* Na;
3218           next_base_len++;
3219           next_base.SetValue(next_base_len, curr_base.Value(j + 2));
3220           if (i + 1 == nr) { // top
3221             Na = uv_et[next_base_len - 1].node;
3222             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na);
3223             u = uv_et[next_base_len - 1].u;
3224             v = uv_et[next_base_len - 1].v;
3225           }
3226           else {
3227             {
3228               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3229               int nearest_node_j = (int)rel;
3230               rel -= nearest_node_j;
3231               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3232               double u1 = quad->uv_grid[ij].u;
3233               double v1 = quad->uv_grid[ij].v;
3234               double u2 = quad->uv_grid[ij + 1].u;
3235               double v2 = quad->uv_grid[ij + 1].v;
3236               double duj = (u2 - u1) * rel;
3237               double dvj = (v2 - v1) * rel;
3238               u = u1 + duj;
3239               v = v1 + dvj;
3240             }
3241             gp_Pnt P = S->Value(u,v);
3242             SMDS_MeshNode* Na1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3243             meshDS->SetNodeOnFace(Na1, geomFaceID, u, v);
3244             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na1);
3245             Na = Na1;
3246           }
3247           next_par_u.SetValue(next_base_len, u);
3248           next_par_v.SetValue(next_base_len, v);
3249
3250           // b (i + 1, j + 4)
3251           const SMDS_MeshNode* Nb;
3252           next_base_len++;
3253           next_base.SetValue(next_base_len, curr_base.Value(j + 4));
3254           if (i + 1 == nr) { // top
3255             Nb = uv_et[next_base_len - 1].node;
3256             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb);
3257             u = uv_et[next_base_len - 1].u;
3258             v = uv_et[next_base_len - 1].v;
3259           }
3260           else if (j + 4 == curr_base_len) { // right
3261             Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
3262             u = uv_er[i].u;
3263             v = uv_er[i].v;
3264           }
3265           else {
3266             {
3267               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3268               int nearest_node_j = (int)rel;
3269               rel -= nearest_node_j;
3270               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3271               double u1 = quad->uv_grid[ij].u;
3272               double v1 = quad->uv_grid[ij].v;
3273               double u2 = quad->uv_grid[ij + 1].u;
3274               double v2 = quad->uv_grid[ij + 1].v;
3275               double duj = (u2 - u1) * rel;
3276               double dvj = (v2 - v1) * rel;
3277               u = u1 + duj;
3278               v = v1 + dvj;
3279             }
3280             gp_Pnt P = S->Value(u,v);
3281             SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3282             meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v);
3283             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1);
3284             Nb = Nb1;
3285           }
3286           next_par_u.SetValue(next_base_len, u);
3287           next_par_v.SetValue(next_base_len, v);
3288
3289           // c
3290           u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 2)) / 2.0;
3291           v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 2)) / 2.0;
3292           gp_Pnt P = S->Value(u,v);
3293           SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z());
3294           meshDS->SetNodeOnFace(Nc, geomFaceID, u, v);
3295
3296           // d
3297           u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 1)) / 2.0;
3298           v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 1)) / 2.0;
3299           P = S->Value(u,v);
3300           SMDS_MeshNode* Nd = meshDS->AddNode(P.X(), P.Y(), P.Z());
3301           meshDS->SetNodeOnFace(Nd, geomFaceID, u, v);
3302
3303           // e
3304           u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len)) / 2.0;
3305           v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len)) / 2.0;
3306           P = S->Value(u,v);
3307           SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z());
3308           meshDS->SetNodeOnFace(Ne, geomFaceID, u, v);
3309
3310           // Faces
3311           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i),
3312                                               NodesBRD.Value(curr_base.Value(j + 1), i),
3313                                               Nc,
3314                                               NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1));
3315           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3316
3317           SMDS_MeshFace* F2 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i),
3318                                               NodesBRD.Value(curr_base.Value(j + 2), i),
3319                                               Nd, Nc);
3320           if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID);
3321
3322           SMDS_MeshFace* F3 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i),
3323                                               NodesBRD.Value(curr_base.Value(j + 3), i),
3324                                               Ne, Nd);
3325           if (F3) meshDS->SetMeshElementOnShape(F3, geomFaceID);
3326
3327           SMDS_MeshFace* F4 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 3), i),
3328                                               NodesBRD.Value(curr_base.Value(j + 4), i),
3329                                               Nb, Ne);
3330           if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID);
3331
3332           SMDS_MeshFace* F5 = myHelper->AddFace(Nc, Nd, Na,
3333                                               NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1));
3334           if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID);
3335
3336           SMDS_MeshFace* F6 = myHelper->AddFace(Nd, Ne, Nb, Na);
3337           if (F6) meshDS->SetMeshElementOnShape(F6, geomFaceID);
3338
3339           j += 4;
3340
3341           // not reduced middle elements
3342           if (icol < nb_col) {
3343             if (remainder > 0 && icol == nb_col - 1 && i > remainder / 2)
3344               // pass middle elements before stopped short "column"
3345               break;
3346
3347             int free_add = free_middle;
3348             if (remainder > 0 && icol == nb_col - 1)
3349               // next "column" is short
3350               free_add -= (nr - 1) - (remainder / 2);
3351
3352             for (int imiddle = 1; imiddle <= free_add; imiddle++) {
3353               // f (i + 1, j + imiddle)
3354               const SMDS_MeshNode* Nf;
3355               double u,v;
3356               next_base.SetValue(++next_base_len, curr_base.Value(j + imiddle));
3357               if (i + 1 == nr) { // top
3358                 Nf = uv_et[next_base_len - 1].node;
3359                 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
3360                 u = uv_et[next_base_len - 1].u;
3361                 v = uv_et[next_base_len - 1].v;
3362               }
3363               else if (j + imiddle == curr_base_len) { // right
3364                 Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
3365                 u = uv_er[i].u;
3366                 v = uv_er[i].v;
3367               }
3368               else {
3369                 {
3370                   double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3371                   int nearest_node_j = (int)rel;
3372                   rel -= nearest_node_j;
3373                   int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3374                   double u1 = quad->uv_grid[ij].u;
3375                   double v1 = quad->uv_grid[ij].v;
3376                   double u2 = quad->uv_grid[ij + 1].u;
3377                   double v2 = quad->uv_grid[ij + 1].v;
3378                   double duj = (u2 - u1) * rel;
3379                   double dvj = (v2 - v1) * rel;
3380                   u = u1 + duj;
3381                   v = v1 + dvj;
3382                 }
3383                 gp_Pnt P = S->Value(u,v);
3384                 SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3385                 meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
3386                 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
3387                 Nf = Nf1;
3388               }
3389               next_par_u.SetValue(next_base_len, u);
3390               next_par_v.SetValue(next_base_len, v);
3391               SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j - 1 + imiddle), i),
3392                                                   NodesBRD.Value(curr_base.Value(j + imiddle), i),
3393                                                   NodesBRD.Value(next_base.Value(next_base_len), i + 1),
3394                                                   NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3395               if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3396             }
3397             j += free_add;
3398           }
3399         }
3400
3401         // not reduced right elements
3402         for (; j < curr_base_len; j++) {
3403           // f (i + 1, j + 1)
3404           const SMDS_MeshNode* Nf;
3405           double u,v;
3406           next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
3407           if (i + 1 == nr) { // top
3408             Nf = uv_et[next_base_len - 1].node;
3409             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
3410             u = uv_et[next_base_len - 1].u;
3411             v = uv_et[next_base_len - 1].v;
3412           }
3413           else if (j + 1 == curr_base_len) { // right
3414             Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
3415             u = uv_er[i].u;
3416             v = uv_er[i].v;
3417           }
3418           else {
3419             {
3420               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3421               int nearest_node_j = (int)rel;
3422               rel -= nearest_node_j;
3423               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3424               double u1 = quad->uv_grid[ij].u;
3425               double v1 = quad->uv_grid[ij].v;
3426               double u2 = quad->uv_grid[ij + 1].u;
3427               double v2 = quad->uv_grid[ij + 1].v;
3428               double duj = (u2 - u1) * rel;
3429               double dvj = (v2 - v1) * rel;
3430               u = u1 + duj;
3431               v = v1 + dvj;
3432             }
3433             gp_Pnt P = S->Value(u,v);
3434             SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3435             meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
3436             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
3437             Nf = Nf1;
3438           }
3439           next_par_u.SetValue(next_base_len, u);
3440           next_par_v.SetValue(next_base_len, v);
3441           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
3442                                               NodesBRD.Value(curr_base.Value(j + 1), i),
3443                                               NodesBRD.Value(next_base.Value(next_base_len), i + 1),
3444                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3445           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3446         }
3447
3448         curr_base_len = next_base_len;
3449         curr_base = next_base;
3450         curr_par_u = next_par_u;
3451         curr_par_v = next_par_v;
3452         next_base_len = 0;
3453       }
3454     } // end "linear" simple reduce "42"
3455     else if (is_lin_31) {
3456       // "linear" simple reduce "31": 2->6->10->14
3457       //
3458       //  .-----------------------------.-----------------------------. nr
3459       //  |     \                 /     |     \                 /     |
3460       //  |         .---------.         |         .---------.         |
3461       //  |         |         |         |         |         |         |
3462       //  .---------.---------.---------.---------.---------.---------.
3463       //  |        / \       / \        |        / \       / \        |
3464       //  |       /   .-----.   \       |       /   .-----.   \       | i
3465       //  |      /    |     |    \      |      /    |     |    \      |
3466       //  .-----.-----.-----.-----.-----.-----.-----.-----.-----.-----.
3467       //  |    /     / \   / \     \    |    /     / \   / \     \    |
3468       //  |   /     /   .-.   \     \   |   /     /   .-.   \     \   |
3469       //  |  /     /   /   \   \     \  |  /     /   /   \   \     \  |
3470       //  .--.----.---.-----.---.-----.-.--.----.---.-----.---.-----.-. 1
3471       //  1                             j                             nb
3472
3473       int delta_all = nb - nt;
3474       int delta_one_col = (nr - 1) * 2;
3475       int nb_col = delta_all / delta_one_col;
3476       int remainder = delta_all - nb_col * delta_one_col;
3477       if (remainder > 0) {
3478         nb_col++;
3479       }
3480       int free_left = ((nt - 1) - nb_col) / 2;
3481       free_left += nr - 2;
3482       int free_middle = (nr - 2) * 2;
3483       if (remainder > 0 && nb_col == 1) {
3484         int nb_rows_short_col = remainder / 2;
3485         int nb_rows_thrown = (nr - 1) - nb_rows_short_col;
3486         free_left -= nb_rows_thrown;
3487       }
3488
3489       for (i = 1; i < nr; i++, free_middle -= 2, free_left -= 1) { // layer by layer
3490         // left
3491         NodesBRD.SetValue(1, i+1, uv_el[i].node);
3492         next_base.SetValue(++next_base_len, 1);
3493         // right
3494         NodesBRD.SetValue(nb, i+1, uv_er[i].node);
3495
3496         // left
3497         next_par_u.SetValue(next_base_len, uv_el[i].u);
3498         next_par_v.SetValue(next_base_len, uv_el[i].v);
3499
3500         // to calculate normalized parameter, we must know number of points in next layer
3501         int nb_next = curr_base_len - nb_col * 2;
3502         if (remainder > 0 && i > remainder / 2)
3503           // take into account short "column"
3504           nb_next += 2;
3505         if (nb_next < nt) nb_next = nt;
3506
3507         // not reduced left elements
3508         for (j = 1; j <= free_left; j++) {
3509           // f (i + 1, j + 1)
3510           const SMDS_MeshNode* Nf;
3511           double u,v;
3512           next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
3513           if (i + 1 == nr) { // top
3514             Nf = uv_et[next_base_len - 1].node;
3515             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
3516             u = uv_et[next_base_len - 1].u;
3517             v = uv_et[next_base_len - 1].v;
3518           }
3519           else {
3520             {
3521               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3522               int nearest_node_j = (int)rel;
3523               rel -= nearest_node_j;
3524               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3525               double u1 = quad->uv_grid[ij].u;
3526               double v1 = quad->uv_grid[ij].v;
3527               double u2 = quad->uv_grid[ij + 1].u;
3528               double v2 = quad->uv_grid[ij + 1].v;
3529               double duj = (u2 - u1) * rel;
3530               double dvj = (v2 - v1) * rel;
3531               u = u1 + duj;
3532               v = v1 + dvj;
3533             }
3534             gp_Pnt P = S->Value(u,v);
3535             SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3536             meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
3537             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
3538             Nf = Nf1;
3539           }
3540           next_par_u.SetValue(next_base_len, u);
3541           next_par_v.SetValue(next_base_len, v);
3542           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
3543                                               NodesBRD.Value(curr_base.Value(j + 1), i),
3544                                               NodesBRD.Value(next_base.Value(next_base_len), i + 1),
3545                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3546           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3547         }
3548
3549         for (int icol = 1; icol <= nb_col; icol++) {
3550
3551           if (remainder > 0 && icol == nb_col && i > remainder / 2)
3552             // stop short "column"
3553             break;
3554
3555           // add one "H": nodes b,c,e and faces 1,2,4,5
3556           //
3557           //  .---------b i + 1
3558           //  |\   5   /|
3559           //  | \     / |
3560           //  |  c---e  |
3561           //  |1 |2  |4 |
3562           //  |  |   |  |
3563           //  .--.---.--. i
3564           //
3565           //  j j+1 j+2 j+3
3566
3567           double u,v;
3568
3569           // b (i + 1, j + 3)
3570           const SMDS_MeshNode* Nb;
3571           next_base_len++;
3572           next_base.SetValue(next_base_len, curr_base.Value(j + 3));
3573           if (i + 1 == nr) { // top
3574             Nb = uv_et[next_base_len - 1].node;
3575             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb);
3576             u = uv_et[next_base_len - 1].u;
3577             v = uv_et[next_base_len - 1].v;
3578           }
3579           else if (j + 3 == curr_base_len) { // right
3580             Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
3581             u = uv_er[i].u;
3582             v = uv_er[i].v;
3583           }
3584           else {
3585             {
3586               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3587               int nearest_node_j = (int)rel;
3588               rel -= nearest_node_j;
3589               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3590               double u1 = quad->uv_grid[ij].u;
3591               double v1 = quad->uv_grid[ij].v;
3592               double u2 = quad->uv_grid[ij + 1].u;
3593               double v2 = quad->uv_grid[ij + 1].v;
3594               double duj = (u2 - u1) * rel;
3595               double dvj = (v2 - v1) * rel;
3596               u = u1 + duj;
3597               v = v1 + dvj;
3598             }
3599             gp_Pnt P = S->Value(u,v);
3600             SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3601             meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v);
3602             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1);
3603             Nb = Nb1;
3604           }
3605           next_par_u.SetValue(next_base_len, u);
3606           next_par_v.SetValue(next_base_len, v);
3607
3608           // c and d
3609           double u1 = (curr_par_u.Value(j) + next_par_u.Value(next_base_len - 1)) / 2.0;
3610           double u2 = (curr_par_u.Value(j + 3) + next_par_u.Value(next_base_len)) / 2.0;
3611           double u3 = (u2 - u1) / 3.0;
3612
3613           double v1 = (curr_par_v.Value(j) + next_par_v.Value(next_base_len - 1)) / 2.0;
3614           double v2 = (curr_par_v.Value(j + 3) + next_par_v.Value(next_base_len)) / 2.0;
3615           double v3 = (v2 - v1) / 3.0;
3616
3617           // c
3618           u = u1 + u3;
3619           v = v1 + v3;
3620           gp_Pnt P = S->Value(u,v);
3621           SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z());
3622           meshDS->SetNodeOnFace(Nc, geomFaceID, u, v);
3623
3624           // e
3625           u = u1 + u3 + u3;
3626           v = v1 + v3 + v3;
3627           P = S->Value(u,v);
3628           SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z());
3629           meshDS->SetNodeOnFace(Ne, geomFaceID, u, v);
3630
3631           // Faces
3632           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i),
3633                                               NodesBRD.Value(curr_base.Value(j + 1), i),
3634                                               Nc,
3635                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3636           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3637
3638           SMDS_MeshFace* F2 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i),
3639                                               NodesBRD.Value(curr_base.Value(j + 2), i),
3640                                               Ne, Nc);
3641           if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID);
3642
3643           SMDS_MeshFace* F4 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i),
3644                                               NodesBRD.Value(curr_base.Value(j + 3), i),
3645                                               Nb, Ne);
3646           if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID);
3647
3648           SMDS_MeshFace* F5 = myHelper->AddFace(Nc, Ne, Nb,
3649                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3650           if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID);
3651
3652           j += 3;
3653
3654           // not reduced middle elements
3655           if (icol < nb_col) {
3656             if (remainder > 0 && icol == nb_col - 1 && i > remainder / 2)
3657               // pass middle elements before stopped short "column"
3658               break;
3659
3660             int free_add = free_middle;
3661             if (remainder > 0 && icol == nb_col - 1)
3662               // next "column" is short
3663               free_add -= (nr - 1) - (remainder / 2);
3664
3665             for (int imiddle = 1; imiddle <= free_add; imiddle++) {
3666               // f (i + 1, j + imiddle)
3667               const SMDS_MeshNode* Nf;
3668               double u,v;
3669               next_base.SetValue(++next_base_len, curr_base.Value(j + imiddle));
3670               if (i + 1 == nr) { // top
3671                 Nf = uv_et[next_base_len - 1].node;
3672                 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
3673                 u = uv_et[next_base_len - 1].u;
3674                 v = uv_et[next_base_len - 1].v;
3675               }
3676               else if (j + imiddle == curr_base_len) { // right
3677                 Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
3678                 u = uv_er[i].u;
3679                 v = uv_er[i].v;
3680               }
3681               else {
3682                 {
3683                   double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3684                   int nearest_node_j = (int)rel;
3685                   rel -= nearest_node_j;
3686                   int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3687                   double u1 = quad->uv_grid[ij].u;
3688                   double v1 = quad->uv_grid[ij].v;
3689                   double u2 = quad->uv_grid[ij + 1].u;
3690                   double v2 = quad->uv_grid[ij + 1].v;
3691                   double duj = (u2 - u1) * rel;
3692                   double dvj = (v2 - v1) * rel;
3693                   u = u1 + duj;
3694                   v = v1 + dvj;
3695                 }
3696                 gp_Pnt P = S->Value(u,v);
3697                 SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3698                 meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
3699                 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
3700                 Nf = Nf1;
3701               }
3702               next_par_u.SetValue(next_base_len, u);
3703               next_par_v.SetValue(next_base_len, v);
3704               SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j - 1 + imiddle), i),
3705                                                   NodesBRD.Value(curr_base.Value(j + imiddle), i),
3706                                                   NodesBRD.Value(next_base.Value(next_base_len), i + 1),
3707                                                   NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3708               if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3709             }
3710             j += free_add;
3711           }
3712         }
3713
3714         // not reduced right elements
3715         for (; j < curr_base_len; j++) {
3716           // f (i + 1, j + 1)
3717           const SMDS_MeshNode* Nf;
3718           double u,v;
3719           next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
3720           if (i + 1 == nr) { // top
3721             Nf = uv_et[next_base_len - 1].node;
3722             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
3723             u = uv_et[next_base_len - 1].u;
3724             v = uv_et[next_base_len - 1].v;
3725           }
3726           else if (j + 1 == curr_base_len) { // right
3727             Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
3728             u = uv_er[i].u;
3729             v = uv_er[i].v;
3730           }
3731           else {
3732             {
3733               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3734               int nearest_node_j = (int)rel;
3735               rel -= nearest_node_j;
3736               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3737               double u1 = quad->uv_grid[ij].u;
3738               double v1 = quad->uv_grid[ij].v;
3739               double u2 = quad->uv_grid[ij + 1].u;
3740               double v2 = quad->uv_grid[ij + 1].v;
3741               double duj = (u2 - u1) * rel;
3742               double dvj = (v2 - v1) * rel;
3743               u = u1 + duj;
3744               v = v1 + dvj;
3745             }
3746             gp_Pnt P = S->Value(u,v);
3747             SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3748             meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
3749             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
3750             Nf = Nf1;
3751           }
3752           next_par_u.SetValue(next_base_len, u);
3753           next_par_v.SetValue(next_base_len, v);
3754           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
3755                                                 NodesBRD.Value(curr_base.Value(j + 1), i),
3756                                                 NodesBRD.Value(next_base.Value(next_base_len), i + 1),
3757                                                 NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3758           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3759         }
3760
3761         curr_base_len = next_base_len;
3762         curr_base = next_base;
3763         curr_par_u = next_par_u;
3764         curr_par_v = next_par_v;
3765         next_base_len = 0;
3766       }
3767     } // end "linear" simple reduce "31"
3768     else {
3769     }
3770   } // end Simple Reduce implementation
3771
3772   bool isOk = true;
3773   return isOk;
3774 }
3775
3776 //================================================================================
3777 namespace // data for smoothing
3778 {
3779   struct TSmoothNode;
3780   // --------------------------------------------------------------------------------
3781   /*!
3782    * \brief Structure used to check validity of node position after smoothing.
3783    *        It holds two nodes connected to a smoothed node and belonging to
3784    *        one mesh face
3785    */
3786   struct TTriangle
3787   {
3788     TSmoothNode* _n1;
3789     TSmoothNode* _n2;
3790     TTriangle( TSmoothNode* n1=0, TSmoothNode* n2=0 ): _n1(n1), _n2(n2) {}
3791
3792     inline bool IsForward( gp_UV uv ) const;
3793   };
3794   // --------------------------------------------------------------------------------
3795   /*!
3796    * \brief Data of a smoothed node
3797    */
3798   struct TSmoothNode
3799   {
3800     gp_XY _uv;
3801     vector< TTriangle > _triangles; // if empty, then node is not movable
3802   };
3803   // --------------------------------------------------------------------------------
3804   inline bool TTriangle::IsForward( gp_UV uv ) const
3805   {
3806     gp_Vec2d v1( uv, _n1->_uv ), v2( uv, _n2->_uv );
3807     double d = v1 ^ v2;
3808     return d > 1e-100;
3809   }
3810 }
3811
3812 //================================================================================
3813 /*!
3814  * \brief Set UV of nodes on degenerated VERTEXes in the middle of degenerated EDGE
3815  *
3816  * WARNING: this method must be called AFTER retrieving UVPtStruct's from quad
3817  */
3818 //================================================================================
3819
3820 void StdMeshers_Quadrangle_2D::UpdateDegenUV(FaceQuadStruct* quad)
3821 {
3822   for ( unsigned i = 0; i < quad->side.size(); ++i )
3823   {
3824     StdMeshers_FaceSide* side = quad->side[i];
3825     const vector<UVPtStruct>& uvVec = side->GetUVPtStruct();
3826
3827     // find which end of the side is on degenerated shape
3828     int degenInd = -1;
3829     if ( myHelper->IsDegenShape( uvVec[0].node->getshapeId() ))
3830       degenInd = 0;
3831     else if ( myHelper->IsDegenShape( uvVec.back().node->getshapeId() ))
3832       degenInd = uvVec.size() - 1;
3833     else
3834       continue;
3835
3836     // find another side sharing the degenerated shape
3837     bool isPrev = ( degenInd == 0 );
3838     if ( i >= TOP_SIDE )
3839       isPrev = !isPrev;
3840     int i2 = ( isPrev ? ( i + 3 ) : ( i + 1 )) % 4;
3841     StdMeshers_FaceSide* side2 = quad->side[ i2 ];
3842     const vector<UVPtStruct>& uvVec2 = side2->GetUVPtStruct();
3843     int degenInd2 = -1;
3844     if ( uvVec[ degenInd ].node == uvVec2[0].node )
3845       degenInd2 = 0;
3846     else if ( uvVec[ degenInd ].node == uvVec2.back().node )
3847       degenInd2 = uvVec2.size() - 1;
3848     else
3849       throw SALOME_Exception( LOCALIZED( "Logical error" ));
3850
3851     // move UV in the middle
3852     uvPtStruct& uv1 = const_cast<uvPtStruct&>( uvVec [ degenInd  ]);
3853     uvPtStruct& uv2 = const_cast<uvPtStruct&>( uvVec2[ degenInd2 ]);
3854     uv1.u = uv2.u = 0.5 * ( uv1.u + uv2.u );
3855     uv1.v = uv2.v = 0.5 * ( uv1.v + uv2.v );
3856   }
3857 }
3858
3859 //================================================================================
3860 /*!
3861  * \brief Perform smoothing of 2D elements on a FACE with ignored degenerated EDGE
3862  */
3863 //================================================================================
3864
3865 void StdMeshers_Quadrangle_2D::Smooth (FaceQuadStruct* quad)
3866 {
3867   if ( !myNeedSmooth ) return;
3868
3869   // Get nodes to smooth
3870
3871   typedef map< const SMDS_MeshNode*, TSmoothNode, TIDCompare > TNo2SmooNoMap;
3872   TNo2SmooNoMap smooNoMap;
3873
3874   const TopoDS_Face& geomFace = TopoDS::Face( myHelper->GetSubShape() );
3875   SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
3876   SMESHDS_SubMesh* fSubMesh = meshDS->MeshElements( geomFace );
3877   SMDS_NodeIteratorPtr nIt = fSubMesh->GetNodes();
3878   while ( nIt->more() ) // loop on nodes bound to a FACE
3879   {
3880     const SMDS_MeshNode* node = nIt->next();
3881     TSmoothNode & sNode = smooNoMap[ node ];
3882     sNode._uv = myHelper->GetNodeUV( geomFace, node );
3883
3884     // set sNode._triangles
3885     SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator( SMDSAbs_Face );
3886     while ( fIt->more() )
3887     {
3888       const SMDS_MeshElement* face = fIt->next();
3889       const int nbN = face->NbCornerNodes();
3890       const int nInd = face->GetNodeIndex( node );
3891       const int prevInd = myHelper->WrapIndex( nInd - 1, nbN );
3892       const int nextInd = myHelper->WrapIndex( nInd + 1, nbN );
3893       const SMDS_MeshNode* prevNode = face->GetNode( prevInd );
3894       const SMDS_MeshNode* nextNode = face->GetNode( nextInd );
3895       sNode._triangles.push_back( TTriangle( & smooNoMap[ prevNode ],
3896                                              & smooNoMap[ nextNode ]));
3897     }
3898   }
3899   // set _uv of smooth nodes on FACE boundary
3900   for ( unsigned i = 0; i < quad->side.size(); ++i )
3901   {
3902     const vector<UVPtStruct>& uvVec = quad->side[i]->GetUVPtStruct();
3903     for ( unsigned j = 0; j < uvVec.size(); ++j )
3904     {
3905       TSmoothNode & sNode = smooNoMap[ uvVec[j].node ];
3906       sNode._uv.SetCoord( uvVec[j].u, uvVec[j].v );
3907     }
3908   }
3909
3910   // define refernce orientation in 2D
3911   TNo2SmooNoMap::iterator n2sn = smooNoMap.begin();
3912   for ( ; n2sn != smooNoMap.end(); ++n2sn )
3913     if ( !n2sn->second._triangles.empty() )
3914       break;
3915   if ( n2sn == smooNoMap.end() ) return;
3916   const TSmoothNode & sampleNode = n2sn->second;
3917   const bool refForward = ( sampleNode._triangles[0].IsForward( sampleNode._uv ));
3918
3919   // Smoothing
3920
3921   for ( int iLoop = 0; iLoop < 5; ++iLoop )
3922   {
3923     for ( n2sn = smooNoMap.begin(); n2sn != smooNoMap.end(); ++n2sn )
3924     {
3925       TSmoothNode& sNode = n2sn->second;
3926       if ( sNode._triangles.empty() )
3927         continue; // not movable node
3928
3929       // compute a new UV
3930       gp_XY newUV (0,0);
3931       for ( unsigned i = 0; i < sNode._triangles.size(); ++i )
3932         newUV += sNode._triangles[i]._n1->_uv;
3933       newUV /= sNode._triangles.size();
3934
3935       // check validity of the newUV
3936       bool isValid = true;
3937       for ( unsigned i = 0; i < sNode._triangles.size() && isValid; ++i )
3938         isValid = ( sNode._triangles[i].IsForward( newUV ) == refForward );
3939
3940       if ( isValid )
3941         sNode._uv = newUV;
3942     }
3943   }
3944
3945   // Set new XYZ to the smoothed nodes
3946
3947   Handle(Geom_Surface) surface = BRep_Tool::Surface( geomFace );
3948
3949   for ( n2sn = smooNoMap.begin(); n2sn != smooNoMap.end(); ++n2sn )
3950   {
3951     TSmoothNode& sNode = n2sn->second;
3952     if ( sNode._triangles.empty() )
3953       continue; // not movable node
3954
3955     SMDS_MeshNode* node = const_cast< SMDS_MeshNode*>( n2sn->first );
3956     gp_Pnt xyz = surface->Value( sNode._uv.X(), sNode._uv.Y() );
3957     meshDS->MoveNode( node, xyz.X(), xyz.Y(), xyz.Z() );
3958
3959     // store the new UV
3960     node->SetPosition( SMDS_PositionPtr( new SMDS_FacePosition( sNode._uv.X(), sNode._uv.Y() )));
3961   }
3962
3963   // Move medium nodes in quadratic mesh
3964   if ( _quadraticMesh )
3965   {
3966     const TLinkNodeMap& links = myHelper->GetTLinkNodeMap();
3967     TLinkNodeMap::const_iterator linkIt = links.begin();
3968     for ( ; linkIt != links.end(); ++linkIt )
3969     {
3970       const SMESH_TLink& link = linkIt->first;
3971       SMDS_MeshNode*     node = const_cast< SMDS_MeshNode*>( linkIt->second );
3972
3973       if ( node->getshapeId() != myHelper->GetSubShapeID() )
3974         continue; // medium node is on EDGE or VERTEX
3975
3976       gp_XY uv1 = myHelper->GetNodeUV( geomFace, link.node1(), node );
3977       gp_XY uv2 = myHelper->GetNodeUV( geomFace, link.node2(), node );
3978
3979       gp_XY uv  = myHelper->GetMiddleUV( surface, uv1, uv2 );
3980       node->SetPosition( SMDS_PositionPtr( new SMDS_FacePosition( uv.X(), uv.Y() )));
3981
3982       gp_Pnt xyz = surface->Value( uv.X(), uv.Y() );
3983       meshDS->MoveNode( node, xyz.X(), xyz.Y(), xyz.Z() );
3984     }
3985   }
3986 }