src/StdMeshers/StdMeshers_Quadrangle_2D.cxx

   1 //  Copyright (C) 2007-2008  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 //  SMESH SMESH : implementaion of SMESH idl descriptions
  23 //  File   : StdMeshers_Quadrangle_2D.cxx
  24 //           Moved here from SMESH_Quadrangle_2D.cxx
  25 //  Author : Paul RASCLE, EDF
  26 //  Module : SMESH
  27 //
  28 #include "StdMeshers_Quadrangle_2D.hxx"
  29
  30 #include "StdMeshers_FaceSide.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 <BRepTools_WireExplorer.hxx>
  45 #include <BRep_Tool.hxx>
  46 #include <Geom_Surface.hxx>
  47 #include <NCollection_DefineArray2.hxx>
  48 #include <Precision.hxx>
  49 #include <TColStd_SequenceOfReal.hxx>
  50 #include <TColgp_SequenceOfXY.hxx>
  51 #include <TopExp.hxx>
  52 #include <TopTools_ListIteratorOfListOfShape.hxx>
  53 #include <TopoDS.hxx>
  54
  55 #include "utilities.h"
  56 #include "Utils_ExceptHandlers.hxx"
  57
  58 #ifndef StdMeshers_Array2OfNode_HeaderFile
  59 #define StdMeshers_Array2OfNode_HeaderFile
  60 typedef const SMDS_MeshNode* SMDS_MeshNodePtr;
  61 DEFINE_BASECOLLECTION (StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr)
  62 DEFINE_ARRAY2(StdMeshers_Array2OfNode,
  63               StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr)
  64 #endif
  65
  66 using namespace std;
  67
  68 typedef gp_XY gp_UV;
  69 typedef SMESH_Comment TComm;
  70
  71 //=============================================================================
  72 /*!
  73  *
  74  */
  75 //=============================================================================
  76
  77 StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, SMESH_Gen* gen)
  78      : SMESH_2D_Algo(hypId, studyId, gen)
  79 {
  80   MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D");
  81   _name = "Quadrangle_2D";
  82   _shapeType = (1 << TopAbs_FACE);
  83   _compatibleHypothesis.push_back("QuadranglePreference");
  84   _compatibleHypothesis.push_back("TrianglePreference");
  85   myTool = 0;
  86 }
  87
  88 //=============================================================================
  89 /*!
  90  *
  91  */
  92 //=============================================================================
  93
  94 StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D()
  95 {
  96   MESSAGE("StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D");
  97 }
  98
  99 //=============================================================================
 100 /*!
 101  *
 102  */
 103 //=============================================================================
 104
 105 bool StdMeshers_Quadrangle_2D::CheckHypothesis
 106                          (SMESH_Mesh&                          aMesh,
 107                           const TopoDS_Shape&                  aShape,
 108                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 109 {
 110   bool isOk = true;
 111   aStatus = SMESH_Hypothesis::HYP_OK;
 112
 113
 114   const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape, false);
 115   const SMESHDS_Hypothesis *theHyp = 0;
 116
 117   if(hyps.size() > 0){
 118     theHyp = *hyps.begin();
 119     if(strcmp("QuadranglePreference", theHyp->GetName()) == 0) {
 120       myQuadranglePreference= true;
 121       myTrianglePreference= false;
 122     }
 123     else if(strcmp("TrianglePreference", theHyp->GetName()) == 0){
 124       myQuadranglePreference= false;
 125       myTrianglePreference= true;
 126     }
 127   }
 128   else {
 129     myQuadranglePreference = false;
 130     myTrianglePreference = false;
 131   }
 132   return isOk;
 133 }
 134
 135 //=============================================================================
 136 /*!
 137  *
 138  */
 139 //=============================================================================
 140
 141 bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh,
 142                                         const TopoDS_Shape& aShape)// throw (SALOME_Exception)
 143 {
 144   // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside
 145   //Unexpect aCatchSalomeException);
 146
 147   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 148   aMesh.GetSubMesh(aShape);
 149
 150   SMESH_MesherHelper helper(aMesh);
 151   myTool = &helper;
 152
 153   _quadraticMesh = myTool->IsQuadraticSubMesh(aShape);
 154
 155   FaceQuadStruct *quad = CheckNbEdges( aMesh, aShape );
 156   std::auto_ptr<FaceQuadStruct> quadDeleter( quad ); // to delete quad at exit from Compute()
 157   if (!quad)
 158     return false;
 159
 160   if(myQuadranglePreference) {
 161     int n1 = quad->side[0]->NbPoints();
 162     int n2 = quad->side[1]->NbPoints();
 163     int n3 = quad->side[2]->NbPoints();
 164     int n4 = quad->side[3]->NbPoints();
 165     int nfull = n1+n2+n3+n4;
 166     int ntmp = nfull/2;
 167     ntmp = ntmp*2;
 168     if( nfull==ntmp && ( (n1!=n3) || (n2!=n4) ) ) {
 169       // special path for using only quandrangle faces
 170       bool ok = ComputeQuadPref(aMesh, aShape, quad);
 171       return ok;
 172     }
 173   }
 174
 175   // set normalized grid on unit square in parametric domain
 176
 177   if (!SetNormalizedGrid(aMesh, aShape, quad))
 178     return false;
 179
 180   // --- compute 3D values on points, store points & quadrangles
 181
 182   int nbdown  = quad->side[0]->NbPoints();
 183   int nbup    = quad->side[2]->NbPoints();
 184
 185   int nbright = quad->side[1]->NbPoints();
 186   int nbleft  = quad->side[3]->NbPoints();
 187
 188   int nbhoriz  = Min(nbdown, nbup);
 189   int nbvertic = Min(nbright, nbleft);
 190
 191   const TopoDS_Face& F = TopoDS::Face(aShape);
 192   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
 193
 194   // internal mesh nodes
 195   int i, j, geomFaceID = meshDS->ShapeToIndex( F );
 196   for (i = 1; i < nbhoriz - 1; i++) {
 197     for (j = 1; j < nbvertic - 1; j++) {
 198       int ij = j * nbhoriz + i;
 199       double u = quad->uv_grid[ij].u;
 200       double v = quad->uv_grid[ij].v;
 201       gp_Pnt P = S->Value(u, v);
 202       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 203       meshDS->SetNodeOnFace(node, geomFaceID, u, v);
 204       quad->uv_grid[ij].node = node;
 205     }
 206   }
 207
 208   // mesh faces
 209
 210   //             [2]
 211   //      --.--.--.--.--.--  nbvertic
 212   //     |                 | ^
 213   //     |                 | ^
 214   // [3] |                 | ^ j  [1]
 215   //     |                 | ^
 216   //     |                 | ^
 217   //      ---.----.----.---  0
 218   //     0 > > > > > > > > nbhoriz
 219   //              i
 220   //             [0]
 221
 222   i = 0;
 223   int ilow = 0;
 224   int iup = nbhoriz - 1;
 225   if (quad->isEdgeOut[3]) { ilow++; } else { if (quad->isEdgeOut[1]) iup--; }
 226
 227   int jlow = 0;
 228   int jup = nbvertic - 1;
 229   if (quad->isEdgeOut[0]) { jlow++; } else { if (quad->isEdgeOut[2]) jup--; }
 230
 231   // regular quadrangles
 232   for (i = ilow; i < iup; i++) {
 233     for (j = jlow; j < jup; j++) {
 234       const SMDS_MeshNode *a, *b, *c, *d;
 235       a = quad->uv_grid[j * nbhoriz + i].node;
 236       b = quad->uv_grid[j * nbhoriz + i + 1].node;
 237       c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node;
 238       d = quad->uv_grid[(j + 1) * nbhoriz + i].node;
 239       SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
 240       meshDS->SetMeshElementOnShape(face, geomFaceID);
 241     }
 242   }
 243
 244   const vector<UVPtStruct>& uv_e0 = quad->side[0]->GetUVPtStruct(true,0 );
 245   const vector<UVPtStruct>& uv_e1 = quad->side[1]->GetUVPtStruct(false,1);
 246   const vector<UVPtStruct>& uv_e2 = quad->side[2]->GetUVPtStruct(true,1 );
 247   const vector<UVPtStruct>& uv_e3 = quad->side[3]->GetUVPtStruct(false,0);
 248
 249   if ( uv_e0.empty() || uv_e1.empty() || uv_e2.empty() || uv_e3.empty() )
 250     return error( COMPERR_BAD_INPUT_MESH );
 251
 252   double eps = Precision::Confusion();
 253
 254   // Boundary quadrangles
 255
 256   if (quad->isEdgeOut[0]) {
 257     // Down edge is out
 258     //
 259     // |___|___|___|___|___|___|
 260     // |   |   |   |   |   |   |
 261     // |___|___|___|___|___|___|
 262     // |   |   |   |   |   |   |
 263     // |___|___|___|___|___|___| __ first row of the regular grid
 264     // .  .  .  .  .  .  .  .  . __ down edge nodes
 265     //
 266     // >->->->->->->->->->->->-> -- direction of processing
 267
 268     int g = 0; // number of last processed node in the regular grid
 269
 270     // number of last node of the down edge to be processed
 271     int stop = nbdown - 1;
 272     // if right edge is out, we will stop at a node, previous to the last one
 273     if (quad->isEdgeOut[1]) stop--;
 274
 275     // for each node of the down edge find nearest node
 276     // in the first row of the regular grid and link them
 277     for (i = 0; i < stop; i++) {
 278       const SMDS_MeshNode *a, *b, *c, *d;
 279       a = uv_e0[i].node;
 280       b = uv_e0[i + 1].node;
 281       gp_Pnt pb (b->X(), b->Y(), b->Z());
 282
 283       // find node c in the regular grid, which will be linked with node b
 284       int near = g;
 285       if (i == stop - 1) {
 286         // right bound reached, link with the rightmost node
 287         near = iup;
 288         c = quad->uv_grid[nbhoriz + iup].node;
 289       }
 290       else {
 291         // find in the grid node c, nearest to the b
 292         double mind = RealLast();
 293         for (int k = g; k <= iup; k++) {
 294
 295           const SMDS_MeshNode *nk;
 296           if (k < ilow) // this can be, if left edge is out
 297             nk = uv_e3[1].node; // get node from the left edge
 298           else
 299             nk = quad->uv_grid[nbhoriz + k].node; // get one of middle nodes
 300
 301           gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 302           double dist = pb.Distance(pnk);
 303           if (dist < mind - eps) {
 304             c = nk;
 305             near = k;
 306             mind = dist;
 307           } else {
 308             break;
 309           }
 310         }
 311       }
 312
 313       if (near == g) { // make triangle
 314         //SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
 315         SMDS_MeshFace* face = myTool->AddFace(a, b, c);
 316         meshDS->SetMeshElementOnShape(face, geomFaceID);
 317       }
 318       else { // make quadrangle
 319         if (near - 1 < ilow)
 320           d = uv_e3[1].node;
 321         else
 322           d = quad->uv_grid[nbhoriz + near - 1].node;
 323         //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 324
 325         if(!myTrianglePreference){
 326           SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
 327           meshDS->SetMeshElementOnShape(face, geomFaceID);
 328         }
 329         else {
 330           SplitQuad(meshDS, geomFaceID, a, b, c, d);
 331         }
 332
 333         // if node d is not at position g - make additional triangles
 334         if (near - 1 > g) {
 335           for (int k = near - 1; k > g; k--) {
 336             c = quad->uv_grid[nbhoriz + k].node;
 337             if (k - 1 < ilow)
 338               d = uv_e3[1].node;
 339             else
 340               d = quad->uv_grid[nbhoriz + k - 1].node;
 341             //SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
 342             SMDS_MeshFace* face = myTool->AddFace(a, c, d);
 343             meshDS->SetMeshElementOnShape(face, geomFaceID);
 344           }
 345         }
 346         g = near;
 347       }
 348     }
 349   } else {
 350     if (quad->isEdgeOut[2]) {
 351       // Up edge is out
 352       //
 353       // <-<-<-<-<-<-<-<-<-<-<-<-< -- direction of processing
 354       //
 355       // .  .  .  .  .  .  .  .  . __ up edge nodes
 356       //  ___ ___ ___ ___ ___ ___  __ first row of the regular grid
 357       // |   |   |   |   |   |   |
 358       // |___|___|___|___|___|___|
 359       // |   |   |   |   |   |   |
 360       // |___|___|___|___|___|___|
 361       // |   |   |   |   |   |   |
 362
 363       int g = nbhoriz - 1; // last processed node in the regular grid
 364
 365       int stop = 0;
 366       // if left edge is out, we will stop at a second node
 367       if (quad->isEdgeOut[3]) stop++;
 368
 369       // for each node of the up edge find nearest node
 370       // in the first row of the regular grid and link them
 371       for (i = nbup - 1; i > stop; i--) {
 372         const SMDS_MeshNode *a, *b, *c, *d;
 373         a = uv_e2[i].node;
 374         b = uv_e2[i - 1].node;
 375         gp_Pnt pb (b->X(), b->Y(), b->Z());
 376
 377         // find node c in the grid, which will be linked with node b
 378         int near = g;
 379         if (i == stop + 1) { // left bound reached, link with the leftmost node
 380           c = quad->uv_grid[nbhoriz*(nbvertic - 2) + ilow].node;
 381           near = ilow;
 382         } else {
 383           // find node c in the grid, nearest to the b
 384           double mind = RealLast();
 385           for (int k = g; k >= ilow; k--) {
 386             const SMDS_MeshNode *nk;
 387             if (k > iup)
 388               nk = uv_e1[nbright - 2].node;
 389             else
 390               nk = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
 391             gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 392             double dist = pb.Distance(pnk);
 393             if (dist < mind - eps) {
 394               c = nk;
 395               near = k;
 396               mind = dist;
 397             } else {
 398               break;
 399             }
 400           }
 401         }
 402
 403         if (near == g) { // make triangle
 404           //SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
 405           SMDS_MeshFace* face = myTool->AddFace(a, b, c);
 406           meshDS->SetMeshElementOnShape(face, geomFaceID);
 407         }
 408         else { // make quadrangle
 409           if (near + 1 > iup)
 410             d = uv_e1[nbright - 2].node;
 411           else
 412             d = quad->uv_grid[nbhoriz*(nbvertic - 2) + near + 1].node;
 413           //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 414           if(!myTrianglePreference){
 415             SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
 416             meshDS->SetMeshElementOnShape(face, geomFaceID);
 417           }
 418           else {
 419             SplitQuad(meshDS, geomFaceID, a, b, c, d);
 420           }
 421
 422           if (near + 1 < g) { // if d not is at g - make additional triangles
 423             for (int k = near + 1; k < g; k++) {
 424               c = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
 425               if (k + 1 > iup)
 426                 d = uv_e1[nbright - 2].node;
 427               else
 428                 d = quad->uv_grid[nbhoriz*(nbvertic - 2) + k + 1].node;
 429               //SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
 430               SMDS_MeshFace* face = myTool->AddFace(a, c, d);
 431               meshDS->SetMeshElementOnShape(face, geomFaceID);
 432             }
 433           }
 434           g = near;
 435         }
 436       }
 437     }
 438   }
 439
 440   // right or left boundary quadrangles
 441   if (quad->isEdgeOut[1]) {
 442 //    MESSAGE("right edge is out");
 443     int g = 0; // last processed node in the grid
 444     int stop = nbright - 1;
 445     if (quad->isEdgeOut[2]) stop--;
 446     for (i = 0; i < stop; i++) {
 447       const SMDS_MeshNode *a, *b, *c, *d;
 448       a = uv_e1[i].node;
 449       b = uv_e1[i + 1].node;
 450       gp_Pnt pb (b->X(), b->Y(), b->Z());
 451
 452       // find node c in the grid, nearest to the b
 453       int near = g;
 454       if (i == stop - 1) { // up bondary reached
 455         c = quad->uv_grid[nbhoriz*(jup + 1) - 2].node;
 456         near = jup;
 457       } else {
 458         double mind = RealLast();
 459         for (int k = g; k <= jup; k++) {
 460           const SMDS_MeshNode *nk;
 461           if (k < jlow)
 462             nk = uv_e0[nbdown - 2].node;
 463           else
 464             nk = quad->uv_grid[nbhoriz*(k + 1) - 2].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 = meshDS->AddFace(a, b, c);
 479         SMDS_MeshFace* face = myTool->AddFace(a, b, c);
 480         meshDS->SetMeshElementOnShape(face, geomFaceID);
 481       }
 482       else { // make quadrangle
 483         if (near - 1 < jlow)
 484           d = uv_e0[nbdown - 2].node;
 485         else
 486           d = quad->uv_grid[nbhoriz*near - 2].node;
 487         //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 488
 489         if(!myTrianglePreference){
 490           SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
 491           meshDS->SetMeshElementOnShape(face, geomFaceID);
 492         }
 493         else {
 494           SplitQuad(meshDS, geomFaceID, a, b, c, d);
 495         }
 496
 497         if (near - 1 > g) { // if d not is at g - make additional triangles
 498           for (int k = near - 1; k > g; k--) {
 499             c = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
 500             if (k - 1 < jlow)
 501               d = uv_e0[nbdown - 2].node;
 502             else
 503               d = quad->uv_grid[nbhoriz*k - 2].node;
 504             //SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
 505             SMDS_MeshFace* face = myTool->AddFace(a, c, d);
 506             meshDS->SetMeshElementOnShape(face, geomFaceID);
 507           }
 508         }
 509         g = near;
 510       }
 511     }
 512   } else {
 513     if (quad->isEdgeOut[3]) {
 514 //      MESSAGE("left edge is out");
 515       int g = nbvertic - 1; // last processed node in the grid
 516       int stop = 0;
 517       if (quad->isEdgeOut[0]) stop++;
 518       for (i = nbleft - 1; i > stop; i--) {
 519         const SMDS_MeshNode *a, *b, *c, *d;
 520         a = uv_e3[i].node;
 521         b = uv_e3[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) { // down bondary reached
 527           c = quad->uv_grid[nbhoriz*jlow + 1].node;
 528           near = jlow;
 529         } else {
 530           double mind = RealLast();
 531           for (int k = g; k >= jlow; k--) {
 532             const SMDS_MeshNode *nk;
 533             if (k > jup)
 534               nk = uv_e2[1].node;
 535             else
 536               nk = quad->uv_grid[nbhoriz*k + 1].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 = meshDS->AddFace(a, b, c);
 551           SMDS_MeshFace* face = myTool->AddFace(a, b, c);
 552           meshDS->SetMeshElementOnShape(face, geomFaceID);
 553         }
 554         else { // make quadrangle
 555           if (near + 1 > jup)
 556             d = uv_e2[1].node;
 557           else
 558             d = quad->uv_grid[nbhoriz*(near + 1) + 1].node;
 559           //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 560           if(!myTrianglePreference){
 561             SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
 562             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].node;
 571               if (k + 1 > jup)
 572                 d = uv_e2[1].node;
 573               else
 574                 d = quad->uv_grid[nbhoriz*(k + 1) + 1].node;
 575               //SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
 576               SMDS_MeshFace* face = myTool->AddFace(a, c, d);
 577               meshDS->SetMeshElementOnShape(face, geomFaceID);
 578             }
 579           }
 580           g = near;
 581         }
 582       }
 583     }
 584   }
 585
 586   bool isOk = true;
 587   return isOk;
 588 }
 589
 590 //================================================================================
 591 /*!
 592  * \brief Return true if only two given edges meat at their common vertex
 593  */
 594 //================================================================================
 595
 596 static bool twoEdgesMeatAtVertex(const TopoDS_Edge& e1,
 597                                  const TopoDS_Edge& e2,
 598                                  SMESH_Mesh &       mesh)
 599 {
 600   TopoDS_Vertex v;
 601   if ( !TopExp::CommonVertex( e1, e2, v ))
 602     return false;
 603   TopTools_ListIteratorOfListOfShape ancestIt( mesh.GetAncestors( v ));
 604   for ( ; ancestIt.More() ; ancestIt.Next() )
 605     if ( ancestIt.Value().ShapeType() == TopAbs_EDGE )
 606       if ( !e1.IsSame( ancestIt.Value() ) && !e2.IsSame( ancestIt.Value() ))
 607         return false;
 608   return true;
 609 }
 610
 611 //=============================================================================
 612 /*!
 613  *
 614  */
 615 //=============================================================================
 616
 617 FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh &         aMesh,
 618                                                        const TopoDS_Shape & aShape)
 619   //throw(SALOME_Exception)
 620 {
 621   const TopoDS_Face & F = TopoDS::Face(aShape);
 622   const bool ignoreMediumNodes = _quadraticMesh;
 623
 624   // verify 1 wire only, with 4 edges
 625   TopoDS_Vertex V;
 626   list< TopoDS_Edge > edges;
 627   list< int > nbEdgesInWire;
 628   int nbWire = SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
 629   if (nbWire != 1) {
 630     error(COMPERR_BAD_SHAPE, TComm("Wrong number of wires: ") << nbWire);
 631     return 0;
 632   }
 633   FaceQuadStruct* quad = new FaceQuadStruct;
 634   quad->uv_grid = 0;
 635   quad->side.reserve(nbEdgesInWire.front());
 636
 637   int nbSides = 0;
 638   list< TopoDS_Edge >::iterator edgeIt = edges.begin();
 639   if ( nbEdgesInWire.front() == 4 ) { // exactly 4 edges
 640     for ( ; edgeIt != edges.end(); ++edgeIt, nbSides++ )
 641       quad->side.push_back( new StdMeshers_FaceSide(F, *edgeIt, &aMesh,
 642                                                     nbSides<TOP_SIDE, ignoreMediumNodes));
 643   }
 644   else if ( nbEdgesInWire.front() > 4 ) { // more than 4 edges - try to unite some
 645     list< TopoDS_Edge > sideEdges;
 646     while ( !edges.empty()) {
 647       sideEdges.clear();
 648       sideEdges.splice( sideEdges.end(), edges, edges.begin()); // edges.front() -> sideEdges.end()
 649       bool sameSide = true;
 650       while ( !edges.empty() && sameSide ) {
 651         sameSide = SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() );
 652         if ( sameSide )
 653           sideEdges.splice( sideEdges.end(), edges, edges.begin());
 654       }
 655       if ( nbSides == 0 ) { // go backward from the first edge
 656         sameSide = true;
 657         while ( !edges.empty() && sameSide ) {
 658           sameSide = SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() );
 659           if ( sameSide )
 660             sideEdges.splice( sideEdges.begin(), edges, --edges.end());
 661         }
 662       }
 663       quad->side.push_back( new StdMeshers_FaceSide(F, sideEdges, &aMesh,
 664                                                     nbSides<TOP_SIDE, ignoreMediumNodes));
 665       ++nbSides;
 666     }
 667     // issue 20222. Try to unite only edges shared by two same faces
 668     if (nbSides < 4) {
 669       // delete found sides
 670       { FaceQuadStruct cleaner( *quad ); }
 671       quad->side.clear();
 672       quad->side.reserve(nbEdgesInWire.front());
 673       nbSides = 0;
 674
 675       SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
 676       while ( !edges.empty()) {
 677         sideEdges.clear();
 678         sideEdges.splice( sideEdges.end(), edges, edges.begin());
 679         bool sameSide = true;
 680         while ( !edges.empty() && sameSide ) {
 681           sameSide =
 682             SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() ) &&
 683             twoEdgesMeatAtVertex( sideEdges.back(), edges.front(), aMesh );
 684           if ( sameSide )
 685             sideEdges.splice( sideEdges.end(), edges, edges.begin());
 686         }
 687         if ( nbSides == 0 ) { // go backward from the first edge
 688           sameSide = true;
 689           while ( !edges.empty() && sameSide ) {
 690             sameSide =
 691               SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() ) &&
 692               twoEdgesMeatAtVertex( sideEdges.front(), edges.back(), aMesh );
 693             if ( sameSide )
 694               sideEdges.splice( sideEdges.begin(), edges, --edges.end());
 695           }
 696         }
 697         quad->side.push_back( new StdMeshers_FaceSide(F, sideEdges, &aMesh,
 698                                                       nbSides<TOP_SIDE, ignoreMediumNodes));
 699         ++nbSides;
 700       }
 701     }
 702   }
 703   if (nbSides != 4) {
 704 #ifdef _DEBUG_
 705     MESSAGE ( "StdMeshers_Quadrangle_2D. Edge IDs of " << nbSides << " sides:\n" );
 706     for ( int i = 0; i < nbSides; ++i ) {
 707       MESSAGE ( " ( " );
 708       for ( int e = 0; e < quad->side[i]->NbEdges(); ++e )
 709         MESSAGE ( myTool->GetMeshDS()->ShapeToIndex( quad->side[i]->Edge( e )) << " " );
 710       MESSAGE ( ")\n" );
 711     }
 712     //cout << endl;
 713 #endif
 714     if ( !nbSides )
 715       nbSides = nbEdgesInWire.front();
 716     error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides);
 717     delete quad;
 718     quad = 0;
 719   }
 720
 721   return quad;
 722 }
 723
 724 //=============================================================================
 725 /*!
 726  *  CheckAnd2Dcompute
 727  */
 728 //=============================================================================
 729
 730 FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute
 731                            (SMESH_Mesh &         aMesh,
 732                             const TopoDS_Shape & aShape,
 733                             const bool           CreateQuadratic) //throw(SALOME_Exception)
 734 {
 735   _quadraticMesh = CreateQuadratic;
 736
 737   FaceQuadStruct *quad = CheckNbEdges(aMesh, aShape);
 738
 739   if(!quad) return 0;
 740
 741   // set normalized grid on unit square in parametric domain
 742   bool stat = SetNormalizedGrid(aMesh, aShape, quad);
 743   if(!stat) {
 744     if(!quad)
 745       delete quad;
 746     quad = 0;
 747   }
 748
 749   return quad;
 750 }
 751
 752 //=============================================================================
 753 /*!
 754  *
 755  */
 756 //=============================================================================
 757
 758 faceQuadStruct::~faceQuadStruct()
 759 {
 760   for (int i = 0; i < side.size(); i++) {
 761     if (side[i])     delete side[i];
 762   }
 763   if (uv_grid)       delete [] uv_grid;
 764 }
 765
 766 namespace {
 767   inline const vector<UVPtStruct>& GetUVPtStructIn(FaceQuadStruct* quad, int i, int nbSeg)
 768   {
 769     bool   isXConst   = ( i == BOTTOM_SIDE || i == TOP_SIDE );
 770     double constValue = ( i == BOTTOM_SIDE || i == LEFT_SIDE ) ? 0 : 1;
 771     return
 772       quad->isEdgeOut[i] ?
 773       quad->side[i]->SimulateUVPtStruct(nbSeg,isXConst,constValue) :
 774       quad->side[i]->GetUVPtStruct(isXConst,constValue);
 775   }
 776 }
 777
 778 //=============================================================================
 779 /*!
 780  *
 781  */
 782 //=============================================================================
 783
 784 bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh,
 785                                                   const TopoDS_Shape& aShape,
 786                                                   FaceQuadStruct* & quad) //throw (SALOME_Exception)
 787 {
 788   // Algorithme décrit dans "Génération automatique de maillages"
 789   // P.L. GEORGE, MASSON, § 6.4.1 p. 84-85
 790   // traitement dans le domaine paramétrique 2d u,v
 791   // transport - projection sur le carré unité
 792
 793 //  MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid");
 794 //  const TopoDS_Face& F = TopoDS::Face(aShape);
 795
 796   // 1 --- find orientation of the 4 edges, by test on extrema
 797
 798   //      max             min                    0     x1     1
 799   //     |<----north-2-------^                a3 -------------> a2
 800   //     |                   |                   ^1          1^
 801   //    west-3            east-1 =right          |            |
 802   //     |                   |         ==>       |            |
 803   //  y0 |                   | y1                |            |
 804   //     |                   |                   |0          0|
 805   //     v----south-0-------->                a0 -------------> a1
 806   //      min             max                    0     x0     1
 807   //             =down
 808   //
 809
 810   // 3 --- 2D normalized values on unit square [0..1][0..1]
 811
 812   int nbhoriz  = Min(quad->side[0]->NbPoints(), quad->side[2]->NbPoints());
 813   int nbvertic = Min(quad->side[1]->NbPoints(), quad->side[3]->NbPoints());
 814
 815   quad->isEdgeOut[0] = (quad->side[0]->NbPoints() > quad->side[2]->NbPoints());
 816   quad->isEdgeOut[1] = (quad->side[1]->NbPoints() > quad->side[3]->NbPoints());
 817   quad->isEdgeOut[2] = (quad->side[2]->NbPoints() > quad->side[0]->NbPoints());
 818   quad->isEdgeOut[3] = (quad->side[3]->NbPoints() > quad->side[1]->NbPoints());
 819
 820   UVPtStruct *uv_grid = quad->uv_grid = new UVPtStruct[nbvertic * nbhoriz];
 821
 822   const vector<UVPtStruct>& uv_e0 = GetUVPtStructIn( quad, 0, nbhoriz - 1 );
 823   const vector<UVPtStruct>& uv_e1 = GetUVPtStructIn( quad, 1, nbvertic - 1 );
 824   const vector<UVPtStruct>& uv_e2 = GetUVPtStructIn( quad, 2, nbhoriz - 1 );
 825   const vector<UVPtStruct>& uv_e3 = GetUVPtStructIn( quad, 3, nbvertic - 1 );
 826
 827   if ( uv_e0.empty() || uv_e1.empty() || uv_e2.empty() || uv_e3.empty() )
 828     //return error( "Can't find nodes on sides");
 829     return error( COMPERR_BAD_INPUT_MESH );
 830
 831   // nodes Id on "in" edges
 832   if (! quad->isEdgeOut[0]) {
 833     int j = 0;
 834     for (int i = 0; i < nbhoriz; i++) { // down
 835       int ij = j * nbhoriz + i;
 836       uv_grid[ij].node = uv_e0[i].node;
 837     }
 838   }
 839   if (! quad->isEdgeOut[1]) {
 840     int i = nbhoriz - 1;
 841     for (int j = 0; j < nbvertic; j++) { // right
 842       int ij = j * nbhoriz + i;
 843       uv_grid[ij].node = uv_e1[j].node;
 844     }
 845   }
 846   if (! quad->isEdgeOut[2]) {
 847     int j = nbvertic - 1;
 848     for (int i = 0; i < nbhoriz; i++) { // up
 849       int ij = j * nbhoriz + i;
 850       uv_grid[ij].node = uv_e2[i].node;
 851     }
 852   }
 853   if (! quad->isEdgeOut[3]) {
 854     int i = 0;
 855     for (int j = 0; j < nbvertic; j++) { // left
 856       int ij = j * nbhoriz + i;
 857       uv_grid[ij].node = uv_e3[j].node;
 858     }
 859   }
 860
 861   // normalized 2d values on grid
 862   for (int i = 0; i < nbhoriz; i++)
 863   {
 864     for (int j = 0; j < nbvertic; j++)
 865     {
 866       int ij = j * nbhoriz + i;
 867       // --- droite i cste : x = x0 + y(x1-x0)
 868       double x0 = uv_e0[i].normParam;   // bas - sud
 869       double x1 = uv_e2[i].normParam;   // haut - nord
 870       // --- droite j cste : y = y0 + x(y1-y0)
 871       double y0 = uv_e3[j].normParam;   // gauche-ouest
 872       double y1 = uv_e1[j].normParam;   // droite - est
 873       // --- intersection : x=x0+(y0+x(y1-y0))(x1-x0)
 874       double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
 875       double y = y0 + x * (y1 - y0);
 876       uv_grid[ij].x = x;
 877       uv_grid[ij].y = y;
 878       //MESSAGE("-xy-01 "<<x0<<" "<<x1<<" "<<y0<<" "<<y1);
 879       //MESSAGE("-xy-norm "<<i<<" "<<j<<" "<<x<<" "<<y);
 880     }
 881   }
 882
 883   // 4 --- projection on 2d domain (u,v)
 884   gp_UV a0( uv_e0.front().u, uv_e0.front().v );
 885   gp_UV a1( uv_e0.back().u,  uv_e0.back().v );
 886   gp_UV a2( uv_e2.back().u,  uv_e2.back().v );
 887   gp_UV a3( uv_e2.front().u, uv_e2.front().v );
 888
 889   for (int i = 0; i < nbhoriz; i++)
 890   {
 891     for (int j = 0; j < nbvertic; j++)
 892     {
 893       int ij = j * nbhoriz + i;
 894       double x = uv_grid[ij].x;
 895       double y = uv_grid[ij].y;
 896       double param_0 = uv_e0[0].normParam + x * (uv_e0.back().normParam - uv_e0[0].normParam); // sud
 897       double param_2 = uv_e2[0].normParam + x * (uv_e2.back().normParam - uv_e2[0].normParam); // nord
 898       double param_1 = uv_e1[0].normParam + y * (uv_e1.back().normParam - uv_e1[0].normParam); // est
 899       double param_3 = uv_e3[0].normParam + y * (uv_e3.back().normParam - uv_e3[0].normParam); // ouest
 900
 901       //MESSAGE("params "<<param_0<<" "<<param_1<<" "<<param_2<<" "<<param_3);
 902       gp_UV p0 = quad->side[0]->Value2d(param_0).XY();
 903       gp_UV p1 = quad->side[1]->Value2d(param_1).XY();
 904       gp_UV p2 = quad->side[2]->Value2d(param_2).XY();
 905       gp_UV p3 = quad->side[3]->Value2d(param_3).XY();
 906
 907       gp_UV uv = (1 - y) * p0 + x * p1 + y * p2 + (1 - x) * p3;
 908       uv -= (1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3;
 909
 910       uv_grid[ij].u = uv.X();
 911       uv_grid[ij].v = uv.Y();
 912     }
 913   }
 914   return true;
 915 }
 916
 917 //=======================================================================
 918 //function : ShiftQuad
 919 //purpose  : auxilary function for ComputeQuadPref
 920 //=======================================================================
 921
 922 static void ShiftQuad(FaceQuadStruct* quad, const int num, bool)
 923 {
 924   StdMeshers_FaceSide* side[4] = { quad->side[0], quad->side[1], quad->side[2], quad->side[3] };
 925   for (int i = BOTTOM_SIDE; i < NB_SIDES; ++i ) {
 926     int id = ( i + num ) % NB_SIDES;
 927     bool wasForward = ( i < TOP_SIDE );
 928     bool newForward = ( id < TOP_SIDE );
 929     if ( wasForward != newForward )
 930       side[ i ]->Reverse();
 931     quad->side[ id ] = side[ i ];
 932   }
 933 }
 934
 935 //=======================================================================
 936 //function : CalcUV
 937 //purpose  : auxilary function for ComputeQuadPref
 938 //=======================================================================
 939
 940 static gp_UV CalcUV(double x0, double x1, double y0, double y1,
 941                     FaceQuadStruct* quad,
 942                     const gp_UV& a0, const gp_UV& a1,
 943                     const gp_UV& a2, const gp_UV& a3)
 944 {
 945   const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0 );
 946   const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
 947   const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1 );
 948   const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
 949
 950   double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
 951   double y = y0 + x * (y1 - y0);
 952
 953   double param_b = uv_eb[0].normParam + x * (uv_eb.back().normParam - uv_eb[0].normParam);
 954   double param_t = uv_et[0].normParam + x * (uv_et.back().normParam - uv_et[0].normParam);
 955   double param_r = uv_er[0].normParam + y * (uv_er.back().normParam - uv_er[0].normParam);
 956   double param_l = uv_el[0].normParam + y * (uv_el.back().normParam - uv_el[0].normParam);
 957
 958   gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(param_b).XY();
 959   gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(param_r).XY();
 960   gp_UV p2 = quad->side[TOP_SIDE   ]->Value2d(param_t).XY();
 961   gp_UV p3 = quad->side[LEFT_SIDE  ]->Value2d(param_l).XY();
 962
 963   gp_UV uv = p0 * (1 - y) + p1 * x + p2 * y + p3 * (1 - x);
 964
 965   uv -= (1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3;
 966
 967   return uv;
 968 }
 969
 970 //=======================================================================
 971 //function : CalcUV2
 972 //purpose  : auxilary function for ComputeQuadPref
 973 //=======================================================================
 974
 975 static gp_UV CalcUV2(double x, double y,
 976                      FaceQuadStruct* quad,
 977                      const gp_UV& a0, const gp_UV& a1,
 978                      const gp_UV& a2, const gp_UV& a3)
 979 {
 980   const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0 );
 981   const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
 982   const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1 );
 983   const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
 984
 985   //double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
 986   //double y = y0 + x * (y1 - y0);
 987
 988   double param_b = uv_eb[0].normParam + x * (uv_eb.back().normParam - uv_eb[0].normParam);
 989   double param_t = uv_et[0].normParam + x * (uv_et.back().normParam - uv_et[0].normParam);
 990   double param_r = uv_er[0].normParam + y * (uv_er.back().normParam - uv_er[0].normParam);
 991   double param_l = uv_el[0].normParam + y * (uv_el.back().normParam - uv_el[0].normParam);
 992
 993   gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(param_b).XY();
 994   gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(param_r).XY();
 995   gp_UV p2 = quad->side[TOP_SIDE   ]->Value2d(param_t).XY();
 996   gp_UV p3 = quad->side[LEFT_SIDE  ]->Value2d(param_l).XY();
 997
 998   gp_UV uv = p0 * (1 - y) + p1 * x + p2 * y + p3 * (1 - x);
 999
1000   uv -= (1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3;
1001
1002   return uv;
1003 }
1004
1005
1006 //=======================================================================
1007 /*!
1008  * Create only quandrangle faces
1009  */
1010 //=======================================================================
1011
1012 bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh &        aMesh,
1013                                                 const TopoDS_Shape& aShape,
1014                                                 FaceQuadStruct*     quad)
1015 {
1016   // Auxilary key in order to keep old variant
1017   // of meshing after implementation new variant
1018   // for bug 0016220 from Mantis.
1019   bool OldVersion = false;
1020
1021   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
1022   const TopoDS_Face& F = TopoDS::Face(aShape);
1023   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
1024 //  const TopoDS_Wire& W = BRepTools::OuterWire(F);
1025   bool WisF = true;
1026 //   if(W.Orientation()==TopAbs_FORWARD)
1027 //     WisF = true;
1028   //if(WisF) cout<<"W is FORWARD"<<endl;
1029   //else cout<<"W is REVERSED"<<endl;
1030 //   bool FisF = (F.Orientation()==TopAbs_FORWARD);
1031 //   if(!FisF) WisF = !WisF;
1032 //  WisF = FisF;
1033   int i,j,geomFaceID = meshDS->ShapeToIndex( F );
1034
1035   int nb = quad->side[0]->NbPoints();
1036   int nr = quad->side[1]->NbPoints();
1037   int nt = quad->side[2]->NbPoints();
1038   int nl = quad->side[3]->NbPoints();
1039   int dh = abs(nb-nt);
1040   int dv = abs(nr-nl);
1041
1042   if( dh>=dv ) {
1043     if( nt>nb ) {
1044       // it is a base case => not shift quad but me be replacement is need
1045       ShiftQuad(quad,0,WisF);
1046     }
1047     else {
1048       // we have to shift quad on 2
1049       ShiftQuad(quad,2,WisF);
1050     }
1051   }
1052   else {
1053     if( nr>nl ) {
1054       // we have to shift quad on 1
1055       ShiftQuad(quad,1,WisF);
1056     }
1057     else {
1058       // we have to shift quad on 3
1059       ShiftQuad(quad,3,WisF);
1060     }
1061   }
1062
1063   nb = quad->side[0]->NbPoints();
1064   nr = quad->side[1]->NbPoints();
1065   nt = quad->side[2]->NbPoints();
1066   nl = quad->side[3]->NbPoints();
1067   dh = abs(nb-nt);
1068   dv = abs(nr-nl);
1069   int nbh  = Max(nb,nt);
1070   int nbv = Max(nr,nl);
1071   int addh = 0;
1072   int addv = 0;
1073
1074   // ----------- Old version ---------------
1075   // orientation of face and 3 main domain for future faces
1076   //       0   top    1
1077   //      1------------1
1078   //       |   |  |   |
1079   //       |   |  |   |
1080   //       | L |  | R |
1081   //  left |   |  |   | rigth
1082   //       |  /    \  |
1083   //       | /  C   \ |
1084   //       |/        \|
1085   //      0------------0
1086   //       0  bottom  1
1087
1088   // ----------- New version ---------------
1089   // orientation of face and 3 main domain for future faces
1090   //       0   top    1
1091   //      1------------1
1092   //       |  |____|  |
1093   //       |  /    \  |
1094   //       | /  C   \ |
1095   //  left |/________\| rigth
1096   //       |          |
1097   //       |          |
1098   //       |          |
1099   //      0------------0
1100   //       0  bottom  1
1101
1102   if(dh>dv) {
1103     addv = (dh-dv)/2;
1104     nbv = nbv + addv;
1105   }
1106   else { // dv>=dh
1107     addh = (dv-dh)/2;
1108     nbh = nbh + addh;
1109   }
1110
1111   const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0 );
1112   const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
1113   const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1 );
1114   const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
1115
1116   // arrays for normalized params
1117   //cout<<"Dump B:"<<endl;
1118   TColStd_SequenceOfReal npb, npr, npt, npl;
1119   for(i=0; i<nb; i++) {
1120     npb.Append(uv_eb[i].normParam);
1121     //cout<<"i="<<i<<" par="<<uv_eb[i].normParam<<" npar="<<uv_eb[i].normParam;
1122     //const SMDS_MeshNode* N = uv_eb[i].node;
1123     //cout<<" node("<<N->X()<<","<<N->Y()<<","<<N->Z()<<")"<<endl;
1124   }
1125   for(i=0; i<nr; i++) {
1126     npr.Append(uv_er[i].normParam);
1127   }
1128   for(i=0; i<nt; i++) {
1129     npt.Append(uv_et[i].normParam);
1130   }
1131   for(i=0; i<nl; i++) {
1132     npl.Append(uv_el[i].normParam);
1133   }
1134
1135   int dl,dr;
1136   if(OldVersion) {
1137     // add some params to right and left after the first param
1138     // insert to right
1139     dr = nbv - nr;
1140     double dpr = (npr.Value(2) - npr.Value(1))/(dr+1);
1141     for(i=1; i<=dr; i++) {
1142       npr.InsertAfter(1,npr.Value(2)-dpr);
1143     }
1144     // insert to left
1145     dl = nbv - nl;
1146     dpr = (npl.Value(2) - npl.Value(1))/(dl+1);
1147     for(i=1; i<=dl; i++) {
1148       npl.InsertAfter(1,npl.Value(2)-dpr);
1149     }
1150   }
1151   //cout<<"npb:";
1152   //for(i=1; i<=npb.Length(); i++) {
1153   //  cout<<" "<<npb.Value(i);
1154   //}
1155   //cout<<endl;
1156
1157   gp_XY a0( uv_eb.front().u, uv_eb.front().v );
1158   gp_XY a1( uv_eb.back().u,  uv_eb.back().v );
1159   gp_XY a2( uv_et.back().u,  uv_et.back().v );
1160   gp_XY a3( uv_et.front().u, uv_et.front().v );
1161   //cout<<" a0("<<a0.X()<<","<<a0.Y()<<")"<<" a1("<<a1.X()<<","<<a1.Y()<<")"
1162   //    <<" a2("<<a2.X()<<","<<a2.Y()<<")"<<" a3("<<a3.X()<<","<<a3.Y()<<")"<<endl;
1163
1164   int nnn = Min(nr,nl);
1165   // auxilary sequence of XY for creation nodes
1166   // in the bottom part of central domain
1167   // it's length must be == nbv-nnn-1
1168   TColgp_SequenceOfXY UVL;
1169   TColgp_SequenceOfXY UVR;
1170
1171   if(OldVersion) {
1172     // step1: create faces for left domain
1173     StdMeshers_Array2OfNode NodesL(1,dl+1,1,nl);
1174     // add left nodes
1175     for(j=1; j<=nl; j++)
1176       NodesL.SetValue(1,j,uv_el[j-1].node);
1177     if(dl>0) {
1178       // add top nodes
1179       for(i=1; i<=dl; i++)
1180         NodesL.SetValue(i+1,nl,uv_et[i].node);
1181       // create and add needed nodes
1182       TColgp_SequenceOfXY UVtmp;
1183       for(i=1; i<=dl; i++) {
1184         double x0 = npt.Value(i+1);
1185         double x1 = x0;
1186         // diagonal node
1187         double y0 = npl.Value(i+1);
1188         double y1 = npr.Value(i+1);
1189         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1190         gp_Pnt P = S->Value(UV.X(),UV.Y());
1191         SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1192         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1193         NodesL.SetValue(i+1,1,N);
1194         if(UVL.Length()<nbv-nnn-1) UVL.Append(UV);
1195         // internal nodes
1196         for(j=2; j<nl; j++) {
1197           double y0 = npl.Value(dl+j);
1198           double y1 = npr.Value(dl+j);
1199           gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1200           gp_Pnt P = S->Value(UV.X(),UV.Y());
1201           SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1202           meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1203           NodesL.SetValue(i+1,j,N);
1204           if( i==dl ) UVtmp.Append(UV);
1205         }
1206       }
1207       for(i=1; i<=UVtmp.Length() && UVL.Length()<nbv-nnn-1; i++) {
1208         UVL.Append(UVtmp.Value(i));
1209       }
1210       //cout<<"Dump NodesL:"<<endl;
1211       //for(i=1; i<=dl+1; i++) {
1212       //  cout<<"i="<<i;
1213       //  for(j=1; j<=nl; j++) {
1214       //    cout<<" ("<<NodesL.Value(i,j)->X()<<","<<NodesL.Value(i,j)->Y()<<","<<NodesL.Value(i,j)->Z()<<")";
1215       //  }
1216       //  cout<<endl;
1217       //}
1218       // create faces
1219       for(i=1; i<=dl; i++) {
1220         for(j=1; j<nl; j++) {
1221           if(WisF) {
1222             SMDS_MeshFace* F =
1223               myTool->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j),
1224                               NodesL.Value(i+1,j+1), NodesL.Value(i,j+1));
1225             meshDS->SetMeshElementOnShape(F, geomFaceID);
1226           }
1227           else {
1228             SMDS_MeshFace* F =
1229               myTool->AddFace(NodesL.Value(i,j), NodesL.Value(i,j+1),
1230                               NodesL.Value(i+1,j+1), NodesL.Value(i+1,j));
1231             meshDS->SetMeshElementOnShape(F, geomFaceID);
1232           }
1233         }
1234       }
1235     }
1236     else {
1237       // fill UVL using c2d
1238       for(i=1; i<npl.Length() && UVL.Length()<nbv-nnn-1; i++) {
1239         UVL.Append( gp_UV ( uv_el[i].u, uv_el[i].v ));
1240       }
1241     }
1242
1243     // step2: create faces for right domain
1244     StdMeshers_Array2OfNode NodesR(1,dr+1,1,nr);
1245     // add right nodes
1246     for(j=1; j<=nr; j++)
1247       NodesR.SetValue(1,j,uv_er[nr-j].node);
1248     if(dr>0) {
1249       // add top nodes
1250       for(i=1; i<=dr; i++)
1251         NodesR.SetValue(i+1,1,uv_et[nt-1-i].node);
1252       // create and add needed nodes
1253       TColgp_SequenceOfXY UVtmp;
1254       for(i=1; i<=dr; i++) {
1255         double x0 = npt.Value(nt-i);
1256         double x1 = x0;
1257         // diagonal node
1258         double y0 = npl.Value(i+1);
1259         double y1 = npr.Value(i+1);
1260         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1261         gp_Pnt P = S->Value(UV.X(),UV.Y());
1262         SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1263         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1264         NodesR.SetValue(i+1,nr,N);
1265         if(UVR.Length()<nbv-nnn-1) UVR.Append(UV);
1266         // internal nodes
1267         for(j=2; j<nr; j++) {
1268           double y0 = npl.Value(nbv-j+1);
1269           double y1 = npr.Value(nbv-j+1);
1270           gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1271           gp_Pnt P = S->Value(UV.X(),UV.Y());
1272           SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1273           meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1274           NodesR.SetValue(i+1,j,N);
1275           if( i==dr ) UVtmp.Prepend(UV);
1276         }
1277       }
1278       for(i=1; i<=UVtmp.Length() && UVR.Length()<nbv-nnn-1; i++) {
1279         UVR.Append(UVtmp.Value(i));
1280       }
1281       // create faces
1282       for(i=1; i<=dr; i++) {
1283         for(j=1; j<nr; j++) {
1284           if(WisF) {
1285             SMDS_MeshFace* F =
1286               myTool->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j),
1287                               NodesR.Value(i+1,j+1), NodesR.Value(i,j+1));
1288             meshDS->SetMeshElementOnShape(F, geomFaceID);
1289           }
1290           else {
1291             SMDS_MeshFace* F =
1292               myTool->AddFace(NodesR.Value(i,j), NodesR.Value(i,j+1),
1293                               NodesR.Value(i+1,j+1), NodesR.Value(i+1,j));
1294             meshDS->SetMeshElementOnShape(F, geomFaceID);
1295           }
1296         }
1297       }
1298     }
1299     else {
1300       // fill UVR using c2d
1301       for(i=1; i<npr.Length() && UVR.Length()<nbv-nnn-1; i++) {
1302         UVR.Append( gp_UV( uv_er[i].u, uv_er[i].v ));
1303       }
1304     }
1305
1306     // step3: create faces for central domain
1307     StdMeshers_Array2OfNode NodesC(1,nb,1,nbv);
1308     // add first string using NodesL
1309     for(i=1; i<=dl+1; i++)
1310       NodesC.SetValue(1,i,NodesL(i,1));
1311     for(i=2; i<=nl; i++)
1312       NodesC.SetValue(1,dl+i,NodesL(dl+1,i));
1313     // add last string using NodesR
1314     for(i=1; i<=dr+1; i++)
1315       NodesC.SetValue(nb,i,NodesR(i,nr));
1316     for(i=1; i<nr; i++)
1317       NodesC.SetValue(nb,dr+i+1,NodesR(dr+1,nr-i));
1318     // add top nodes (last columns)
1319     for(i=dl+2; i<nbh-dr; i++)
1320       NodesC.SetValue(i-dl,nbv,uv_et[i-1].node);
1321     // add bottom nodes (first columns)
1322     for(i=2; i<nb; i++)
1323       NodesC.SetValue(i,1,uv_eb[i-1].node);
1324
1325     // create and add needed nodes
1326     // add linear layers
1327     for(i=2; i<nb; i++) {
1328       double x0 = npt.Value(dl+i);
1329       double x1 = x0;
1330       for(j=1; j<nnn; j++) {
1331         double y0 = npl.Value(nbv-nnn+j);
1332         double y1 = npr.Value(nbv-nnn+j);
1333         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1334         gp_Pnt P = S->Value(UV.X(),UV.Y());
1335         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1336         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1337         NodesC.SetValue(i,nbv-nnn+j,N);
1338       }
1339     }
1340     // add diagonal layers
1341     //cout<<"UVL.Length()="<<UVL.Length()<<" UVR.Length()="<<UVR.Length()<<endl;
1342     //cout<<"Dump UVL:"<<endl;
1343     //for(i=1; i<=UVL.Length(); i++) {
1344     //  cout<<" ("<<UVL.Value(i).X()<<","<<UVL.Value(i).Y()<<")";
1345     //}
1346     //cout<<endl;
1347     for(i=1; i<nbv-nnn; i++) {
1348       double du = UVR.Value(i).X() - UVL.Value(i).X();
1349       double dv = UVR.Value(i).Y() - UVL.Value(i).Y();
1350       for(j=2; j<nb; j++) {
1351         double u = UVL.Value(i).X() + du*npb.Value(j);
1352         double v = UVL.Value(i).Y() + dv*npb.Value(j);
1353         gp_Pnt P = S->Value(u,v);
1354         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1355         meshDS->SetNodeOnFace(N, geomFaceID, u, v);
1356         NodesC.SetValue(j,i+1,N);
1357       }
1358     }
1359     // create faces
1360     for(i=1; i<nb; i++) {
1361       for(j=1; j<nbv; j++) {
1362         if(WisF) {
1363           SMDS_MeshFace* F =
1364             myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
1365                             NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
1366           meshDS->SetMeshElementOnShape(F, geomFaceID);
1367         }
1368         else {
1369           SMDS_MeshFace* F =
1370             myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1),
1371                             NodesC.Value(i+1,j+1), NodesC.Value(i+1,j));
1372           meshDS->SetMeshElementOnShape(F, geomFaceID);
1373         }
1374       }
1375     }
1376   }
1377
1378   else { // New version (!OldVersion)
1379     // step1: create faces for bottom rectangle domain
1380     StdMeshers_Array2OfNode NodesBRD(1,nb,1,nnn-1);
1381     // fill UVL and UVR using c2d
1382     for(j=0; j<nb; j++) {
1383       NodesBRD.SetValue(j+1,1,uv_eb[j].node);
1384     }
1385     for(i=1; i<nnn-1; i++) {
1386       NodesBRD.SetValue(1,i+1,uv_el[i].node);
1387       NodesBRD.SetValue(nb,i+1,uv_er[i].node);
1388       double du = uv_er[i].u - uv_el[i].u;
1389       double dv = uv_er[i].v - uv_el[i].v;
1390       for(j=2; j<nb; j++) {
1391         double u = uv_el[i].u + du*npb.Value(j);
1392         double v = uv_el[i].v + dv*npb.Value(j);
1393         gp_Pnt P = S->Value(u,v);
1394         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1395         meshDS->SetNodeOnFace(N, geomFaceID, u, v);
1396         NodesBRD.SetValue(j,i+1,N);
1397
1398       }
1399     }
1400     for(j=1; j<nnn-1; j++) {
1401       for(i=1; i<nb; i++) {
1402         if(WisF) {
1403           SMDS_MeshFace* F =
1404             myTool->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i+1,j),
1405                             NodesBRD.Value(i+1,j+1), NodesBRD.Value(i,j+1));
1406           meshDS->SetMeshElementOnShape(F, geomFaceID);
1407         }
1408         else {
1409           SMDS_MeshFace* F =
1410             myTool->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i,j+1),
1411                             NodesBRD.Value(i+1,j+1), NodesBRD.Value(i+1,j));
1412           meshDS->SetMeshElementOnShape(F, geomFaceID);
1413         }
1414       }
1415     }
1416
1417     int drl = abs(nr-nl);
1418     // create faces for region C
1419     StdMeshers_Array2OfNode NodesC(1,nb,1,drl+1+addv);
1420     // add nodes from previous region
1421     for(j=1; j<=nb; j++) {
1422       NodesC.SetValue(j,1,NodesBRD.Value(j,nnn-1));
1423     }
1424     if( (drl+addv) > 0 ) {
1425       int n1,n2;
1426       if(nr>nl) {
1427         n1 = 1;
1428         n2 = drl + 1;
1429         TColgp_SequenceOfXY UVtmp;
1430         double drparam = npr.Value(nr) - npr.Value(nnn-1);
1431         double dlparam = npl.Value(nnn) - npl.Value(nnn-1);
1432         double y0,y1;
1433         for(i=1; i<=drl; i++) {
1434           // add existed nodes from right edge
1435           NodesC.SetValue(nb,i+1,uv_er[nnn+i-2].node);
1436           //double dtparam = npt.Value(i+1);
1437           y1 = npr.Value(nnn+i-1); // param on right edge
1438           double dpar = (y1 - npr.Value(nnn-1))/drparam;
1439           y0 = npl.Value(nnn-1) + dpar*dlparam; // param on left edge
1440           double dy = y1 - y0;
1441           for(j=1; j<nb; j++) {
1442             double x = npt.Value(i+1) + npb.Value(j)*(1-npt.Value(i+1));
1443             double y = y0 + dy*x;
1444             gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1445             gp_Pnt P = S->Value(UV.X(),UV.Y());
1446             SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1447             meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1448             NodesC.SetValue(j,i+1,N);
1449           }
1450         }
1451         double dy0 = (1-y0)/(addv+1);
1452         double dy1 = (1-y1)/(addv+1);
1453         for(i=1; i<=addv; i++) {
1454           double yy0 = y0 + dy0*i;
1455           double yy1 = y1 + dy1*i;
1456           double dyy = yy1 - yy0;
1457           for(j=1; j<=nb; j++) {
1458             double x = npt.Value(i+1+drl) +
1459               npb.Value(j) * ( npt.Value(nt-i) - npt.Value(i+1+drl) );
1460             double y = yy0 + dyy*x;
1461             gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1462             gp_Pnt P = S->Value(UV.X(),UV.Y());
1463             SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1464             meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1465             NodesC.SetValue(j,i+drl+1,N);
1466           }
1467         }
1468       }
1469       else { // nr<nl
1470         n2 = 1;
1471         n1 = drl + 1;
1472         TColgp_SequenceOfXY UVtmp;
1473         double dlparam = npl.Value(nl) - npl.Value(nnn-1);
1474         double drparam = npr.Value(nnn) - npr.Value(nnn-1);
1475         double y0 = npl.Value(nnn-1);
1476         double y1 = npr.Value(nnn-1);
1477         for(i=1; i<=drl; i++) {
1478           // add existed nodes from right edge
1479           NodesC.SetValue(1,i+1,uv_el[nnn+i-2].node);
1480           y0 = npl.Value(nnn+i-1); // param on left edge
1481           double dpar = (y0 - npl.Value(nnn-1))/dlparam;
1482           y1 = npr.Value(nnn-1) + dpar*drparam; // param on right edge
1483           double dy = y1 - y0;
1484           for(j=2; j<=nb; j++) {
1485             double x = npb.Value(j)*npt.Value(nt-i);
1486             double y = y0 + dy*x;
1487             gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1488             gp_Pnt P = S->Value(UV.X(),UV.Y());
1489             SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1490             meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1491             NodesC.SetValue(j,i+1,N);
1492           }
1493         }
1494         double dy0 = (1-y0)/(addv+1);
1495         double dy1 = (1-y1)/(addv+1);
1496         for(i=1; i<=addv; i++) {
1497           double yy0 = y0 + dy0*i;
1498           double yy1 = y1 + dy1*i;
1499           double dyy = yy1 - yy0;
1500           for(j=1; j<=nb; j++) {
1501             double x = npt.Value(i+1) +
1502               npb.Value(j) * ( npt.Value(nt-i-drl) - npt.Value(i+1) );
1503             double y = yy0 + dyy*x;
1504             gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1505             gp_Pnt P = S->Value(UV.X(),UV.Y());
1506             SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1507             meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1508             NodesC.SetValue(j,i+drl+1,N);
1509           }
1510         }
1511       }
1512       // create faces
1513       for(j=1; j<=drl+addv; j++) {
1514         for(i=1; i<nb; i++) {
1515           if(WisF) {
1516             SMDS_MeshFace* F =
1517               myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
1518                               NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
1519             meshDS->SetMeshElementOnShape(F, geomFaceID);
1520           }
1521           else {
1522             SMDS_MeshFace* F =
1523               myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1),
1524                               NodesC.Value(i+1,j+1), NodesC.Value(i+1,j));
1525             meshDS->SetMeshElementOnShape(F, geomFaceID);
1526           }
1527         }
1528       } // end nr<nl
1529
1530       StdMeshers_Array2OfNode NodesLast(1,nt,1,2);
1531       for(i=1; i<=nt; i++) {
1532         NodesLast.SetValue(i,2,uv_et[i-1].node);
1533       }
1534       int nnn=0;
1535       for(i=n1; i<drl+addv+1; i++) {
1536         nnn++;
1537         NodesLast.SetValue(nnn,1,NodesC.Value(1,i));
1538       }
1539       for(i=1; i<=nb; i++) {
1540         nnn++;
1541         NodesLast.SetValue(nnn,1,NodesC.Value(i,drl+addv+1));
1542       }
1543       for(i=drl+addv; i>=n2; i--) {
1544         nnn++;
1545         NodesLast.SetValue(nnn,1,NodesC.Value(nb,i));
1546       }
1547       for(i=1; i<nt; i++) {
1548         if(WisF) {
1549           SMDS_MeshFace* F =
1550             myTool->AddFace(NodesLast.Value(i,1), NodesLast.Value(i+1,1),
1551                             NodesLast.Value(i+1,2), NodesLast.Value(i,2));
1552           meshDS->SetMeshElementOnShape(F, geomFaceID);
1553         }
1554         else {
1555           SMDS_MeshFace* F =
1556             myTool->AddFace(NodesLast.Value(i,1), NodesLast.Value(i,2),
1557                             NodesLast.Value(i+1,2), NodesLast.Value(i+1,2));
1558           meshDS->SetMeshElementOnShape(F, geomFaceID);
1559         }
1560       }
1561     } // if( (drl+addv) > 0 )
1562
1563   } // end new version implementation
1564
1565   bool isOk = true;
1566   return isOk;
1567 }
1568
1569 //=============================================================================
1570 /*! Split quadrangle in to 2 triangles by smallest diagonal
1571  *
1572  */
1573 //=============================================================================
1574 void StdMeshers_Quadrangle_2D::SplitQuad(SMESHDS_Mesh *theMeshDS,
1575                                     int theFaceID,
1576                                     const SMDS_MeshNode* theNode1,
1577                                     const SMDS_MeshNode* theNode2,
1578                                     const SMDS_MeshNode* theNode3,
1579                                     const SMDS_MeshNode* theNode4)
1580 {
1581   gp_Pnt a(theNode1->X(),theNode1->Y(),theNode1->Z());
1582   gp_Pnt b(theNode2->X(),theNode2->Y(),theNode2->Z());
1583   gp_Pnt c(theNode3->X(),theNode3->Y(),theNode3->Z());
1584   gp_Pnt d(theNode4->X(),theNode4->Y(),theNode4->Z());
1585   SMDS_MeshFace* face;
1586   if(a.Distance(c) > b.Distance(d)){
1587     face = myTool->AddFace(theNode2, theNode4 , theNode1);
1588     theMeshDS->SetMeshElementOnShape(face, theFaceID );
1589     face = myTool->AddFace(theNode2, theNode3, theNode4);
1590     theMeshDS->SetMeshElementOnShape(face, theFaceID );
1591
1592   }
1593   else{
1594     face = myTool->AddFace(theNode1, theNode2 ,theNode3);
1595     theMeshDS->SetMeshElementOnShape(face, theFaceID );
1596     face = myTool->AddFace(theNode1, theNode3, theNode4);
1597     theMeshDS->SetMeshElementOnShape(face, theFaceID );
1598   }
1599 }