src/StdMeshers/StdMeshers_Quadrangle_2D.cxx

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