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     int nl1 = nl;
2195
2196     if (nr == nl) {
2197       if (nb < nt) {
2198         nt1 = nb;
2199         nb1 = nt;
2200       }
2201     }
2202     else if (nb == nt) {
2203       nl1 = nb; // and == nt
2204       nr1 = nb; // and == nt
2205       if (nl < nr) {
2206         nt1 = nl;
2207         nb1 = nr;
2208       }
2209       else {
2210         nt1 = nr;
2211         nb1 = nl;
2212       }
2213     }
2214     else {
2215       return false;
2216     }
2217
2218     // number of rows and columns
2219     int nrows = nr1 - 1; // and also == nl1 - 1
2220     int ncol_top = nt1 - 1;
2221     int ncol_bot = nb1 - 1;
2222     // maximum number of bottom elements for "tree" simple reduce 3->1
2223     int max_tree31 = ncol_top * pow(3.0, nrows);
2224     if (ncol_bot > max_tree31)
2225       MultipleReduce = true;
2226   }
2227
2228   if (MultipleReduce) { // == ComputeQuadPref QUAD_QUADRANGLE_PREF_REVERSED
2229     //==================================================
2230     int dh = abs(nb-nt);
2231     int dv = abs(nr-nl);
2232
2233     if (dh >= dv) {
2234       if (nt > nb) {
2235         // it is a base case => not shift quad but may be replacement is need
2236         ShiftQuad(quad,0,true);
2237       }
2238       else {
2239         // we have to shift quad on 2
2240         ShiftQuad(quad,2,true);
2241       }
2242     }
2243     else {
2244       if (nr > nl) {
2245         // we have to shift quad on 1
2246         ShiftQuad(quad,1,true);
2247       }
2248       else {
2249         // we have to shift quad on 3
2250         ShiftQuad(quad,3,true);
2251       }
2252     }
2253
2254     nb = quad->side[0]->NbPoints();
2255     nr = quad->side[1]->NbPoints();
2256     nt = quad->side[2]->NbPoints();
2257     nl = quad->side[3]->NbPoints();
2258     dh = abs(nb-nt);
2259     dv = abs(nr-nl);
2260     int nbh  = Max(nb,nt);
2261     int nbv = Max(nr,nl);
2262     int addh = 0;
2263     int addv = 0;
2264
2265     if (dh>dv) {
2266       addv = (dh-dv)/2;
2267       nbv = nbv + addv;
2268     }
2269     else { // dv>=dh
2270       addh = (dv-dh)/2;
2271       nbh = nbh + addh;
2272     }
2273
2274     const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
2275     const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
2276     const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1);
2277     const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
2278
2279     if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl)
2280       return error(COMPERR_BAD_INPUT_MESH);
2281
2282     if ( myNeedSmooth )
2283       UpdateDegenUV( quad );
2284
2285     // arrays for normalized params
2286     TColStd_SequenceOfReal npb, npr, npt, npl;
2287     for (j = 0; j < nb; j++) {
2288       npb.Append(uv_eb[j].normParam);
2289     }
2290     for (i = 0; i < nr; i++) {
2291       npr.Append(uv_er[i].normParam);
2292     }
2293     for (j = 0; j < nt; j++) {
2294       npt.Append(uv_et[j].normParam);
2295     }
2296     for (i = 0; i < nl; i++) {
2297       npl.Append(uv_el[i].normParam);
2298     }
2299
2300     int dl,dr;
2301     // orientation of face and 3 main domain for future faces
2302     //       0   top    1
2303     //      1------------1
2304     //       |   |  |   |
2305     //       |   |  |   |
2306     //       | L |  | R |
2307     //  left |   |  |   | rigth
2308     //       |  /    \  |
2309     //       | /  C   \ |
2310     //       |/        \|
2311     //      0------------0
2312     //       0  bottom  1
2313
2314     // add some params to right and left after the first param
2315     // insert to right
2316     dr = nbv - nr;
2317     double dpr = (npr.Value(2) - npr.Value(1))/(dr+1);
2318     for (i=1; i<=dr; i++) {
2319       npr.InsertAfter(1,npr.Value(2)-dpr);
2320     }
2321     // insert to left
2322     dl = nbv - nl;
2323     dpr = (npl.Value(2) - npl.Value(1))/(dl+1);
2324     for (i=1; i<=dl; i++) {
2325       npl.InsertAfter(1,npl.Value(2)-dpr);
2326     }
2327
2328     gp_XY a0 (uv_eb.front().u, uv_eb.front().v);
2329     gp_XY a1 (uv_eb.back().u,  uv_eb.back().v);
2330     gp_XY a2 (uv_et.back().u,  uv_et.back().v);
2331     gp_XY a3 (uv_et.front().u, uv_et.front().v);
2332
2333     int nnn = Min(nr,nl);
2334     // auxilary sequence of XY for creation nodes
2335     // in the bottom part of central domain
2336     // it's length must be == nbv-nnn-1
2337     TColgp_SequenceOfXY UVL;
2338     TColgp_SequenceOfXY UVR;
2339     //==================================================
2340
2341     // step1: create faces for left domain
2342     StdMeshers_Array2OfNode NodesL(1,dl+1,1,nl);
2343     // add left nodes
2344     for (j=1; j<=nl; j++)
2345       NodesL.SetValue(1,j,uv_el[j-1].node);
2346     if (dl>0) {
2347       // add top nodes
2348       for (i=1; i<=dl; i++)
2349         NodesL.SetValue(i+1,nl,uv_et[i].node);
2350       // create and add needed nodes
2351       TColgp_SequenceOfXY UVtmp;
2352       for (i=1; i<=dl; i++) {
2353         double x0 = npt.Value(i+1);
2354         double x1 = x0;
2355         // diagonal node
2356         double y0 = npl.Value(i+1);
2357         double y1 = npr.Value(i+1);
2358         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
2359         gp_Pnt P = S->Value(UV.X(),UV.Y());
2360         SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2361         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
2362         NodesL.SetValue(i+1,1,N);
2363         if (UVL.Length()<nbv-nnn-1) UVL.Append(UV);
2364         // internal nodes
2365         for (j=2; j<nl; j++) {
2366           double y0 = npl.Value(dl+j);
2367           double y1 = npr.Value(dl+j);
2368           gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
2369           gp_Pnt P = S->Value(UV.X(),UV.Y());
2370           SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2371           meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
2372           NodesL.SetValue(i+1,j,N);
2373           if (i==dl) UVtmp.Append(UV);
2374         }
2375       }
2376       for (i=1; i<=UVtmp.Length() && UVL.Length()<nbv-nnn-1; i++) {
2377         UVL.Append(UVtmp.Value(i));
2378       }
2379       // create faces
2380       for (i=1; i<=dl; i++) {
2381         for (j=1; j<nl; j++) {
2382             SMDS_MeshFace* F =
2383               myHelper->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j),
2384                               NodesL.Value(i+1,j+1), NodesL.Value(i,j+1));
2385             if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
2386         }
2387       }
2388     }
2389     else {
2390       // fill UVL using c2d
2391       for (i=1; i<npl.Length() && UVL.Length()<nbv-nnn-1; i++) {
2392         UVL.Append(gp_UV (uv_el[i].u, uv_el[i].v));
2393       }
2394     }
2395
2396     // step2: create faces for right domain
2397     StdMeshers_Array2OfNode NodesR(1,dr+1,1,nr);
2398     // add right nodes
2399     for (j=1; j<=nr; j++)
2400       NodesR.SetValue(1,j,uv_er[nr-j].node);
2401     if (dr>0) {
2402       // add top nodes
2403       for (i=1; i<=dr; i++)
2404         NodesR.SetValue(i+1,1,uv_et[nt-1-i].node);
2405       // create and add needed nodes
2406       TColgp_SequenceOfXY UVtmp;
2407       for (i=1; i<=dr; i++) {
2408         double x0 = npt.Value(nt-i);
2409         double x1 = x0;
2410         // diagonal node
2411         double y0 = npl.Value(i+1);
2412         double y1 = npr.Value(i+1);
2413         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
2414         gp_Pnt P = S->Value(UV.X(),UV.Y());
2415         SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2416         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
2417         NodesR.SetValue(i+1,nr,N);
2418         if (UVR.Length()<nbv-nnn-1) UVR.Append(UV);
2419         // internal nodes
2420         for (j=2; j<nr; j++) {
2421           double y0 = npl.Value(nbv-j+1);
2422           double y1 = npr.Value(nbv-j+1);
2423           gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
2424           gp_Pnt P = S->Value(UV.X(),UV.Y());
2425           SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2426           meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
2427           NodesR.SetValue(i+1,j,N);
2428           if (i==dr) UVtmp.Prepend(UV);
2429         }
2430       }
2431       for (i=1; i<=UVtmp.Length() && UVR.Length()<nbv-nnn-1; i++) {
2432         UVR.Append(UVtmp.Value(i));
2433       }
2434       // create faces
2435       for (i=1; i<=dr; i++) {
2436         for (j=1; j<nr; j++) {
2437             SMDS_MeshFace* F =
2438               myHelper->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j),
2439                               NodesR.Value(i+1,j+1), NodesR.Value(i,j+1));
2440             if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
2441         }
2442       }
2443     }
2444     else {
2445       // fill UVR using c2d
2446       for (i=1; i<npr.Length() && UVR.Length()<nbv-nnn-1; i++) {
2447         UVR.Append(gp_UV(uv_er[i].u, uv_er[i].v));
2448       }
2449     }
2450
2451     // step3: create faces for central domain
2452     StdMeshers_Array2OfNode NodesC(1,nb,1,nbv);
2453     // add first line using NodesL
2454     for (i=1; i<=dl+1; i++)
2455       NodesC.SetValue(1,i,NodesL(i,1));
2456     for (i=2; i<=nl; i++)
2457       NodesC.SetValue(1,dl+i,NodesL(dl+1,i));
2458     // add last line using NodesR
2459     for (i=1; i<=dr+1; i++)
2460       NodesC.SetValue(nb,i,NodesR(i,nr));
2461     for (i=1; i<nr; i++)
2462       NodesC.SetValue(nb,dr+i+1,NodesR(dr+1,nr-i));
2463     // add top nodes (last columns)
2464     for (i=dl+2; i<nbh-dr; i++)
2465       NodesC.SetValue(i-dl,nbv,uv_et[i-1].node);
2466     // add bottom nodes (first columns)
2467     for (i=2; i<nb; i++)
2468       NodesC.SetValue(i,1,uv_eb[i-1].node);
2469
2470     // create and add needed nodes
2471     // add linear layers
2472     for (i=2; i<nb; i++) {
2473       double x0 = npt.Value(dl+i);
2474       double x1 = x0;
2475       for (j=1; j<nnn; j++) {
2476         double y0 = npl.Value(nbv-nnn+j);
2477         double y1 = npr.Value(nbv-nnn+j);
2478         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
2479         gp_Pnt P = S->Value(UV.X(),UV.Y());
2480         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2481         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
2482         NodesC.SetValue(i,nbv-nnn+j,N);
2483       }
2484     }
2485     // add diagonal layers
2486     for (i=1; i<nbv-nnn; i++) {
2487       double du = UVR.Value(i).X() - UVL.Value(i).X();
2488       double dv = UVR.Value(i).Y() - UVL.Value(i).Y();
2489       for (j=2; j<nb; j++) {
2490         double u = UVL.Value(i).X() + du*npb.Value(j);
2491         double v = UVL.Value(i).Y() + dv*npb.Value(j);
2492         gp_Pnt P = S->Value(u,v);
2493         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2494         meshDS->SetNodeOnFace(N, geomFaceID, u, v);
2495         NodesC.SetValue(j,i+1,N);
2496       }
2497     }
2498     // create faces
2499     for (i=1; i<nb; i++) {
2500       for (j=1; j<nbv; j++) {
2501           SMDS_MeshFace* F =
2502             myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
2503                             NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
2504           if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
2505       }
2506     }
2507     // TODO ???
2508   } // end Multiple Reduce implementation
2509   else { // Simple Reduce (!MultipleReduce)
2510     //=========================================================
2511     if (nr == nl) {
2512       if (nt < nb) {
2513         // it is a base case => not shift quad
2514         //ShiftQuad(quad,0,true);
2515       }
2516       else {
2517         // we have to shift quad on 2
2518         ShiftQuad(quad,2,true);
2519       }
2520     }
2521     else {
2522       if (nl > nr) {
2523         // we have to shift quad on 1
2524         ShiftQuad(quad,1,true);
2525       }
2526       else {
2527         // we have to shift quad on 3
2528         ShiftQuad(quad,3,true);
2529       }
2530     }
2531
2532     nb = quad->side[0]->NbPoints();
2533     nr = quad->side[1]->NbPoints();
2534     nt = quad->side[2]->NbPoints();
2535     nl = quad->side[3]->NbPoints();
2536
2537     // number of rows and columns
2538     int nrows = nr - 1; // and also == nl - 1
2539     int ncol_top = nt - 1;
2540     int ncol_bot = nb - 1;
2541     int npair_top = ncol_top / 2;
2542     // maximum number of bottom elements for "linear" simple reduce 4->2
2543     int max_lin = ncol_top + npair_top * 2 * nrows;
2544     // maximum number of bottom elements for "linear" simple reduce 4->2
2545     int max_lin31 = ncol_top + ncol_top * 2 * nrows;
2546     // maximum number of bottom elements for "tree" simple reduce 4->2
2547     int max_tree42 = npair_top * pow(2.0, nrows + 1);
2548     if (ncol_top > npair_top * 2) {
2549       int delta = ncol_bot - max_tree42;
2550       for (int irow = 1; irow < nrows; irow++) {
2551         int nfour = delta / 4;
2552         delta -= nfour * 2;
2553       }
2554       if (delta <= (ncol_top - npair_top * 2))
2555         max_tree42 = ncol_bot;
2556     }
2557     // maximum number of bottom elements for "tree" simple reduce 3->1
2558     //int max_tree31 = ncol_top * pow(3.0, nrows);
2559     bool is_lin_31 = false;
2560     bool is_lin_42 = false;
2561     bool is_tree_31 = false;
2562     bool is_tree_42 = false;
2563     if (ncol_bot > max_lin) {
2564       if (ncol_bot <= max_lin31) {
2565         is_lin_31 = true;
2566         max_lin = max_lin31;
2567       }
2568     }
2569     else {
2570       // if ncol_bot is a 3*n or not 2*n
2571       if ((ncol_bot/3)*3 == ncol_bot || (ncol_bot/2)*2 != ncol_bot) {
2572         is_lin_31 = true;
2573         max_lin = max_lin31;
2574       }
2575       else {
2576         is_lin_42 = true;
2577       }
2578     }
2579     if (ncol_bot > max_lin) { // not "linear"
2580       is_tree_31 = (ncol_bot > max_tree42);
2581       if (ncol_bot <= max_tree42) {
2582         if ((ncol_bot/3)*3 == ncol_bot || (ncol_bot/2)*2 != ncol_bot) {
2583           is_tree_31 = true;
2584         }
2585         else {
2586           is_tree_42 = true;
2587         }
2588       }
2589     }
2590
2591     const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
2592     const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
2593     const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1);
2594     const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
2595
2596     if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl)
2597       return error(COMPERR_BAD_INPUT_MESH);
2598
2599     // arrays for normalized params
2600     TColStd_SequenceOfReal npb, npr, npt, npl;
2601     for (j = 0; j < nb; j++) {
2602       npb.Append(uv_eb[j].normParam);
2603     }
2604     for (i = 0; i < nr; i++) {
2605       npr.Append(uv_er[i].normParam);
2606     }
2607     for (j = 0; j < nt; j++) {
2608       npt.Append(uv_et[j].normParam);
2609     }
2610     for (i = 0; i < nl; i++) {
2611       npl.Append(uv_el[i].normParam);
2612     }
2613
2614     // We will ajust new points to this grid
2615     if (!SetNormalizedGrid(aMesh, aShape, quad))
2616       return false;
2617
2618     // TODO ???
2619     gp_XY a0 (uv_eb.front().u, uv_eb.front().v);
2620     gp_XY a1 (uv_eb.back().u,  uv_eb.back().v);
2621     gp_XY a2 (uv_et.back().u,  uv_et.back().v);
2622     gp_XY a3 (uv_et.front().u, uv_et.front().v);
2623     //=========================================================
2624
2625     TColStd_SequenceOfInteger curr_base, next_base;
2626     TColStd_SequenceOfReal curr_par_u, curr_par_v;
2627     TColStd_SequenceOfReal next_par_u, next_par_v;
2628     StdMeshers_Array2OfNode NodesBRD (1,nb, 1,nr);
2629     for (j = 1; j <= nb; j++) {
2630       NodesBRD.SetValue(j, 1, uv_eb[j - 1].node); // bottom
2631       curr_base.Append(j);
2632       next_base.Append(-1);
2633       curr_par_u.Append(uv_eb[j-1].u);
2634       curr_par_v.Append(uv_eb[j-1].v);
2635       next_par_u.Append(0.);
2636       next_par_v.Append(0.);
2637     }
2638     for (j = 1; j <= nt; j++) {
2639       NodesBRD.SetValue(j, nr, uv_et[j - 1].node); // top
2640     }
2641
2642     int curr_base_len = nb;
2643     int next_base_len = 0;
2644
2645     if (is_tree_42) {
2646       // "tree" simple reduce "42": 2->4->8->16->32->...
2647       //
2648       //  .-------------------------------.-------------------------------. nr
2649       //  |    \                          |                          /    |
2650       //  |         \     .---------------.---------------.     /         |
2651       //  |               |               |               |               |
2652       //  .---------------.---------------.---------------.---------------.
2653       //  | \             |             / | \             |             / |
2654       //  |     \ .-------.-------. /     |     \ .-------.-------. /     |
2655       //  |       |       |       |       |       |       |       |       |
2656       //  .-------.-------.-------.-------.-------.-------.-------.-------. i
2657       //  |\      |      /|\      |      /|\      |      /|\      |      /|
2658       //  |  \.---.---./  |  \.---.---./  |  \.---.---./  |  \.---.---./  |
2659       //  |   |   |   |   |   |   |   |   |   |   |   |   |   |   |   |   |
2660       //  .---.---.---.---.---.---.---.---.---.---.---.---.---.---.---.---.
2661       //  |\  |  /|\  |  /|\  |  /|\  |  /|\  |  /|\  |  /|\  |  /|\  |  /|
2662       //  | .-.-. | .-.-. | .-.-. | .-.-. | .-.-. | .-.-. | .-.-. | .-.-. |
2663       //  | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | |
2664       //  .-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. 1
2665       //  1                               j                               nb
2666
2667       for (i = 1; i < nr; i++) { // layer by layer
2668         // left
2669         NodesBRD.SetValue(1, i+1, uv_el[i].node);
2670         next_base.SetValue(++next_base_len, 1);
2671         // right
2672         NodesBRD.SetValue(nb, i+1, uv_er[i].node);
2673
2674         next_par_u.SetValue(next_base_len, uv_el[i].u);
2675         next_par_v.SetValue(next_base_len, uv_el[i].v);
2676
2677         // to stop reducing, if number of nodes reaches nt
2678         int delta = curr_base_len - nt;
2679
2680         //double du = uv_er[i].u - uv_el[i].u;
2681         //double dv = uv_er[i].v - uv_el[i].v;
2682
2683         // to calculate normalized parameter, we must know number of points in next layer
2684         int nb_four = (curr_base_len - 1) / 4;
2685         int nb_next = nb_four*2 + (curr_base_len - nb_four*4);
2686         if (nb_next < nt) nb_next = nt;
2687
2688         for (j = 1; j + 4 <= curr_base_len && delta > 0; j += 4, delta -= 2) {
2689           // add one "HH": nodes a,b,c,d,e and faces 1,2,3,4,5,6
2690           //
2691           //  .-----a-----b i + 1
2692           //  |\ 5  | 6  /|
2693           //  | \   |   / |
2694           //  |  c--d--e  |
2695           //  |1 |2 |3 |4 |
2696           //  |  |  |  |  |
2697           //  .--.--.--.--. i
2698           //
2699           //  j     j+2   j+4
2700
2701           double u,v;
2702
2703           // a (i + 1, j + 2)
2704           const SMDS_MeshNode* Na;
2705           next_base_len++;
2706           next_base.SetValue(next_base_len, curr_base.Value(j + 2));
2707           if (i + 1 == nr) { // top
2708             Na = uv_et[next_base_len - 1].node;
2709             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na);
2710             u = uv_et[next_base_len - 1].u;
2711             v = uv_et[next_base_len - 1].v;
2712           }
2713           else {
2714             //double norm_par = double(next_base_len - 1)/double(nb_next - 1);
2715             //u = uv_el[i].u + du * norm_par;
2716             //v = uv_el[i].v + dv * norm_par;
2717             {
2718               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
2719               int nearest_node_j = (int)rel;
2720               rel -= nearest_node_j;
2721               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
2722               double u1 = quad->uv_grid[ij].u;
2723               double v1 = quad->uv_grid[ij].v;
2724               double u2 = quad->uv_grid[ij + 1].u;
2725               double v2 = quad->uv_grid[ij + 1].v;
2726               double duj = (u2 - u1) * rel;
2727               double dvj = (v2 - v1) * rel;
2728               u = u1 + duj;
2729               v = v1 + dvj;
2730             }
2731             //u = uv_el[i].u + du*npb.Value(curr_base.Value(j + 2));
2732             //v = uv_el[i].v + dv*npb.Value(curr_base.Value(j + 2));
2733             gp_Pnt P = S->Value(u,v);
2734             SMDS_MeshNode* Na1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
2735             meshDS->SetNodeOnFace(Na1, geomFaceID, u, v);
2736             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na1);
2737             Na = Na1;
2738           }
2739           next_par_u.SetValue(next_base_len, u);
2740           next_par_v.SetValue(next_base_len, v);
2741
2742           // b (i + 1, j + 4)
2743           const SMDS_MeshNode* Nb;
2744           next_base_len++;
2745           next_base.SetValue(next_base_len, curr_base.Value(j + 4));
2746           if (i + 1 == nr) { // top
2747             Nb = uv_et[next_base_len - 1].node;
2748             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb);
2749             u = uv_et[next_base_len - 1].u;
2750             v = uv_et[next_base_len - 1].v;
2751           }
2752           else if (j + 4 == curr_base_len) { // right
2753             Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
2754             u = uv_er[i].u;
2755             v = uv_er[i].v;
2756           }
2757           else {
2758             //double norm_par = double(next_base_len - 1)/double(nb_next - 1);
2759             //u = uv_el[i].u + du * norm_par;
2760             //v = uv_el[i].v + dv * norm_par;
2761             {
2762               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
2763               int nearest_node_j = (int)rel;
2764               rel -= nearest_node_j;
2765               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
2766               double u1 = quad->uv_grid[ij].u;
2767               double v1 = quad->uv_grid[ij].v;
2768               double u2 = quad->uv_grid[ij + 1].u;
2769               double v2 = quad->uv_grid[ij + 1].v;
2770               double duj = (u2 - u1) * rel;
2771               double dvj = (v2 - v1) * rel;
2772               u = u1 + duj;
2773               v = v1 + dvj;
2774             }
2775             //u = uv_el[i].u + du*npb.Value(curr_base.Value(j + 4));
2776             //v = uv_el[i].v + dv*npb.Value(curr_base.Value(j + 4));
2777             gp_Pnt P = S->Value(u,v);
2778             SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
2779             meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v);
2780             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1);
2781             Nb = Nb1;
2782           }
2783           next_par_u.SetValue(next_base_len, u);
2784           next_par_v.SetValue(next_base_len, v);
2785
2786           // c
2787           u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 2)) / 2.0;
2788           v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 2)) / 2.0;
2789           gp_Pnt P = S->Value(u,v);
2790           SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z());
2791           meshDS->SetNodeOnFace(Nc, geomFaceID, u, v);
2792
2793           // d
2794           u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 1)) / 2.0;
2795           v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 1)) / 2.0;
2796           P = S->Value(u,v);
2797           SMDS_MeshNode* Nd = meshDS->AddNode(P.X(), P.Y(), P.Z());
2798           meshDS->SetNodeOnFace(Nd, geomFaceID, u, v);
2799
2800           // e
2801           u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len)) / 2.0;
2802           v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len)) / 2.0;
2803           P = S->Value(u,v);
2804           SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z());
2805           meshDS->SetNodeOnFace(Ne, geomFaceID, u, v);
2806
2807           // Faces
2808           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i),
2809                                               NodesBRD.Value(curr_base.Value(j + 1), i),
2810                                               Nc,
2811                                               NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1));
2812           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
2813
2814           SMDS_MeshFace* F2 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i),
2815                                               NodesBRD.Value(curr_base.Value(j + 2), i),
2816                                               Nd, Nc);
2817           if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID);
2818
2819           SMDS_MeshFace* F3 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i),
2820                                               NodesBRD.Value(curr_base.Value(j + 3), i),
2821                                               Ne, Nd);
2822           if (F3) meshDS->SetMeshElementOnShape(F3, geomFaceID);
2823
2824           SMDS_MeshFace* F4 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 3), i),
2825                                               NodesBRD.Value(curr_base.Value(j + 4), i),
2826                                               Nb, Ne);
2827           if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID);
2828
2829           SMDS_MeshFace* F5 = myHelper->AddFace(Nc, Nd, Na,
2830                                               NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1));
2831           if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID);
2832
2833           SMDS_MeshFace* F6 = myHelper->AddFace(Nd, Ne, Nb, Na);
2834           if (F6) meshDS->SetMeshElementOnShape(F6, geomFaceID);
2835         }
2836
2837         // not reduced side elements (if any)
2838         for (; j < curr_base_len; j++) {
2839           // f (i + 1, j + 1)
2840           const SMDS_MeshNode* Nf;
2841           double u,v;
2842           next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
2843           if (i + 1 == nr) { // top
2844             Nf = uv_et[next_base_len - 1].node;
2845             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
2846             u = uv_et[next_base_len - 1].u;
2847             v = uv_et[next_base_len - 1].v;
2848           }
2849           else if (j + 1 == curr_base_len) { // right
2850             Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
2851             u = uv_er[i].u;
2852             v = uv_er[i].v;
2853           }
2854           else {
2855             //double norm_par = double(next_base_len - 1)/double(nb_next - 1);
2856             //u = uv_el[i].u + du * norm_par;
2857             //v = uv_el[i].v + dv * norm_par;
2858             {
2859               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
2860               int nearest_node_j = (int)rel;
2861               rel -= nearest_node_j;
2862               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
2863               double u1 = quad->uv_grid[ij].u;
2864               double v1 = quad->uv_grid[ij].v;
2865               double u2 = quad->uv_grid[ij + 1].u;
2866               double v2 = quad->uv_grid[ij + 1].v;
2867               double duj = (u2 - u1) * rel;
2868               double dvj = (v2 - v1) * rel;
2869               u = u1 + duj;
2870               v = v1 + dvj;
2871             }
2872             //u = uv_el[i].u + du*npb.Value(curr_base.Value(j + 1));
2873             //v = uv_el[i].v + dv*npb.Value(curr_base.Value(j + 1));
2874             gp_Pnt P = S->Value(u,v);
2875             SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
2876             meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
2877             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
2878             Nf = Nf1;
2879           }
2880           next_par_u.SetValue(next_base_len, u);
2881           next_par_v.SetValue(next_base_len, v);
2882           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
2883                                               NodesBRD.Value(curr_base.Value(j + 1), i),
2884                                               NodesBRD.Value(next_base.Value(next_base_len), i + 1),
2885                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
2886           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
2887         }
2888
2889         curr_base_len = next_base_len;
2890         curr_base = next_base;
2891         curr_par_u = next_par_u;
2892         curr_par_v = next_par_v;
2893         next_base_len = 0;
2894       }
2895     } // end "tree" simple reduce "42"
2896     else if (is_tree_31) {
2897       // "tree" simple reduce "31": 1->3->9->27->...
2898       //
2899       //  .-----------------------------------------------------. nr
2900       //  |        \                                   /        |
2901       //  |                 .-----------------.                 |
2902       //  |                 |                 |                 |
2903       //  .-----------------.-----------------.-----------------.
2904       //  |   \         /   |   \         /   |   \         /   |
2905       //  |     .-----.     |     .-----.     |     .-----.     | i
2906       //  |     |     |     |     |     |     |     |     |     |
2907       //  .-----.-----.-----.-----.-----.-----.-----.-----.-----.
2908       //  |\   /|\   /|\   /|\   /|\   /|\   /|\   /|\   /|\   /|
2909       //  | .-. | .-. | .-. | .-. | .-. | .-. | .-. | .-. | .-. |
2910       //  | | | | | | | | | | | | | | | | | | | | | | | | | | | |
2911       //  .-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. 1
2912       //  1                          j                          nb
2913
2914       for (i = 1; i < nr; i++) { // layer by layer
2915         // left
2916         NodesBRD.SetValue(1, i+1, uv_el[i].node);
2917         next_base.SetValue(++next_base_len, 1);
2918         // right
2919         NodesBRD.SetValue(nb, i+1, uv_er[i].node);
2920
2921         next_par_u.SetValue(next_base_len, uv_el[i].u);
2922         next_par_v.SetValue(next_base_len, uv_el[i].v);
2923
2924         // to stop reducing, if number of nodes reaches nt
2925         int delta = curr_base_len - nt;
2926
2927         // to calculate normalized parameter, we must know number of points in next layer
2928         int nb_three = (curr_base_len - 1) / 3;
2929         int nb_next = nb_three + (curr_base_len - nb_three*3);
2930         if (nb_next < nt) nb_next = nt;
2931
2932         for (j = 1; j + 3 <= curr_base_len && delta > 0; j += 3, delta -= 2) {
2933           // add one "H": nodes b,c,e and faces 1,2,4,5
2934           //
2935           //  .---------b i + 1
2936           //  |\   5   /|
2937           //  | \     / |
2938           //  |  c---e  |
2939           //  |1 |2  |4 |
2940           //  |  |   |  |
2941           //  .--.---.--. i
2942           //
2943           //  j j+1 j+2 j+3
2944
2945           double u,v;
2946
2947           // b (i + 1, j + 3)
2948           const SMDS_MeshNode* Nb;
2949           next_base_len++;
2950           next_base.SetValue(next_base_len, curr_base.Value(j + 3));
2951           if (i + 1 == nr) { // top
2952             Nb = uv_et[next_base_len - 1].node;
2953             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb);
2954             u = uv_et[next_base_len - 1].u;
2955             v = uv_et[next_base_len - 1].v;
2956           }
2957           else if (j + 3 == curr_base_len) { // right
2958             Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
2959             u = uv_er[i].u;
2960             v = uv_er[i].v;
2961           }
2962           else {
2963             {
2964               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
2965               int nearest_node_j = (int)rel;
2966               rel -= nearest_node_j;
2967               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
2968               double u1 = quad->uv_grid[ij].u;
2969               double v1 = quad->uv_grid[ij].v;
2970               double u2 = quad->uv_grid[ij + 1].u;
2971               double v2 = quad->uv_grid[ij + 1].v;
2972               double duj = (u2 - u1) * rel;
2973               double dvj = (v2 - v1) * rel;
2974               u = u1 + duj;
2975               v = v1 + dvj;
2976             }
2977             gp_Pnt P = S->Value(u,v);
2978             SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
2979             meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v);
2980             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1);
2981             Nb = Nb1;
2982           }
2983           next_par_u.SetValue(next_base_len, u);
2984           next_par_v.SetValue(next_base_len, v);
2985
2986           // c and e
2987           double u1 = (curr_par_u.Value(j) + next_par_u.Value(next_base_len - 1)) / 2.0;
2988           double u2 = (curr_par_u.Value(j + 3) + next_par_u.Value(next_base_len)) / 2.0;
2989           double u3 = (u2 - u1) / 3.0;
2990
2991           double v1 = (curr_par_v.Value(j) + next_par_v.Value(next_base_len - 1)) / 2.0;
2992           double v2 = (curr_par_v.Value(j + 3) + next_par_v.Value(next_base_len)) / 2.0;
2993           double v3 = (v2 - v1) / 3.0;
2994
2995           // c
2996           u = u1 + u3;
2997           v = v1 + v3;
2998           gp_Pnt P = S->Value(u,v);
2999           SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z());
3000           meshDS->SetNodeOnFace(Nc, geomFaceID, u, v);
3001
3002           // e
3003           u = u1 + u3 + u3;
3004           v = v1 + v3 + v3;
3005           P = S->Value(u,v);
3006           SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z());
3007           meshDS->SetNodeOnFace(Ne, geomFaceID, u, v);
3008
3009           // Faces
3010           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i),
3011                                               NodesBRD.Value(curr_base.Value(j + 1), i),
3012                                               Nc,
3013                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3014           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3015
3016           SMDS_MeshFace* F2 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i),
3017                                               NodesBRD.Value(curr_base.Value(j + 2), i),
3018                                               Ne, Nc);
3019           if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID);
3020
3021           SMDS_MeshFace* F4 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i),
3022                                               NodesBRD.Value(curr_base.Value(j + 3), i),
3023                                               Nb, Ne);
3024           if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID);
3025
3026           SMDS_MeshFace* F5 = myHelper->AddFace(Nc, Ne, Nb,
3027                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3028           if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID);
3029         }
3030
3031         // not reduced side elements (if any)
3032         for (; j < curr_base_len; j++) {
3033           // f (i + 1, j + 1)
3034           const SMDS_MeshNode* Nf;
3035           double u,v;
3036           next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
3037           if (i + 1 == nr) { // top
3038             Nf = uv_et[next_base_len - 1].node;
3039             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
3040             u = uv_et[next_base_len - 1].u;
3041             v = uv_et[next_base_len - 1].v;
3042           }
3043           else if (j + 1 == curr_base_len) { // right
3044             Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
3045             u = uv_er[i].u;
3046             v = uv_er[i].v;
3047           }
3048           else {
3049             {
3050               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3051               int nearest_node_j = (int)rel;
3052               rel -= nearest_node_j;
3053               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3054               double u1 = quad->uv_grid[ij].u;
3055               double v1 = quad->uv_grid[ij].v;
3056               double u2 = quad->uv_grid[ij + 1].u;
3057               double v2 = quad->uv_grid[ij + 1].v;
3058               double duj = (u2 - u1) * rel;
3059               double dvj = (v2 - v1) * rel;
3060               u = u1 + duj;
3061               v = v1 + dvj;
3062             }
3063             gp_Pnt P = S->Value(u,v);
3064             SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3065             meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
3066             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
3067             Nf = Nf1;
3068           }
3069           next_par_u.SetValue(next_base_len, u);
3070           next_par_v.SetValue(next_base_len, v);
3071           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
3072                                               NodesBRD.Value(curr_base.Value(j + 1), i),
3073                                               NodesBRD.Value(next_base.Value(next_base_len), i + 1),
3074                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3075           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3076         }
3077
3078         curr_base_len = next_base_len;
3079         curr_base = next_base;
3080         curr_par_u = next_par_u;
3081         curr_par_v = next_par_v;
3082         next_base_len = 0;
3083       }
3084     } // end "tree" simple reduce "31"
3085     else if (is_lin_42) {
3086       // "linear" simple reduce "42": 4->8->12->16
3087       //
3088       //  .---------------.---------------.---------------.---------------. nr
3089       //  | \             |             / | \             |             / |
3090       //  |     \ .-------.-------. /     |     \ .-------.-------. /     |
3091       //  |       |       |       |       |       |       |       |       |
3092       //  .-------.-------.-------.-------.-------.-------.-------.-------.
3093       //  |      / \      |      / \      |      / \      |      / \      |
3094       //  |     /   \.----.----./   \     |     /   \.----.----./   \     | i
3095       //  |     /    |    |    |    \     |     /    |    |    |    \     |
3096       //  .-----.----.----.----.----.-----.-----.----.----.----.----.-----.
3097       //  |     /   / \   |  /  \   \     |     /   / \   |  /  \   \     |
3098       //  |    /   /    .-.-.    \   \    |    /   /    .-.-.    \   \    |
3099       //  |   /   /    /  |  \    \   \   |   /   /    /  |  \    \   \   |
3100       //  .---.---.---.---.---.---.---.---.---.---.---.---.---.---.---.---. 1
3101       //  1                               j                               nb
3102
3103       // nt = 5, nb = 7, nr = 4
3104       //int delta_all = 2;
3105       //int delta_one_col = 6;
3106       //int nb_col = 0;
3107       //int remainder = 2;
3108       //if (remainder > 0) nb_col++;
3109       //nb_col = 1;
3110       //int free_left = 1;
3111       //free_left += 2;
3112       //int free_middle = 4;
3113
3114       int delta_all = nb - nt;
3115       int delta_one_col = (nr - 1) * 2;
3116       int nb_col = delta_all / delta_one_col;
3117       int remainder = delta_all - nb_col * delta_one_col;
3118       if (remainder > 0) {
3119         nb_col++;
3120       }
3121       int free_left = ((nt - 1) - nb_col * 2) / 2;
3122       free_left += nr - 2;
3123       int free_middle = (nr - 2) * 2;
3124       if (remainder > 0 && nb_col == 1) {
3125         int nb_rows_short_col = remainder / 2;
3126         int nb_rows_thrown = (nr - 1) - nb_rows_short_col;
3127         free_left -= nb_rows_thrown;
3128       }
3129
3130       // nt = 5, nb = 17, nr = 4
3131       //int delta_all = 12;
3132       //int delta_one_col = 6;
3133       //int nb_col = 2;
3134       //int remainder = 0;
3135       //int free_left = 2;
3136       //int free_middle = 4;
3137
3138       for (i = 1; i < nr; i++, free_middle -= 2, free_left -= 1) { // layer by layer
3139         // left
3140         NodesBRD.SetValue(1, i+1, uv_el[i].node);
3141         next_base.SetValue(++next_base_len, 1);
3142         // right
3143         NodesBRD.SetValue(nb, i+1, uv_er[i].node);
3144
3145         // left
3146         next_par_u.SetValue(next_base_len, uv_el[i].u);
3147         next_par_v.SetValue(next_base_len, uv_el[i].v);
3148
3149         // to calculate normalized parameter, we must know number of points in next layer
3150         int nb_next = curr_base_len - nb_col * 2;
3151         if (remainder > 0 && i > remainder / 2)
3152           // take into account short "column"
3153           nb_next += 2;
3154         if (nb_next < nt) nb_next = nt;
3155
3156         // not reduced left elements
3157         for (j = 1; j <= free_left; j++) {
3158           // f (i + 1, j + 1)
3159           const SMDS_MeshNode* Nf;
3160           double u,v;
3161           next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
3162           if (i + 1 == nr) { // top
3163             Nf = uv_et[next_base_len - 1].node;
3164             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
3165             u = uv_et[next_base_len - 1].u;
3166             v = uv_et[next_base_len - 1].v;
3167           }
3168           else {
3169             {
3170               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3171               int nearest_node_j = (int)rel;
3172               rel -= nearest_node_j;
3173               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3174               double u1 = quad->uv_grid[ij].u;
3175               double v1 = quad->uv_grid[ij].v;
3176               double u2 = quad->uv_grid[ij + 1].u;
3177               double v2 = quad->uv_grid[ij + 1].v;
3178               double duj = (u2 - u1) * rel;
3179               double dvj = (v2 - v1) * rel;
3180               u = u1 + duj;
3181               v = v1 + dvj;
3182             }
3183             gp_Pnt P = S->Value(u,v);
3184             SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3185             meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
3186             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
3187             Nf = Nf1;
3188           }
3189           next_par_u.SetValue(next_base_len, u);
3190           next_par_v.SetValue(next_base_len, v);
3191           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
3192                                               NodesBRD.Value(curr_base.Value(j + 1), i),
3193                                               NodesBRD.Value(next_base.Value(next_base_len), i + 1),
3194                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3195           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3196         }
3197
3198         for (int icol = 1; icol <= nb_col; icol++) {
3199
3200           if (remainder > 0 && icol == nb_col && i > remainder / 2)
3201             // stop short "column"
3202             break;
3203
3204           // add one "HH": nodes a,b,c,d,e and faces 1,2,3,4,5,6
3205           //
3206           //  .-----a-----b i + 1
3207           //  |\ 5  | 6  /|
3208           //  | \   |   / |
3209           //  |  c--d--e  |
3210           //  |1 |2 |3 |4 |
3211           //  |  |  |  |  |
3212           //  .--.--.--.--. i
3213           //
3214           //  j     j+2   j+4
3215
3216           double u,v;
3217
3218           // a (i + 1, j + 2)
3219           const SMDS_MeshNode* Na;
3220           next_base_len++;
3221           next_base.SetValue(next_base_len, curr_base.Value(j + 2));
3222           if (i + 1 == nr) { // top
3223             Na = uv_et[next_base_len - 1].node;
3224             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na);
3225             u = uv_et[next_base_len - 1].u;
3226             v = uv_et[next_base_len - 1].v;
3227           }
3228           else {
3229             {
3230               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3231               int nearest_node_j = (int)rel;
3232               rel -= nearest_node_j;
3233               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3234               double u1 = quad->uv_grid[ij].u;
3235               double v1 = quad->uv_grid[ij].v;
3236               double u2 = quad->uv_grid[ij + 1].u;
3237               double v2 = quad->uv_grid[ij + 1].v;
3238               double duj = (u2 - u1) * rel;
3239               double dvj = (v2 - v1) * rel;
3240               u = u1 + duj;
3241               v = v1 + dvj;
3242             }
3243             gp_Pnt P = S->Value(u,v);
3244             SMDS_MeshNode* Na1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3245             meshDS->SetNodeOnFace(Na1, geomFaceID, u, v);
3246             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na1);
3247             Na = Na1;
3248           }
3249           next_par_u.SetValue(next_base_len, u);
3250           next_par_v.SetValue(next_base_len, v);
3251
3252           // b (i + 1, j + 4)
3253           const SMDS_MeshNode* Nb;
3254           next_base_len++;
3255           next_base.SetValue(next_base_len, curr_base.Value(j + 4));
3256           if (i + 1 == nr) { // top
3257             Nb = uv_et[next_base_len - 1].node;
3258             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb);
3259             u = uv_et[next_base_len - 1].u;
3260             v = uv_et[next_base_len - 1].v;
3261           }
3262           else if (j + 4 == curr_base_len) { // right
3263             Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
3264             u = uv_er[i].u;
3265             v = uv_er[i].v;
3266           }
3267           else {
3268             {
3269               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3270               int nearest_node_j = (int)rel;
3271               rel -= nearest_node_j;
3272               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3273               double u1 = quad->uv_grid[ij].u;
3274               double v1 = quad->uv_grid[ij].v;
3275               double u2 = quad->uv_grid[ij + 1].u;
3276               double v2 = quad->uv_grid[ij + 1].v;
3277               double duj = (u2 - u1) * rel;
3278               double dvj = (v2 - v1) * rel;
3279               u = u1 + duj;
3280               v = v1 + dvj;
3281             }
3282             gp_Pnt P = S->Value(u,v);
3283             SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3284             meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v);
3285             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1);
3286             Nb = Nb1;
3287           }
3288           next_par_u.SetValue(next_base_len, u);
3289           next_par_v.SetValue(next_base_len, v);
3290
3291           // c
3292           u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 2)) / 2.0;
3293           v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 2)) / 2.0;
3294           gp_Pnt P = S->Value(u,v);
3295           SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z());
3296           meshDS->SetNodeOnFace(Nc, geomFaceID, u, v);
3297
3298           // d
3299           u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 1)) / 2.0;
3300           v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 1)) / 2.0;
3301           P = S->Value(u,v);
3302           SMDS_MeshNode* Nd = meshDS->AddNode(P.X(), P.Y(), P.Z());
3303           meshDS->SetNodeOnFace(Nd, geomFaceID, u, v);
3304
3305           // e
3306           u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len)) / 2.0;
3307           v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len)) / 2.0;
3308           P = S->Value(u,v);
3309           SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z());
3310           meshDS->SetNodeOnFace(Ne, geomFaceID, u, v);
3311
3312           // Faces
3313           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i),
3314                                               NodesBRD.Value(curr_base.Value(j + 1), i),
3315                                               Nc,
3316                                               NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1));
3317           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3318
3319           SMDS_MeshFace* F2 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i),
3320                                               NodesBRD.Value(curr_base.Value(j + 2), i),
3321                                               Nd, Nc);
3322           if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID);
3323
3324           SMDS_MeshFace* F3 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i),
3325                                               NodesBRD.Value(curr_base.Value(j + 3), i),
3326                                               Ne, Nd);
3327           if (F3) meshDS->SetMeshElementOnShape(F3, geomFaceID);
3328
3329           SMDS_MeshFace* F4 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 3), i),
3330                                               NodesBRD.Value(curr_base.Value(j + 4), i),
3331                                               Nb, Ne);
3332           if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID);
3333
3334           SMDS_MeshFace* F5 = myHelper->AddFace(Nc, Nd, Na,
3335                                               NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1));
3336           if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID);
3337
3338           SMDS_MeshFace* F6 = myHelper->AddFace(Nd, Ne, Nb, Na);
3339           if (F6) meshDS->SetMeshElementOnShape(F6, geomFaceID);
3340
3341           j += 4;
3342
3343           // not reduced middle elements
3344           if (icol < nb_col) {
3345             if (remainder > 0 && icol == nb_col - 1 && i > remainder / 2)
3346               // pass middle elements before stopped short "column"
3347               break;
3348
3349             int free_add = free_middle;
3350             if (remainder > 0 && icol == nb_col - 1)
3351               // next "column" is short
3352               free_add -= (nr - 1) - (remainder / 2);
3353
3354             for (int imiddle = 1; imiddle <= free_add; imiddle++) {
3355               // f (i + 1, j + imiddle)
3356               const SMDS_MeshNode* Nf;
3357               double u,v;
3358               next_base.SetValue(++next_base_len, curr_base.Value(j + imiddle));
3359               if (i + 1 == nr) { // top
3360                 Nf = uv_et[next_base_len - 1].node;
3361                 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
3362                 u = uv_et[next_base_len - 1].u;
3363                 v = uv_et[next_base_len - 1].v;
3364               }
3365               else if (j + imiddle == curr_base_len) { // right
3366                 Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
3367                 u = uv_er[i].u;
3368                 v = uv_er[i].v;
3369               }
3370               else {
3371                 {
3372                   double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3373                   int nearest_node_j = (int)rel;
3374                   rel -= nearest_node_j;
3375                   int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3376                   double u1 = quad->uv_grid[ij].u;
3377                   double v1 = quad->uv_grid[ij].v;
3378                   double u2 = quad->uv_grid[ij + 1].u;
3379                   double v2 = quad->uv_grid[ij + 1].v;
3380                   double duj = (u2 - u1) * rel;
3381                   double dvj = (v2 - v1) * rel;
3382                   u = u1 + duj;
3383                   v = v1 + dvj;
3384                 }
3385                 gp_Pnt P = S->Value(u,v);
3386                 SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3387                 meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
3388                 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
3389                 Nf = Nf1;
3390               }
3391               next_par_u.SetValue(next_base_len, u);
3392               next_par_v.SetValue(next_base_len, v);
3393               SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j - 1 + imiddle), i),
3394                                                   NodesBRD.Value(curr_base.Value(j + imiddle), i),
3395                                                   NodesBRD.Value(next_base.Value(next_base_len), i + 1),
3396                                                   NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3397               if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3398             }
3399             j += free_add;
3400           }
3401         }
3402
3403         // not reduced right elements
3404         for (; j < curr_base_len; j++) {
3405           // f (i + 1, j + 1)
3406           const SMDS_MeshNode* Nf;
3407           double u,v;
3408           next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
3409           if (i + 1 == nr) { // top
3410             Nf = uv_et[next_base_len - 1].node;
3411             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
3412             u = uv_et[next_base_len - 1].u;
3413             v = uv_et[next_base_len - 1].v;
3414           }
3415           else if (j + 1 == curr_base_len) { // right
3416             Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
3417             u = uv_er[i].u;
3418             v = uv_er[i].v;
3419           }
3420           else {
3421             {
3422               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3423               int nearest_node_j = (int)rel;
3424               rel -= nearest_node_j;
3425               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3426               double u1 = quad->uv_grid[ij].u;
3427               double v1 = quad->uv_grid[ij].v;
3428               double u2 = quad->uv_grid[ij + 1].u;
3429               double v2 = quad->uv_grid[ij + 1].v;
3430               double duj = (u2 - u1) * rel;
3431               double dvj = (v2 - v1) * rel;
3432               u = u1 + duj;
3433               v = v1 + dvj;
3434             }
3435             gp_Pnt P = S->Value(u,v);
3436             SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3437             meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
3438             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
3439             Nf = Nf1;
3440           }
3441           next_par_u.SetValue(next_base_len, u);
3442           next_par_v.SetValue(next_base_len, v);
3443           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
3444                                               NodesBRD.Value(curr_base.Value(j + 1), i),
3445                                               NodesBRD.Value(next_base.Value(next_base_len), i + 1),
3446                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3447           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3448         }
3449
3450         curr_base_len = next_base_len;
3451         curr_base = next_base;
3452         curr_par_u = next_par_u;
3453         curr_par_v = next_par_v;
3454         next_base_len = 0;
3455       }
3456     } // end "linear" simple reduce "42"
3457     else if (is_lin_31) {
3458       // "linear" simple reduce "31": 2->6->10->14
3459       //
3460       //  .-----------------------------.-----------------------------. nr
3461       //  |     \                 /     |     \                 /     |
3462       //  |         .---------.         |         .---------.         |
3463       //  |         |         |         |         |         |         |
3464       //  .---------.---------.---------.---------.---------.---------.
3465       //  |        / \       / \        |        / \       / \        |
3466       //  |       /   .-----.   \       |       /   .-----.   \       | i
3467       //  |      /    |     |    \      |      /    |     |    \      |
3468       //  .-----.-----.-----.-----.-----.-----.-----.-----.-----.-----.
3469       //  |    /     / \   / \     \    |    /     / \   / \     \    |
3470       //  |   /     /   .-.   \     \   |   /     /   .-.   \     \   |
3471       //  |  /     /   /   \   \     \  |  /     /   /   \   \     \  |
3472       //  .--.----.---.-----.---.-----.-.--.----.---.-----.---.-----.-. 1
3473       //  1                             j                             nb
3474
3475       int delta_all = nb - nt;
3476       int delta_one_col = (nr - 1) * 2;
3477       int nb_col = delta_all / delta_one_col;
3478       int remainder = delta_all - nb_col * delta_one_col;
3479       if (remainder > 0) {
3480         nb_col++;
3481       }
3482       int free_left = ((nt - 1) - nb_col) / 2;
3483       free_left += nr - 2;
3484       int free_middle = (nr - 2) * 2;
3485       if (remainder > 0 && nb_col == 1) {
3486         int nb_rows_short_col = remainder / 2;
3487         int nb_rows_thrown = (nr - 1) - nb_rows_short_col;
3488         free_left -= nb_rows_thrown;
3489       }
3490
3491       for (i = 1; i < nr; i++, free_middle -= 2, free_left -= 1) { // layer by layer
3492         // left
3493         NodesBRD.SetValue(1, i+1, uv_el[i].node);
3494         next_base.SetValue(++next_base_len, 1);
3495         // right
3496         NodesBRD.SetValue(nb, i+1, uv_er[i].node);
3497
3498         // left
3499         next_par_u.SetValue(next_base_len, uv_el[i].u);
3500         next_par_v.SetValue(next_base_len, uv_el[i].v);
3501
3502         // to calculate normalized parameter, we must know number of points in next layer
3503         int nb_next = curr_base_len - nb_col * 2;
3504         if (remainder > 0 && i > remainder / 2)
3505           // take into account short "column"
3506           nb_next += 2;
3507         if (nb_next < nt) nb_next = nt;
3508
3509         // not reduced left elements
3510         for (j = 1; j <= free_left; j++) {
3511           // f (i + 1, j + 1)
3512           const SMDS_MeshNode* Nf;
3513           double u,v;
3514           next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
3515           if (i + 1 == nr) { // top
3516             Nf = uv_et[next_base_len - 1].node;
3517             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
3518             u = uv_et[next_base_len - 1].u;
3519             v = uv_et[next_base_len - 1].v;
3520           }
3521           else {
3522             {
3523               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3524               int nearest_node_j = (int)rel;
3525               rel -= nearest_node_j;
3526               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3527               double u1 = quad->uv_grid[ij].u;
3528               double v1 = quad->uv_grid[ij].v;
3529               double u2 = quad->uv_grid[ij + 1].u;
3530               double v2 = quad->uv_grid[ij + 1].v;
3531               double duj = (u2 - u1) * rel;
3532               double dvj = (v2 - v1) * rel;
3533               u = u1 + duj;
3534               v = v1 + dvj;
3535             }
3536             gp_Pnt P = S->Value(u,v);
3537             SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3538             meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
3539             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
3540             Nf = Nf1;
3541           }
3542           next_par_u.SetValue(next_base_len, u);
3543           next_par_v.SetValue(next_base_len, v);
3544           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
3545                                               NodesBRD.Value(curr_base.Value(j + 1), i),
3546                                               NodesBRD.Value(next_base.Value(next_base_len), i + 1),
3547                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3548           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3549         }
3550
3551         for (int icol = 1; icol <= nb_col; icol++) {
3552
3553           if (remainder > 0 && icol == nb_col && i > remainder / 2)
3554             // stop short "column"
3555             break;
3556
3557           // add one "H": nodes b,c,e and faces 1,2,4,5
3558           //
3559           //  .---------b i + 1
3560           //  |\   5   /|
3561           //  | \     / |
3562           //  |  c---e  |
3563           //  |1 |2  |4 |
3564           //  |  |   |  |
3565           //  .--.---.--. i
3566           //
3567           //  j j+1 j+2 j+3
3568
3569           double u,v;
3570
3571           // b (i + 1, j + 3)
3572           const SMDS_MeshNode* Nb;
3573           next_base_len++;
3574           next_base.SetValue(next_base_len, curr_base.Value(j + 3));
3575           if (i + 1 == nr) { // top
3576             Nb = uv_et[next_base_len - 1].node;
3577             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb);
3578             u = uv_et[next_base_len - 1].u;
3579             v = uv_et[next_base_len - 1].v;
3580           }
3581           else if (j + 3 == curr_base_len) { // right
3582             Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
3583             u = uv_er[i].u;
3584             v = uv_er[i].v;
3585           }
3586           else {
3587             {
3588               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3589               int nearest_node_j = (int)rel;
3590               rel -= nearest_node_j;
3591               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3592               double u1 = quad->uv_grid[ij].u;
3593               double v1 = quad->uv_grid[ij].v;
3594               double u2 = quad->uv_grid[ij + 1].u;
3595               double v2 = quad->uv_grid[ij + 1].v;
3596               double duj = (u2 - u1) * rel;
3597               double dvj = (v2 - v1) * rel;
3598               u = u1 + duj;
3599               v = v1 + dvj;
3600             }
3601             gp_Pnt P = S->Value(u,v);
3602             SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3603             meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v);
3604             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1);
3605             Nb = Nb1;
3606           }
3607           next_par_u.SetValue(next_base_len, u);
3608           next_par_v.SetValue(next_base_len, v);
3609
3610           // c and d
3611           double u1 = (curr_par_u.Value(j) + next_par_u.Value(next_base_len - 1)) / 2.0;
3612           double u2 = (curr_par_u.Value(j + 3) + next_par_u.Value(next_base_len)) / 2.0;
3613           double u3 = (u2 - u1) / 3.0;
3614
3615           double v1 = (curr_par_v.Value(j) + next_par_v.Value(next_base_len - 1)) / 2.0;
3616           double v2 = (curr_par_v.Value(j + 3) + next_par_v.Value(next_base_len)) / 2.0;
3617           double v3 = (v2 - v1) / 3.0;
3618
3619           // c
3620           u = u1 + u3;
3621           v = v1 + v3;
3622           gp_Pnt P = S->Value(u,v);
3623           SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z());
3624           meshDS->SetNodeOnFace(Nc, geomFaceID, u, v);
3625
3626           // e
3627           u = u1 + u3 + u3;
3628           v = v1 + v3 + v3;
3629           P = S->Value(u,v);
3630           SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z());
3631           meshDS->SetNodeOnFace(Ne, geomFaceID, u, v);
3632
3633           // Faces
3634           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i),
3635                                               NodesBRD.Value(curr_base.Value(j + 1), i),
3636                                               Nc,
3637                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3638           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3639
3640           SMDS_MeshFace* F2 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i),
3641                                               NodesBRD.Value(curr_base.Value(j + 2), i),
3642                                               Ne, Nc);
3643           if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID);
3644
3645           SMDS_MeshFace* F4 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i),
3646                                               NodesBRD.Value(curr_base.Value(j + 3), i),
3647                                               Nb, Ne);
3648           if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID);
3649
3650           SMDS_MeshFace* F5 = myHelper->AddFace(Nc, Ne, Nb,
3651                                               NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3652           if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID);
3653
3654           j += 3;
3655
3656           // not reduced middle elements
3657           if (icol < nb_col) {
3658             if (remainder > 0 && icol == nb_col - 1 && i > remainder / 2)
3659               // pass middle elements before stopped short "column"
3660               break;
3661
3662             int free_add = free_middle;
3663             if (remainder > 0 && icol == nb_col - 1)
3664               // next "column" is short
3665               free_add -= (nr - 1) - (remainder / 2);
3666
3667             for (int imiddle = 1; imiddle <= free_add; imiddle++) {
3668               // f (i + 1, j + imiddle)
3669               const SMDS_MeshNode* Nf;
3670               double u,v;
3671               next_base.SetValue(++next_base_len, curr_base.Value(j + imiddle));
3672               if (i + 1 == nr) { // top
3673                 Nf = uv_et[next_base_len - 1].node;
3674                 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
3675                 u = uv_et[next_base_len - 1].u;
3676                 v = uv_et[next_base_len - 1].v;
3677               }
3678               else if (j + imiddle == curr_base_len) { // right
3679                 Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
3680                 u = uv_er[i].u;
3681                 v = uv_er[i].v;
3682               }
3683               else {
3684                 {
3685                   double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3686                   int nearest_node_j = (int)rel;
3687                   rel -= nearest_node_j;
3688                   int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3689                   double u1 = quad->uv_grid[ij].u;
3690                   double v1 = quad->uv_grid[ij].v;
3691                   double u2 = quad->uv_grid[ij + 1].u;
3692                   double v2 = quad->uv_grid[ij + 1].v;
3693                   double duj = (u2 - u1) * rel;
3694                   double dvj = (v2 - v1) * rel;
3695                   u = u1 + duj;
3696                   v = v1 + dvj;
3697                 }
3698                 gp_Pnt P = S->Value(u,v);
3699                 SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3700                 meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
3701                 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
3702                 Nf = Nf1;
3703               }
3704               next_par_u.SetValue(next_base_len, u);
3705               next_par_v.SetValue(next_base_len, v);
3706               SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j - 1 + imiddle), i),
3707                                                   NodesBRD.Value(curr_base.Value(j + imiddle), i),
3708                                                   NodesBRD.Value(next_base.Value(next_base_len), i + 1),
3709                                                   NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3710               if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3711             }
3712             j += free_add;
3713           }
3714         }
3715
3716         // not reduced right elements
3717         for (; j < curr_base_len; j++) {
3718           // f (i + 1, j + 1)
3719           const SMDS_MeshNode* Nf;
3720           double u,v;
3721           next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
3722           if (i + 1 == nr) { // top
3723             Nf = uv_et[next_base_len - 1].node;
3724             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
3725             u = uv_et[next_base_len - 1].u;
3726             v = uv_et[next_base_len - 1].v;
3727           }
3728           else if (j + 1 == curr_base_len) { // right
3729             Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
3730             u = uv_er[i].u;
3731             v = uv_er[i].v;
3732           }
3733           else {
3734             {
3735               double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
3736               int nearest_node_j = (int)rel;
3737               rel -= nearest_node_j;
3738               int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
3739               double u1 = quad->uv_grid[ij].u;
3740               double v1 = quad->uv_grid[ij].v;
3741               double u2 = quad->uv_grid[ij + 1].u;
3742               double v2 = quad->uv_grid[ij + 1].v;
3743               double duj = (u2 - u1) * rel;
3744               double dvj = (v2 - v1) * rel;
3745               u = u1 + duj;
3746               v = v1 + dvj;
3747             }
3748             gp_Pnt P = S->Value(u,v);
3749             SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
3750             meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
3751             NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
3752             Nf = Nf1;
3753           }
3754           next_par_u.SetValue(next_base_len, u);
3755           next_par_v.SetValue(next_base_len, v);
3756           SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
3757                                                 NodesBRD.Value(curr_base.Value(j + 1), i),
3758                                                 NodesBRD.Value(next_base.Value(next_base_len), i + 1),
3759                                                 NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
3760           if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
3761         }
3762
3763         curr_base_len = next_base_len;
3764         curr_base = next_base;
3765         curr_par_u = next_par_u;
3766         curr_par_v = next_par_v;
3767         next_base_len = 0;
3768       }
3769     } // end "linear" simple reduce "31"
3770     else {
3771     }
3772   } // end Simple Reduce implementation
3773
3774   bool isOk = true;
3775   return isOk;
3776 }
3777
3778 //================================================================================
3779 namespace // data for smoothing
3780 {
3781   struct TSmoothNode;
3782   // --------------------------------------------------------------------------------
3783   /*!
3784    * \brief Structure used to check validity of node position after smoothing.
3785    *        It holds two nodes connected to a smoothed node and belonging to
3786    *        one mesh face
3787    */
3788   struct TTriangle
3789   {
3790     TSmoothNode* _n1;
3791     TSmoothNode* _n2;
3792     TTriangle( TSmoothNode* n1=0, TSmoothNode* n2=0 ): _n1(n1), _n2(n2) {}
3793
3794     inline bool IsForward( gp_UV uv );
3795   };
3796   // --------------------------------------------------------------------------------
3797   /*!
3798    * \brief Data of a smoothed node
3799    */
3800   struct TSmoothNode
3801   {
3802     gp_XY _uv;
3803     vector< TTriangle > _triangles; // if empty, then node is not movable
3804   };
3805   // --------------------------------------------------------------------------------
3806   inline bool TTriangle::IsForward( gp_UV uv )
3807   {
3808     gp_Vec2d v1( uv, _n1->_uv ), v2( uv, _n2->_uv );
3809     double d = v1 ^ v2;
3810     return d > 1e-100;
3811   }
3812 }
3813
3814 //================================================================================
3815 /*!
3816  * \brief Set UV of nodes on degenerated VERTEXes in the middle of degenerated EDGE
3817  *
3818  * WARNING: this method must be called AFTER retrieving UVPtStruct's from quad
3819  */
3820 //================================================================================
3821
3822 void StdMeshers_Quadrangle_2D::UpdateDegenUV(FaceQuadStruct* quad)
3823 {
3824   for ( unsigned i = 0; i < quad->side.size(); ++i )
3825   {
3826     StdMeshers_FaceSide* side = quad->side[i];
3827     const vector<UVPtStruct>& uvVec = side->GetUVPtStruct();
3828
3829     // find which end of the side is on degenerated shape
3830     int degenInd = -1;
3831     if ( myHelper->IsDegenShape( uvVec[0].node->getshapeId() ))
3832       degenInd = 0;
3833     else if ( myHelper->IsDegenShape( uvVec.back().node->getshapeId() ))
3834       degenInd = uvVec.size() - 1;
3835     else
3836       continue;
3837
3838     // find another side sharing the degenerated shape
3839     bool isPrev = ( degenInd == 0 );
3840     if ( i >= TOP_SIDE )
3841       isPrev = !isPrev;
3842     int i2 = ( isPrev ? ( i + 3 ) : ( i + 1 )) % 4;
3843     StdMeshers_FaceSide* side2 = quad->side[ i2 ];
3844     const vector<UVPtStruct>& uvVec2 = side2->GetUVPtStruct();
3845     int degenInd2 = -1;
3846     if ( uvVec[ degenInd ].node == uvVec2[0].node )
3847       degenInd2 = 0;
3848     else if ( uvVec[ degenInd ].node == uvVec2.back().node )
3849       degenInd2 = uvVec2.size() - 1;
3850     else
3851       throw SALOME_Exception( LOCALIZED( "Logical error" ));
3852
3853     // move UV in the middle
3854     uvPtStruct& uv1 = const_cast<uvPtStruct&>( uvVec [ degenInd  ]);
3855     uvPtStruct& uv2 = const_cast<uvPtStruct&>( uvVec2[ degenInd2 ]);
3856     uv1.u = uv2.u = 0.5 * ( uv1.u + uv2.u );
3857     uv1.v = uv2.v = 0.5 * ( uv1.v + uv2.v );
3858   }
3859 }
3860
3861 //================================================================================
3862 /*!
3863  * \brief Perform smoothing of 2D elements on a FACE with ignored degenerated EDGE
3864  */
3865 //================================================================================
3866
3867 void StdMeshers_Quadrangle_2D::Smooth (FaceQuadStruct* quad)
3868 {
3869   if ( !myNeedSmooth ) return;
3870
3871   // Get nodes to smooth
3872
3873   typedef map< const SMDS_MeshNode*, TSmoothNode, TIDCompare > TNo2SmooNoMap;
3874   TNo2SmooNoMap smooNoMap;
3875
3876   const TopoDS_Face& geomFace = TopoDS::Face( myHelper->GetSubShape() );
3877   SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
3878   SMESHDS_SubMesh* fSubMesh = meshDS->MeshElements( geomFace );
3879   SMDS_NodeIteratorPtr nIt = fSubMesh->GetNodes();
3880   while ( nIt->more() ) // loop on nodes bound to a FACE
3881   {
3882     const SMDS_MeshNode* node = nIt->next();
3883     TSmoothNode & sNode = smooNoMap[ node ];
3884     sNode._uv = myHelper->GetNodeUV( geomFace, node );
3885
3886     // set sNode._triangles
3887     SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator( SMDSAbs_Face );
3888     while ( fIt->more() )
3889     {
3890       const SMDS_MeshElement* face = fIt->next();
3891       const int nbN = face->NbCornerNodes();
3892       const int nInd = face->GetNodeIndex( node );
3893       const int prevInd = myHelper->WrapIndex( nInd - 1, nbN );
3894       const int nextInd = myHelper->WrapIndex( nInd + 1, nbN );
3895       const SMDS_MeshNode* prevNode = face->GetNode( prevInd );
3896       const SMDS_MeshNode* nextNode = face->GetNode( nextInd );
3897       sNode._triangles.push_back( TTriangle( & smooNoMap[ prevNode ],
3898                                              & smooNoMap[ nextNode ]));
3899     }
3900   }
3901   // set _uv of smooth nodes on FACE boundary
3902   for ( unsigned i = 0; i < quad->side.size(); ++i )
3903   {
3904     const vector<UVPtStruct>& uvVec = quad->side[i]->GetUVPtStruct();
3905     for ( unsigned j = 0; j < uvVec.size(); ++j )
3906     {
3907       TSmoothNode & sNode = smooNoMap[ uvVec[j].node ];
3908       sNode._uv.SetCoord( uvVec[j].u, uvVec[j].v );
3909     }
3910   }
3911
3912   // Smoothing
3913
3914   for ( int iLoop = 0; iLoop < 5; ++iLoop )
3915   {
3916     TNo2SmooNoMap::iterator n2sn = smooNoMap.begin();
3917     for ( ; n2sn != smooNoMap.end(); ++n2sn )
3918     {
3919       TSmoothNode& sNode = n2sn->second;
3920       if ( sNode._triangles.empty() )
3921         continue; // not movable node
3922
3923       // compute a new UV
3924       gp_XY newUV (0,0);
3925       for ( unsigned i = 0; i < sNode._triangles.size(); ++i )
3926         newUV += sNode._triangles[i]._n1->_uv;
3927       newUV /= sNode._triangles.size();
3928
3929       // check validity of the newUV
3930       bool isValid = true;
3931       for ( unsigned i = 0; i < sNode._triangles.size() && isValid; ++i )
3932         isValid = sNode._triangles[i].IsForward( newUV );
3933
3934       if ( isValid )
3935         sNode._uv = newUV;
3936     }
3937   }
3938
3939   // Set new XYZ to the smoothed nodes
3940
3941   Handle(Geom_Surface) surface = BRep_Tool::Surface( geomFace );
3942
3943   TNo2SmooNoMap::iterator n2sn = smooNoMap.begin();
3944   for ( ; n2sn != smooNoMap.end(); ++n2sn )
3945   {
3946     TSmoothNode& sNode = n2sn->second;
3947     if ( sNode._triangles.empty() )
3948       continue; // not movable node
3949
3950     SMDS_MeshNode* node = const_cast< SMDS_MeshNode*>( n2sn->first );
3951     gp_Pnt xyz = surface->Value( sNode._uv.X(), sNode._uv.Y() );
3952     meshDS->MoveNode( node, xyz.X(), xyz.Y(), xyz.Z() );
3953
3954     // store the new UV
3955     node->SetPosition( SMDS_PositionPtr( new SMDS_FacePosition( sNode._uv.X(), sNode._uv.Y() )));
3956   }
3957
3958   // Move medium nodes in quadratic mesh
3959   if ( _quadraticMesh )
3960   {
3961     const TLinkNodeMap& links = myHelper->GetTLinkNodeMap();
3962     TLinkNodeMap::const_iterator linkIt = links.begin();
3963     for ( ; linkIt != links.end(); ++linkIt )
3964     {
3965       const SMESH_TLink& link = linkIt->first;
3966       SMDS_MeshNode*     node = const_cast< SMDS_MeshNode*>( linkIt->second );
3967
3968       if ( node->getshapeId() != myHelper->GetSubShapeID() )
3969         continue; // medium node is on EDGE or VERTEX
3970
3971       gp_XY uv1 = myHelper->GetNodeUV( geomFace, link.node1(), node );
3972       gp_XY uv2 = myHelper->GetNodeUV( geomFace, link.node2(), node );
3973
3974       gp_XY uv  = myHelper->GetMiddleUV( surface, uv1, uv2 );
3975       node->SetPosition( SMDS_PositionPtr( new SMDS_FacePosition( uv.X(), uv.Y() )));
3976
3977       gp_Pnt xyz = surface->Value( uv.X(), uv.Y() );
3978       meshDS->MoveNode( node, xyz.X(), xyz.Y(), xyz.Z() );
3979     }
3980   }
3981 }