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 /*!
 593  *  Evaluate
 594  */
 595 //=============================================================================
 596
 597 bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh& aMesh,
 598                                         const TopoDS_Shape& aShape,
 599                                         MapShapeNbElems& aResMap)
 600
 601 {
 602   aMesh.GetSubMesh(aShape);
 603
 604   std::vector<int> aNbNodes(4);
 605   bool IsQuadratic = false;
 606   if( !CheckNbEdgesForEvaluate( aMesh, aShape, aResMap, aNbNodes, IsQuadratic ) ) {
 607     std::vector<int> aResVec(17);
 608     for(int i=0; i<17; i++) aResVec[i] = 0;
 609     SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 610     aResMap.insert(std::make_pair(sm,aResVec));
 611     SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 612     smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
 613     return false;
 614   }
 615
 616   if(myQuadranglePreference) {
 617     int n1 = aNbNodes[0];
 618     int n2 = aNbNodes[1];
 619     int n3 = aNbNodes[2];
 620     int n4 = aNbNodes[3];
 621     int nfull = n1+n2+n3+n4;
 622     int ntmp = nfull/2;
 623     ntmp = ntmp*2;
 624     if( nfull==ntmp && ( (n1!=n3) || (n2!=n4) ) ) {
 625       // special path for using only quandrangle faces
 626       return EvaluateQuadPref(aMesh, aShape, aNbNodes, aResMap, IsQuadratic);
 627       //return true;
 628     }
 629   }
 630
 631   int nbdown  = aNbNodes[0];
 632   int nbup    = aNbNodes[2];
 633
 634   int nbright = aNbNodes[1];
 635   int nbleft  = aNbNodes[3];
 636
 637   int nbhoriz  = Min(nbdown, nbup);
 638   int nbvertic = Min(nbright, nbleft);
 639
 640   int dh = Max(nbdown, nbup) - nbhoriz;
 641   int dv = Max(nbright, nbleft) - nbvertic;
 642
 643   int kdh = 0;
 644   if(dh>0) kdh = 1;
 645   int kdv = 0;
 646   if(dv>0) kdv = 1;
 647
 648   int nbNodes = (nbhoriz-2)*(nbvertic-2);
 649   int nbFaces3 = dh + dv + kdh*(nbvertic-1)*2 + kdv*(nbhoriz-1)*2;
 650   if( kdh==1 && kdv==1 ) nbFaces3 -= 2;
 651   int nbFaces4 = (nbhoriz-1-kdh)*(nbvertic-1-kdv);
 652
 653   std::vector<int> aVec(17);
 654   for(int i=0; i<17; i++) aVec[i] = 0;
 655   if(IsQuadratic) {
 656     aVec[4] = nbFaces3;
 657     aVec[6] = nbFaces4;
 658     int nbbndedges = nbdown + nbup + nbright + nbleft -4;
 659     int nbintedges = ( nbFaces4*4 + nbFaces3*3 - nbbndedges ) / 2;
 660     aVec[0] = nbNodes + nbintedges;
 661   }
 662   else {
 663     aVec[0] = nbNodes;
 664     aVec[3] = nbFaces3;
 665     aVec[5] = nbFaces4;
 666   }
 667   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 668   aResMap.insert(std::make_pair(sm,aVec));
 669
 670   return true;
 671 }
 672
 673
 674 //================================================================================
 675 /*!
 676  * \brief Return true if only two given edges meat at their common vertex
 677  */
 678 //================================================================================
 679
 680 static bool twoEdgesMeatAtVertex(const TopoDS_Edge& e1,
 681                                  const TopoDS_Edge& e2,
 682                                  SMESH_Mesh &       mesh)
 683 {
 684   TopoDS_Vertex v;
 685   if ( !TopExp::CommonVertex( e1, e2, v ))
 686     return false;
 687   TopTools_ListIteratorOfListOfShape ancestIt( mesh.GetAncestors( v ));
 688   for ( ; ancestIt.More() ; ancestIt.Next() )
 689     if ( ancestIt.Value().ShapeType() == TopAbs_EDGE )
 690       if ( !e1.IsSame( ancestIt.Value() ) && !e2.IsSame( ancestIt.Value() ))
 691         return false;
 692   return true;
 693 }
 694
 695 //=============================================================================
 696 /*!
 697  *
 698  */
 699 //=============================================================================
 700
 701 FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh &         aMesh,
 702                                                        const TopoDS_Shape & aShape)
 703   //throw(SALOME_Exception)
 704 {
 705   const TopoDS_Face & F = TopoDS::Face(aShape);
 706   const bool ignoreMediumNodes = _quadraticMesh;
 707
 708   // verify 1 wire only, with 4 edges
 709   TopoDS_Vertex V;
 710   list< TopoDS_Edge > edges;
 711   list< int > nbEdgesInWire;
 712   int nbWire = SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
 713   if (nbWire != 1) {
 714     error(COMPERR_BAD_SHAPE, TComm("Wrong number of wires: ") << nbWire);
 715     return 0;
 716   }
 717   FaceQuadStruct* quad = new FaceQuadStruct;
 718   quad->uv_grid = 0;
 719   quad->side.reserve(nbEdgesInWire.front());
 720
 721   int nbSides = 0;
 722   list< TopoDS_Edge >::iterator edgeIt = edges.begin();
 723   if ( nbEdgesInWire.front() == 4 ) { // exactly 4 edges
 724     for ( ; edgeIt != edges.end(); ++edgeIt, nbSides++ )
 725       quad->side.push_back( new StdMeshers_FaceSide(F, *edgeIt, &aMesh,
 726                                                     nbSides<TOP_SIDE, ignoreMediumNodes));
 727   }
 728   else if ( nbEdgesInWire.front() > 4 ) { // more than 4 edges - try to unite some
 729     list< TopoDS_Edge > sideEdges;
 730     while ( !edges.empty()) {
 731       sideEdges.clear();
 732       sideEdges.splice( sideEdges.end(), edges, edges.begin()); // edges.front() -> sideEdges.end()
 733       bool sameSide = true;
 734       while ( !edges.empty() && sameSide ) {
 735         sameSide = SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() );
 736         if ( sameSide )
 737           sideEdges.splice( sideEdges.end(), edges, edges.begin());
 738       }
 739       if ( nbSides == 0 ) { // go backward from the first edge
 740         sameSide = true;
 741         while ( !edges.empty() && sameSide ) {
 742           sameSide = SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() );
 743           if ( sameSide )
 744             sideEdges.splice( sideEdges.begin(), edges, --edges.end());
 745         }
 746       }
 747       quad->side.push_back( new StdMeshers_FaceSide(F, sideEdges, &aMesh,
 748                                                     nbSides<TOP_SIDE, ignoreMediumNodes));
 749       ++nbSides;
 750     }
 751     // issue 20222. Try to unite only edges shared by two same faces
 752     if (nbSides < 4) {
 753       // delete found sides
 754       { FaceQuadStruct cleaner( *quad ); }
 755       quad->side.clear();
 756       quad->side.reserve(nbEdgesInWire.front());
 757       nbSides = 0;
 758
 759       SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
 760       while ( !edges.empty()) {
 761         sideEdges.clear();
 762         sideEdges.splice( sideEdges.end(), edges, edges.begin());
 763         bool sameSide = true;
 764         while ( !edges.empty() && sameSide ) {
 765           sameSide =
 766             SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() ) &&
 767             twoEdgesMeatAtVertex( sideEdges.back(), edges.front(), aMesh );
 768           if ( sameSide )
 769             sideEdges.splice( sideEdges.end(), edges, edges.begin());
 770         }
 771         if ( nbSides == 0 ) { // go backward from the first edge
 772           sameSide = true;
 773           while ( !edges.empty() && sameSide ) {
 774             sameSide =
 775               SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() ) &&
 776               twoEdgesMeatAtVertex( sideEdges.front(), edges.back(), aMesh );
 777             if ( sameSide )
 778               sideEdges.splice( sideEdges.begin(), edges, --edges.end());
 779           }
 780         }
 781         quad->side.push_back( new StdMeshers_FaceSide(F, sideEdges, &aMesh,
 782                                                       nbSides<TOP_SIDE, ignoreMediumNodes));
 783         ++nbSides;
 784       }
 785     }
 786   }
 787   if (nbSides != 4) {
 788 #ifdef _DEBUG_
 789     MESSAGE ( "StdMeshers_Quadrangle_2D. Edge IDs of " << nbSides << " sides:\n" );
 790     for ( int i = 0; i < nbSides; ++i ) {
 791       MESSAGE ( " ( " );
 792       for ( int e = 0; e < quad->side[i]->NbEdges(); ++e )
 793         MESSAGE ( myTool->GetMeshDS()->ShapeToIndex( quad->side[i]->Edge( e )) << " " );
 794       MESSAGE ( ")\n" );
 795     }
 796     //cout << endl;
 797 #endif
 798     if ( !nbSides )
 799       nbSides = nbEdgesInWire.front();
 800     error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides);
 801     delete quad;
 802     quad = 0;
 803   }
 804
 805   return quad;
 806 }
 807
 808
 809 //=============================================================================
 810 /*!
 811  *
 812  */
 813 //=============================================================================
 814
 815 bool StdMeshers_Quadrangle_2D::CheckNbEdgesForEvaluate(SMESH_Mesh& aMesh,
 816                                                        const TopoDS_Shape & aShape,
 817                                                        MapShapeNbElems& aResMap,
 818                                                        std::vector<int>& aNbNodes,
 819                                                        bool& IsQuadratic)
 820
 821 {
 822   const TopoDS_Face & F = TopoDS::Face(aShape);
 823
 824   // verify 1 wire only, with 4 edges
 825   TopoDS_Vertex V;
 826   list< TopoDS_Edge > edges;
 827   list< int > nbEdgesInWire;
 828   int nbWire = SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
 829   if (nbWire != 1) {
 830     return false;
 831   }
 832
 833   aNbNodes.resize(4);
 834
 835   int nbSides = 0;
 836   list< TopoDS_Edge >::iterator edgeIt = edges.begin();
 837   SMESH_subMesh * sm = aMesh.GetSubMesh( *edgeIt );
 838   MapShapeNbElemsItr anIt = aResMap.find(sm);
 839   if(anIt==aResMap.end()) {
 840     return false;
 841   }
 842   std::vector<int> aVec = (*anIt).second;
 843   IsQuadratic = (aVec[2] > aVec[1]);
 844   if ( nbEdgesInWire.front() == 4 ) { // exactly 4 edges
 845     for(; edgeIt != edges.end(); edgeIt++) {
 846       SMESH_subMesh * sm = aMesh.GetSubMesh( *edgeIt );
 847       MapShapeNbElemsItr anIt = aResMap.find(sm);
 848       if(anIt==aResMap.end()) {
 849         return false;
 850       }
 851       std::vector<int> aVec = (*anIt).second;
 852       if(IsQuadratic)
 853         aNbNodes[nbSides] = (aVec[0]-1)/2 + 2;
 854       else
 855         aNbNodes[nbSides] = aVec[0] + 2;
 856       nbSides++;
 857     }
 858   }
 859   else if ( nbEdgesInWire.front() > 4 ) { // more than 4 edges - try to unite some
 860     list< TopoDS_Edge > sideEdges;
 861     while ( !edges.empty()) {
 862       sideEdges.clear();
 863       sideEdges.splice( sideEdges.end(), edges, edges.begin()); // edges.front() -> sideEdges.end()
 864       bool sameSide = true;
 865       while ( !edges.empty() && sameSide ) {
 866         sameSide = SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() );
 867         if ( sameSide )
 868           sideEdges.splice( sideEdges.end(), edges, edges.begin());
 869       }
 870       if ( nbSides == 0 ) { // go backward from the first edge
 871         sameSide = true;
 872         while ( !edges.empty() && sameSide ) {
 873           sameSide = SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() );
 874           if ( sameSide )
 875             sideEdges.splice( sideEdges.begin(), edges, --edges.end());
 876         }
 877       }
 878       list<TopoDS_Edge>::iterator ite = sideEdges.begin();
 879       aNbNodes[nbSides] = 1;
 880       for(; ite!=sideEdges.end(); ite++) {
 881         SMESH_subMesh * sm = aMesh.GetSubMesh( *ite );
 882         MapShapeNbElemsItr anIt = aResMap.find(sm);
 883         if(anIt==aResMap.end()) {
 884           return false;
 885         }
 886         std::vector<int> aVec = (*anIt).second;
 887         if(IsQuadratic)
 888           aNbNodes[nbSides] += (aVec[0]-1)/2 + 1;
 889         else
 890           aNbNodes[nbSides] += aVec[0] + 1;
 891       }
 892       ++nbSides;
 893     }
 894     // issue 20222. Try to unite only edges shared by two same faces
 895     if (nbSides < 4) {
 896       nbSides = 0;
 897       SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
 898       while ( !edges.empty()) {
 899         sideEdges.clear();
 900         sideEdges.splice( sideEdges.end(), edges, edges.begin());
 901         bool sameSide = true;
 902         while ( !edges.empty() && sameSide ) {
 903           sameSide =
 904             SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() ) &&
 905             twoEdgesMeatAtVertex( sideEdges.back(), edges.front(), aMesh );
 906           if ( sameSide )
 907             sideEdges.splice( sideEdges.end(), edges, edges.begin());
 908         }
 909         if ( nbSides == 0 ) { // go backward from the first edge
 910           sameSide = true;
 911           while ( !edges.empty() && sameSide ) {
 912             sameSide =
 913               SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() ) &&
 914               twoEdgesMeatAtVertex( sideEdges.front(), edges.back(), aMesh );
 915             if ( sameSide )
 916               sideEdges.splice( sideEdges.begin(), edges, --edges.end());
 917           }
 918         }
 919         list<TopoDS_Edge>::iterator ite = sideEdges.begin();
 920         aNbNodes[nbSides] = 1;
 921         for(; ite!=sideEdges.end(); ite++) {
 922           SMESH_subMesh * sm = aMesh.GetSubMesh( *ite );
 923           MapShapeNbElemsItr anIt = aResMap.find(sm);
 924           if(anIt==aResMap.end()) {
 925             return false;
 926           }
 927           std::vector<int> aVec = (*anIt).second;
 928           if(IsQuadratic)
 929             aNbNodes[nbSides] += (aVec[0]-1)/2 + 1;
 930           else
 931             aNbNodes[nbSides] += aVec[0] + 1;
 932         }
 933         ++nbSides;
 934       }
 935     }
 936   }
 937   if (nbSides != 4) {
 938     if ( !nbSides )
 939       nbSides = nbEdgesInWire.front();
 940     error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides);
 941     return false;
 942   }
 943
 944   return true;
 945 }
 946
 947
 948 //=============================================================================
 949 /*!
 950  *  CheckAnd2Dcompute
 951  */
 952 //=============================================================================
 953
 954 FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute
 955                            (SMESH_Mesh &         aMesh,
 956                             const TopoDS_Shape & aShape,
 957                             const bool           CreateQuadratic) //throw(SALOME_Exception)
 958 {
 959   _quadraticMesh = CreateQuadratic;
 960
 961   FaceQuadStruct *quad = CheckNbEdges(aMesh, aShape);
 962
 963   if(!quad) return 0;
 964
 965   // set normalized grid on unit square in parametric domain
 966   bool stat = SetNormalizedGrid(aMesh, aShape, quad);
 967   if(!stat) {
 968     if(!quad)
 969       delete quad;
 970     quad = 0;
 971   }
 972
 973   return quad;
 974 }
 975
 976 //=============================================================================
 977 /*!
 978  *
 979  */
 980 //=============================================================================
 981
 982 faceQuadStruct::~faceQuadStruct()
 983 {
 984   for (int i = 0; i < side.size(); i++) {
 985     if (side[i])     delete side[i];
 986   }
 987   if (uv_grid)       delete [] uv_grid;
 988 }
 989
 990 namespace {
 991   inline const vector<UVPtStruct>& GetUVPtStructIn(FaceQuadStruct* quad, int i, int nbSeg)
 992   {
 993     bool   isXConst   = ( i == BOTTOM_SIDE || i == TOP_SIDE );
 994     double constValue = ( i == BOTTOM_SIDE || i == LEFT_SIDE ) ? 0 : 1;
 995     return
 996       quad->isEdgeOut[i] ?
 997       quad->side[i]->SimulateUVPtStruct(nbSeg,isXConst,constValue) :
 998       quad->side[i]->GetUVPtStruct(isXConst,constValue);
 999   }
1000 }
1001
1002 //=============================================================================
1003 /*!
1004  *
1005  */
1006 //=============================================================================
1007
1008 bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh,
1009                                                   const TopoDS_Shape& aShape,
1010                                                   FaceQuadStruct* & quad) //throw (SALOME_Exception)
1011 {
1012   // Algorithme décrit dans "Génération automatique de maillages"
1013   // P.L. GEORGE, MASSON, § 6.4.1 p. 84-85
1014   // traitement dans le domaine paramétrique 2d u,v
1015   // transport - projection sur le carré unité
1016
1017 //  MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid");
1018 //  const TopoDS_Face& F = TopoDS::Face(aShape);
1019
1020   // 1 --- find orientation of the 4 edges, by test on extrema
1021
1022   //      max             min                    0     x1     1
1023   //     |<----north-2-------^                a3 -------------> a2
1024   //     |                   |                   ^1          1^
1025   //    west-3            east-1 =right          |            |
1026   //     |                   |         ==>       |            |
1027   //  y0 |                   | y1                |            |
1028   //     |                   |                   |0          0|
1029   //     v----south-0-------->                a0 -------------> a1
1030   //      min             max                    0     x0     1
1031   //             =down
1032   //
1033
1034   // 3 --- 2D normalized values on unit square [0..1][0..1]
1035
1036   int nbhoriz  = Min(quad->side[0]->NbPoints(), quad->side[2]->NbPoints());
1037   int nbvertic = Min(quad->side[1]->NbPoints(), quad->side[3]->NbPoints());
1038
1039   quad->isEdgeOut[0] = (quad->side[0]->NbPoints() > quad->side[2]->NbPoints());
1040   quad->isEdgeOut[1] = (quad->side[1]->NbPoints() > quad->side[3]->NbPoints());
1041   quad->isEdgeOut[2] = (quad->side[2]->NbPoints() > quad->side[0]->NbPoints());
1042   quad->isEdgeOut[3] = (quad->side[3]->NbPoints() > quad->side[1]->NbPoints());
1043
1044   UVPtStruct *uv_grid = quad->uv_grid = new UVPtStruct[nbvertic * nbhoriz];
1045
1046   const vector<UVPtStruct>& uv_e0 = GetUVPtStructIn( quad, 0, nbhoriz - 1 );
1047   const vector<UVPtStruct>& uv_e1 = GetUVPtStructIn( quad, 1, nbvertic - 1 );
1048   const vector<UVPtStruct>& uv_e2 = GetUVPtStructIn( quad, 2, nbhoriz - 1 );
1049   const vector<UVPtStruct>& uv_e3 = GetUVPtStructIn( quad, 3, nbvertic - 1 );
1050
1051   if ( uv_e0.empty() || uv_e1.empty() || uv_e2.empty() || uv_e3.empty() )
1052     //return error( "Can't find nodes on sides");
1053     return error( COMPERR_BAD_INPUT_MESH );
1054
1055   // nodes Id on "in" edges
1056   if (! quad->isEdgeOut[0]) {
1057     int j = 0;
1058     for (int i = 0; i < nbhoriz; i++) { // down
1059       int ij = j * nbhoriz + i;
1060       uv_grid[ij].node = uv_e0[i].node;
1061     }
1062   }
1063   if (! quad->isEdgeOut[1]) {
1064     int i = nbhoriz - 1;
1065     for (int j = 0; j < nbvertic; j++) { // right
1066       int ij = j * nbhoriz + i;
1067       uv_grid[ij].node = uv_e1[j].node;
1068     }
1069   }
1070   if (! quad->isEdgeOut[2]) {
1071     int j = nbvertic - 1;
1072     for (int i = 0; i < nbhoriz; i++) { // up
1073       int ij = j * nbhoriz + i;
1074       uv_grid[ij].node = uv_e2[i].node;
1075     }
1076   }
1077   if (! quad->isEdgeOut[3]) {
1078     int i = 0;
1079     for (int j = 0; j < nbvertic; j++) { // left
1080       int ij = j * nbhoriz + i;
1081       uv_grid[ij].node = uv_e3[j].node;
1082     }
1083   }
1084
1085   // normalized 2d values on grid
1086   for (int i = 0; i < nbhoriz; i++)
1087   {
1088     for (int j = 0; j < nbvertic; j++)
1089     {
1090       int ij = j * nbhoriz + i;
1091       // --- droite i cste : x = x0 + y(x1-x0)
1092       double x0 = uv_e0[i].normParam;   // bas - sud
1093       double x1 = uv_e2[i].normParam;   // haut - nord
1094       // --- droite j cste : y = y0 + x(y1-y0)
1095       double y0 = uv_e3[j].normParam;   // gauche-ouest
1096       double y1 = uv_e1[j].normParam;   // droite - est
1097       // --- intersection : x=x0+(y0+x(y1-y0))(x1-x0)
1098       double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
1099       double y = y0 + x * (y1 - y0);
1100       uv_grid[ij].x = x;
1101       uv_grid[ij].y = y;
1102       //MESSAGE("-xy-01 "<<x0<<" "<<x1<<" "<<y0<<" "<<y1);
1103       //MESSAGE("-xy-norm "<<i<<" "<<j<<" "<<x<<" "<<y);
1104     }
1105   }
1106
1107   // 4 --- projection on 2d domain (u,v)
1108   gp_UV a0( uv_e0.front().u, uv_e0.front().v );
1109   gp_UV a1( uv_e0.back().u,  uv_e0.back().v );
1110   gp_UV a2( uv_e2.back().u,  uv_e2.back().v );
1111   gp_UV a3( uv_e2.front().u, uv_e2.front().v );
1112
1113   for (int i = 0; i < nbhoriz; i++)
1114   {
1115     for (int j = 0; j < nbvertic; j++)
1116     {
1117       int ij = j * nbhoriz + i;
1118       double x = uv_grid[ij].x;
1119       double y = uv_grid[ij].y;
1120       double param_0 = uv_e0[0].normParam + x * (uv_e0.back().normParam - uv_e0[0].normParam); // sud
1121       double param_2 = uv_e2[0].normParam + x * (uv_e2.back().normParam - uv_e2[0].normParam); // nord
1122       double param_1 = uv_e1[0].normParam + y * (uv_e1.back().normParam - uv_e1[0].normParam); // est
1123       double param_3 = uv_e3[0].normParam + y * (uv_e3.back().normParam - uv_e3[0].normParam); // ouest
1124
1125       //MESSAGE("params "<<param_0<<" "<<param_1<<" "<<param_2<<" "<<param_3);
1126       gp_UV p0 = quad->side[0]->Value2d(param_0).XY();
1127       gp_UV p1 = quad->side[1]->Value2d(param_1).XY();
1128       gp_UV p2 = quad->side[2]->Value2d(param_2).XY();
1129       gp_UV p3 = quad->side[3]->Value2d(param_3).XY();
1130
1131       gp_UV uv = (1 - y) * p0 + x * p1 + y * p2 + (1 - x) * p3;
1132       uv -= (1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3;
1133
1134       uv_grid[ij].u = uv.X();
1135       uv_grid[ij].v = uv.Y();
1136     }
1137   }
1138   return true;
1139 }
1140
1141 //=======================================================================
1142 //function : ShiftQuad
1143 //purpose  : auxilary function for ComputeQuadPref
1144 //=======================================================================
1145
1146 static void ShiftQuad(FaceQuadStruct* quad, const int num, bool)
1147 {
1148   StdMeshers_FaceSide* side[4] = { quad->side[0], quad->side[1], quad->side[2], quad->side[3] };
1149   for (int i = BOTTOM_SIDE; i < NB_SIDES; ++i ) {
1150     int id = ( i + num ) % NB_SIDES;
1151     bool wasForward = ( i < TOP_SIDE );
1152     bool newForward = ( id < TOP_SIDE );
1153     if ( wasForward != newForward )
1154       side[ i ]->Reverse();
1155     quad->side[ id ] = side[ i ];
1156   }
1157 }
1158
1159 //=======================================================================
1160 //function : CalcUV
1161 //purpose  : auxilary function for ComputeQuadPref
1162 //=======================================================================
1163
1164 static gp_UV CalcUV(double x0, double x1, double y0, double y1,
1165                     FaceQuadStruct* quad,
1166                     const gp_UV& a0, const gp_UV& a1,
1167                     const gp_UV& a2, const gp_UV& a3)
1168 {
1169   const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0 );
1170   const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
1171   const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1 );
1172   const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
1173
1174   double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
1175   double y = y0 + x * (y1 - y0);
1176
1177   double param_b = uv_eb[0].normParam + x * (uv_eb.back().normParam - uv_eb[0].normParam);
1178   double param_t = uv_et[0].normParam + x * (uv_et.back().normParam - uv_et[0].normParam);
1179   double param_r = uv_er[0].normParam + y * (uv_er.back().normParam - uv_er[0].normParam);
1180   double param_l = uv_el[0].normParam + y * (uv_el.back().normParam - uv_el[0].normParam);
1181
1182   gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(param_b).XY();
1183   gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(param_r).XY();
1184   gp_UV p2 = quad->side[TOP_SIDE   ]->Value2d(param_t).XY();
1185   gp_UV p3 = quad->side[LEFT_SIDE  ]->Value2d(param_l).XY();
1186
1187   gp_UV uv = p0 * (1 - y) + p1 * x + p2 * y + p3 * (1 - x);
1188
1189   uv -= (1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3;
1190
1191   return uv;
1192 }
1193
1194 //=======================================================================
1195 //function : CalcUV2
1196 //purpose  : auxilary function for ComputeQuadPref
1197 //=======================================================================
1198
1199 static gp_UV CalcUV2(double x, double y,
1200                      FaceQuadStruct* quad,
1201                      const gp_UV& a0, const gp_UV& a1,
1202                      const gp_UV& a2, const gp_UV& a3)
1203 {
1204   const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0 );
1205   const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
1206   const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1 );
1207   const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
1208
1209   //double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
1210   //double y = y0 + x * (y1 - y0);
1211
1212   double param_b = uv_eb[0].normParam + x * (uv_eb.back().normParam - uv_eb[0].normParam);
1213   double param_t = uv_et[0].normParam + x * (uv_et.back().normParam - uv_et[0].normParam);
1214   double param_r = uv_er[0].normParam + y * (uv_er.back().normParam - uv_er[0].normParam);
1215   double param_l = uv_el[0].normParam + y * (uv_el.back().normParam - uv_el[0].normParam);
1216
1217   gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(param_b).XY();
1218   gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(param_r).XY();
1219   gp_UV p2 = quad->side[TOP_SIDE   ]->Value2d(param_t).XY();
1220   gp_UV p3 = quad->side[LEFT_SIDE  ]->Value2d(param_l).XY();
1221
1222   gp_UV uv = p0 * (1 - y) + p1 * x + p2 * y + p3 * (1 - x);
1223
1224   uv -= (1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3;
1225
1226   return uv;
1227 }
1228
1229
1230 //=======================================================================
1231 /*!
1232  * Create only quandrangle faces
1233  */
1234 //=======================================================================
1235
1236 bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh &        aMesh,
1237                                                 const TopoDS_Shape& aShape,
1238                                                 FaceQuadStruct*     quad)
1239 {
1240   // Auxilary key in order to keep old variant
1241   // of meshing after implementation new variant
1242   // for bug 0016220 from Mantis.
1243   bool OldVersion = false;
1244
1245   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
1246   const TopoDS_Face& F = TopoDS::Face(aShape);
1247   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
1248 //  const TopoDS_Wire& W = BRepTools::OuterWire(F);
1249   bool WisF = true;
1250 //   if(W.Orientation()==TopAbs_FORWARD)
1251 //     WisF = true;
1252   //if(WisF) cout<<"W is FORWARD"<<endl;
1253   //else cout<<"W is REVERSED"<<endl;
1254 //   bool FisF = (F.Orientation()==TopAbs_FORWARD);
1255 //   if(!FisF) WisF = !WisF;
1256 //  WisF = FisF;
1257   int i,j,geomFaceID = meshDS->ShapeToIndex( F );
1258
1259   int nb = quad->side[0]->NbPoints();
1260   int nr = quad->side[1]->NbPoints();
1261   int nt = quad->side[2]->NbPoints();
1262   int nl = quad->side[3]->NbPoints();
1263   int dh = abs(nb-nt);
1264   int dv = abs(nr-nl);
1265
1266   if( dh>=dv ) {
1267     if( nt>nb ) {
1268       // it is a base case => not shift quad but me be replacement is need
1269       ShiftQuad(quad,0,WisF);
1270     }
1271     else {
1272       // we have to shift quad on 2
1273       ShiftQuad(quad,2,WisF);
1274     }
1275   }
1276   else {
1277     if( nr>nl ) {
1278       // we have to shift quad on 1
1279       ShiftQuad(quad,1,WisF);
1280     }
1281     else {
1282       // we have to shift quad on 3
1283       ShiftQuad(quad,3,WisF);
1284     }
1285   }
1286
1287   nb = quad->side[0]->NbPoints();
1288   nr = quad->side[1]->NbPoints();
1289   nt = quad->side[2]->NbPoints();
1290   nl = quad->side[3]->NbPoints();
1291   dh = abs(nb-nt);
1292   dv = abs(nr-nl);
1293   int nbh  = Max(nb,nt);
1294   int nbv = Max(nr,nl);
1295   int addh = 0;
1296   int addv = 0;
1297
1298   // ----------- Old version ---------------
1299   // orientation of face and 3 main domain for future faces
1300   //       0   top    1
1301   //      1------------1
1302   //       |   |  |   |
1303   //       |   |  |   |
1304   //       | L |  | R |
1305   //  left |   |  |   | rigth
1306   //       |  /    \  |
1307   //       | /  C   \ |
1308   //       |/        \|
1309   //      0------------0
1310   //       0  bottom  1
1311
1312   // ----------- New version ---------------
1313   // orientation of face and 3 main domain for future faces
1314   //       0   top    1
1315   //      1------------1
1316   //       |  |____|  |
1317   //       |  /    \  |
1318   //       | /  C   \ |
1319   //  left |/________\| rigth
1320   //       |          |
1321   //       |          |
1322   //       |          |
1323   //      0------------0
1324   //       0  bottom  1
1325
1326   if(dh>dv) {
1327     addv = (dh-dv)/2;
1328     nbv = nbv + addv;
1329   }
1330   else { // dv>=dh
1331     addh = (dv-dh)/2;
1332     nbh = nbh + addh;
1333   }
1334
1335   const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0 );
1336   const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
1337   const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1 );
1338   const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
1339
1340   // arrays for normalized params
1341   //cout<<"Dump B:"<<endl;
1342   TColStd_SequenceOfReal npb, npr, npt, npl;
1343   for(i=0; i<nb; i++) {
1344     npb.Append(uv_eb[i].normParam);
1345     //cout<<"i="<<i<<" par="<<uv_eb[i].normParam<<" npar="<<uv_eb[i].normParam;
1346     //const SMDS_MeshNode* N = uv_eb[i].node;
1347     //cout<<" node("<<N->X()<<","<<N->Y()<<","<<N->Z()<<")"<<endl;
1348   }
1349   for(i=0; i<nr; i++) {
1350     npr.Append(uv_er[i].normParam);
1351   }
1352   for(i=0; i<nt; i++) {
1353     npt.Append(uv_et[i].normParam);
1354   }
1355   for(i=0; i<nl; i++) {
1356     npl.Append(uv_el[i].normParam);
1357   }
1358
1359   int dl,dr;
1360   if(OldVersion) {
1361     // add some params to right and left after the first param
1362     // insert to right
1363     dr = nbv - nr;
1364     double dpr = (npr.Value(2) - npr.Value(1))/(dr+1);
1365     for(i=1; i<=dr; i++) {
1366       npr.InsertAfter(1,npr.Value(2)-dpr);
1367     }
1368     // insert to left
1369     dl = nbv - nl;
1370     dpr = (npl.Value(2) - npl.Value(1))/(dl+1);
1371     for(i=1; i<=dl; i++) {
1372       npl.InsertAfter(1,npl.Value(2)-dpr);
1373     }
1374   }
1375   //cout<<"npb:";
1376   //for(i=1; i<=npb.Length(); i++) {
1377   //  cout<<" "<<npb.Value(i);
1378   //}
1379   //cout<<endl;
1380
1381   gp_XY a0( uv_eb.front().u, uv_eb.front().v );
1382   gp_XY a1( uv_eb.back().u,  uv_eb.back().v );
1383   gp_XY a2( uv_et.back().u,  uv_et.back().v );
1384   gp_XY a3( uv_et.front().u, uv_et.front().v );
1385   //cout<<" a0("<<a0.X()<<","<<a0.Y()<<")"<<" a1("<<a1.X()<<","<<a1.Y()<<")"
1386   //    <<" a2("<<a2.X()<<","<<a2.Y()<<")"<<" a3("<<a3.X()<<","<<a3.Y()<<")"<<endl;
1387
1388   int nnn = Min(nr,nl);
1389   // auxilary sequence of XY for creation nodes
1390   // in the bottom part of central domain
1391   // it's length must be == nbv-nnn-1
1392   TColgp_SequenceOfXY UVL;
1393   TColgp_SequenceOfXY UVR;
1394
1395   if(OldVersion) {
1396     // step1: create faces for left domain
1397     StdMeshers_Array2OfNode NodesL(1,dl+1,1,nl);
1398     // add left nodes
1399     for(j=1; j<=nl; j++)
1400       NodesL.SetValue(1,j,uv_el[j-1].node);
1401     if(dl>0) {
1402       // add top nodes
1403       for(i=1; i<=dl; i++)
1404         NodesL.SetValue(i+1,nl,uv_et[i].node);
1405       // create and add needed nodes
1406       TColgp_SequenceOfXY UVtmp;
1407       for(i=1; i<=dl; i++) {
1408         double x0 = npt.Value(i+1);
1409         double x1 = x0;
1410         // diagonal node
1411         double y0 = npl.Value(i+1);
1412         double y1 = npr.Value(i+1);
1413         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1414         gp_Pnt P = S->Value(UV.X(),UV.Y());
1415         SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1416         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1417         NodesL.SetValue(i+1,1,N);
1418         if(UVL.Length()<nbv-nnn-1) UVL.Append(UV);
1419         // internal nodes
1420         for(j=2; j<nl; j++) {
1421           double y0 = npl.Value(dl+j);
1422           double y1 = npr.Value(dl+j);
1423           gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1424           gp_Pnt P = S->Value(UV.X(),UV.Y());
1425           SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1426           meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1427           NodesL.SetValue(i+1,j,N);
1428           if( i==dl ) UVtmp.Append(UV);
1429         }
1430       }
1431       for(i=1; i<=UVtmp.Length() && UVL.Length()<nbv-nnn-1; i++) {
1432         UVL.Append(UVtmp.Value(i));
1433       }
1434       //cout<<"Dump NodesL:"<<endl;
1435       //for(i=1; i<=dl+1; i++) {
1436       //  cout<<"i="<<i;
1437       //  for(j=1; j<=nl; j++) {
1438       //    cout<<" ("<<NodesL.Value(i,j)->X()<<","<<NodesL.Value(i,j)->Y()<<","<<NodesL.Value(i,j)->Z()<<")";
1439       //  }
1440       //  cout<<endl;
1441       //}
1442       // create faces
1443       for(i=1; i<=dl; i++) {
1444         for(j=1; j<nl; j++) {
1445           if(WisF) {
1446             SMDS_MeshFace* F =
1447               myTool->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j),
1448                               NodesL.Value(i+1,j+1), NodesL.Value(i,j+1));
1449             meshDS->SetMeshElementOnShape(F, geomFaceID);
1450           }
1451           else {
1452             SMDS_MeshFace* F =
1453               myTool->AddFace(NodesL.Value(i,j), NodesL.Value(i,j+1),
1454                               NodesL.Value(i+1,j+1), NodesL.Value(i+1,j));
1455             meshDS->SetMeshElementOnShape(F, geomFaceID);
1456           }
1457         }
1458       }
1459     }
1460     else {
1461       // fill UVL using c2d
1462       for(i=1; i<npl.Length() && UVL.Length()<nbv-nnn-1; i++) {
1463         UVL.Append( gp_UV ( uv_el[i].u, uv_el[i].v ));
1464       }
1465     }
1466
1467     // step2: create faces for right domain
1468     StdMeshers_Array2OfNode NodesR(1,dr+1,1,nr);
1469     // add right nodes
1470     for(j=1; j<=nr; j++)
1471       NodesR.SetValue(1,j,uv_er[nr-j].node);
1472     if(dr>0) {
1473       // add top nodes
1474       for(i=1; i<=dr; i++)
1475         NodesR.SetValue(i+1,1,uv_et[nt-1-i].node);
1476       // create and add needed nodes
1477       TColgp_SequenceOfXY UVtmp;
1478       for(i=1; i<=dr; i++) {
1479         double x0 = npt.Value(nt-i);
1480         double x1 = x0;
1481         // diagonal node
1482         double y0 = npl.Value(i+1);
1483         double y1 = npr.Value(i+1);
1484         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1485         gp_Pnt P = S->Value(UV.X(),UV.Y());
1486         SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1487         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1488         NodesR.SetValue(i+1,nr,N);
1489         if(UVR.Length()<nbv-nnn-1) UVR.Append(UV);
1490         // internal nodes
1491         for(j=2; j<nr; j++) {
1492           double y0 = npl.Value(nbv-j+1);
1493           double y1 = npr.Value(nbv-j+1);
1494           gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1495           gp_Pnt P = S->Value(UV.X(),UV.Y());
1496           SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1497           meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1498           NodesR.SetValue(i+1,j,N);
1499           if( i==dr ) UVtmp.Prepend(UV);
1500         }
1501       }
1502       for(i=1; i<=UVtmp.Length() && UVR.Length()<nbv-nnn-1; i++) {
1503         UVR.Append(UVtmp.Value(i));
1504       }
1505       // create faces
1506       for(i=1; i<=dr; i++) {
1507         for(j=1; j<nr; j++) {
1508           if(WisF) {
1509             SMDS_MeshFace* F =
1510               myTool->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j),
1511                               NodesR.Value(i+1,j+1), NodesR.Value(i,j+1));
1512             meshDS->SetMeshElementOnShape(F, geomFaceID);
1513           }
1514           else {
1515             SMDS_MeshFace* F =
1516               myTool->AddFace(NodesR.Value(i,j), NodesR.Value(i,j+1),
1517                               NodesR.Value(i+1,j+1), NodesR.Value(i+1,j));
1518             meshDS->SetMeshElementOnShape(F, geomFaceID);
1519           }
1520         }
1521       }
1522     }
1523     else {
1524       // fill UVR using c2d
1525       for(i=1; i<npr.Length() && UVR.Length()<nbv-nnn-1; i++) {
1526         UVR.Append( gp_UV( uv_er[i].u, uv_er[i].v ));
1527       }
1528     }
1529
1530     // step3: create faces for central domain
1531     StdMeshers_Array2OfNode NodesC(1,nb,1,nbv);
1532     // add first string using NodesL
1533     for(i=1; i<=dl+1; i++)
1534       NodesC.SetValue(1,i,NodesL(i,1));
1535     for(i=2; i<=nl; i++)
1536       NodesC.SetValue(1,dl+i,NodesL(dl+1,i));
1537     // add last string using NodesR
1538     for(i=1; i<=dr+1; i++)
1539       NodesC.SetValue(nb,i,NodesR(i,nr));
1540     for(i=1; i<nr; i++)
1541       NodesC.SetValue(nb,dr+i+1,NodesR(dr+1,nr-i));
1542     // add top nodes (last columns)
1543     for(i=dl+2; i<nbh-dr; i++)
1544       NodesC.SetValue(i-dl,nbv,uv_et[i-1].node);
1545     // add bottom nodes (first columns)
1546     for(i=2; i<nb; i++)
1547       NodesC.SetValue(i,1,uv_eb[i-1].node);
1548
1549     // create and add needed nodes
1550     // add linear layers
1551     for(i=2; i<nb; i++) {
1552       double x0 = npt.Value(dl+i);
1553       double x1 = x0;
1554       for(j=1; j<nnn; j++) {
1555         double y0 = npl.Value(nbv-nnn+j);
1556         double y1 = npr.Value(nbv-nnn+j);
1557         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1558         gp_Pnt P = S->Value(UV.X(),UV.Y());
1559         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1560         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1561         NodesC.SetValue(i,nbv-nnn+j,N);
1562       }
1563     }
1564     // add diagonal layers
1565     //cout<<"UVL.Length()="<<UVL.Length()<<" UVR.Length()="<<UVR.Length()<<endl;
1566     //cout<<"Dump UVL:"<<endl;
1567     //for(i=1; i<=UVL.Length(); i++) {
1568     //  cout<<" ("<<UVL.Value(i).X()<<","<<UVL.Value(i).Y()<<")";
1569     //}
1570     //cout<<endl;
1571     for(i=1; i<nbv-nnn; i++) {
1572       double du = UVR.Value(i).X() - UVL.Value(i).X();
1573       double dv = UVR.Value(i).Y() - UVL.Value(i).Y();
1574       for(j=2; j<nb; j++) {
1575         double u = UVL.Value(i).X() + du*npb.Value(j);
1576         double v = UVL.Value(i).Y() + dv*npb.Value(j);
1577         gp_Pnt P = S->Value(u,v);
1578         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1579         meshDS->SetNodeOnFace(N, geomFaceID, u, v);
1580         NodesC.SetValue(j,i+1,N);
1581       }
1582     }
1583     // create faces
1584     for(i=1; i<nb; i++) {
1585       for(j=1; j<nbv; j++) {
1586         if(WisF) {
1587           SMDS_MeshFace* F =
1588             myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
1589                             NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
1590           meshDS->SetMeshElementOnShape(F, geomFaceID);
1591         }
1592         else {
1593           SMDS_MeshFace* F =
1594             myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1),
1595                             NodesC.Value(i+1,j+1), NodesC.Value(i+1,j));
1596           meshDS->SetMeshElementOnShape(F, geomFaceID);
1597         }
1598       }
1599     }
1600   }
1601
1602   else { // New version (!OldVersion)
1603     // step1: create faces for bottom rectangle domain
1604     StdMeshers_Array2OfNode NodesBRD(1,nb,1,nnn-1);
1605     // fill UVL and UVR using c2d
1606     for(j=0; j<nb; j++) {
1607       NodesBRD.SetValue(j+1,1,uv_eb[j].node);
1608     }
1609     for(i=1; i<nnn-1; i++) {
1610       NodesBRD.SetValue(1,i+1,uv_el[i].node);
1611       NodesBRD.SetValue(nb,i+1,uv_er[i].node);
1612       double du = uv_er[i].u - uv_el[i].u;
1613       double dv = uv_er[i].v - uv_el[i].v;
1614       for(j=2; j<nb; j++) {
1615         double u = uv_el[i].u + du*npb.Value(j);
1616         double v = uv_el[i].v + dv*npb.Value(j);
1617         gp_Pnt P = S->Value(u,v);
1618         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1619         meshDS->SetNodeOnFace(N, geomFaceID, u, v);
1620         NodesBRD.SetValue(j,i+1,N);
1621
1622       }
1623     }
1624     int nbf=0;
1625     for(j=1; j<nnn-1; j++) {
1626       for(i=1; i<nb; i++) {
1627         nbf++;
1628         if(WisF) {
1629           SMDS_MeshFace* F =
1630             myTool->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i+1,j),
1631                             NodesBRD.Value(i+1,j+1), NodesBRD.Value(i,j+1));
1632           meshDS->SetMeshElementOnShape(F, geomFaceID);
1633         }
1634         else {
1635           SMDS_MeshFace* F =
1636             myTool->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i,j+1),
1637                             NodesBRD.Value(i+1,j+1), NodesBRD.Value(i+1,j));
1638           meshDS->SetMeshElementOnShape(F, geomFaceID);
1639         }
1640       }
1641     }
1642     int drl = abs(nr-nl);
1643     // create faces for region C
1644     StdMeshers_Array2OfNode NodesC(1,nb,1,drl+1+addv);
1645     // add nodes from previous region
1646     for(j=1; j<=nb; j++) {
1647       NodesC.SetValue(j,1,NodesBRD.Value(j,nnn-1));
1648     }
1649     if( (drl+addv) > 0 ) {
1650       int n1,n2;
1651       if(nr>nl) {
1652         n1 = 1;
1653         n2 = drl + 1;
1654         TColgp_SequenceOfXY UVtmp;
1655         double drparam = npr.Value(nr) - npr.Value(nnn-1);
1656         double dlparam = npl.Value(nnn) - npl.Value(nnn-1);
1657         double y0,y1;
1658         for(i=1; i<=drl; i++) {
1659           // add existed nodes from right edge
1660           NodesC.SetValue(nb,i+1,uv_er[nnn+i-2].node);
1661           //double dtparam = npt.Value(i+1);
1662           y1 = npr.Value(nnn+i-1); // param on right edge
1663           double dpar = (y1 - npr.Value(nnn-1))/drparam;
1664           y0 = npl.Value(nnn-1) + dpar*dlparam; // param on left edge
1665           double dy = y1 - y0;
1666           for(j=1; j<nb; j++) {
1667             double x = npt.Value(i+1) + npb.Value(j)*(1-npt.Value(i+1));
1668             double y = y0 + dy*x;
1669             gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1670             gp_Pnt P = S->Value(UV.X(),UV.Y());
1671             SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1672             meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1673             NodesC.SetValue(j,i+1,N);
1674           }
1675         }
1676         double dy0 = (1-y0)/(addv+1);
1677         double dy1 = (1-y1)/(addv+1);
1678         for(i=1; i<=addv; i++) {
1679           double yy0 = y0 + dy0*i;
1680           double yy1 = y1 + dy1*i;
1681           double dyy = yy1 - yy0;
1682           for(j=1; j<=nb; j++) {
1683             double x = npt.Value(i+1+drl) +
1684               npb.Value(j) * ( npt.Value(nt-i) - npt.Value(i+1+drl) );
1685             double y = yy0 + dyy*x;
1686             gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1687             gp_Pnt P = S->Value(UV.X(),UV.Y());
1688             SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1689             meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1690             NodesC.SetValue(j,i+drl+1,N);
1691           }
1692         }
1693       }
1694       else { // nr<nl
1695         n2 = 1;
1696         n1 = drl + 1;
1697         TColgp_SequenceOfXY UVtmp;
1698         double dlparam = npl.Value(nl) - npl.Value(nnn-1);
1699         double drparam = npr.Value(nnn) - npr.Value(nnn-1);
1700         double y0 = npl.Value(nnn-1);
1701         double y1 = npr.Value(nnn-1);
1702         for(i=1; i<=drl; i++) {
1703           // add existed nodes from right edge
1704           NodesC.SetValue(1,i+1,uv_el[nnn+i-2].node);
1705           y0 = npl.Value(nnn+i-1); // param on left edge
1706           double dpar = (y0 - npl.Value(nnn-1))/dlparam;
1707           y1 = npr.Value(nnn-1) + dpar*drparam; // param on right edge
1708           double dy = y1 - y0;
1709           for(j=2; j<=nb; j++) {
1710             double x = npb.Value(j)*npt.Value(nt-i);
1711             double y = y0 + dy*x;
1712             gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1713             gp_Pnt P = S->Value(UV.X(),UV.Y());
1714             SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1715             meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1716             NodesC.SetValue(j,i+1,N);
1717           }
1718         }
1719         double dy0 = (1-y0)/(addv+1);
1720         double dy1 = (1-y1)/(addv+1);
1721         for(i=1; i<=addv; i++) {
1722           double yy0 = y0 + dy0*i;
1723           double yy1 = y1 + dy1*i;
1724           double dyy = yy1 - yy0;
1725           for(j=1; j<=nb; j++) {
1726             double x = npt.Value(i+1) +
1727               npb.Value(j) * ( npt.Value(nt-i-drl) - npt.Value(i+1) );
1728             double y = yy0 + dyy*x;
1729             gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1730             gp_Pnt P = S->Value(UV.X(),UV.Y());
1731             SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1732             meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1733             NodesC.SetValue(j,i+drl+1,N);
1734           }
1735         }
1736       }
1737       // create faces
1738       for(j=1; j<=drl+addv; j++) {
1739         for(i=1; i<nb; i++) {
1740           nbf++;
1741           if(WisF) {
1742             SMDS_MeshFace* F =
1743               myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
1744                               NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
1745             meshDS->SetMeshElementOnShape(F, geomFaceID);
1746           }
1747           else {
1748             SMDS_MeshFace* F =
1749               myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1),
1750                               NodesC.Value(i+1,j+1), NodesC.Value(i+1,j));
1751             meshDS->SetMeshElementOnShape(F, geomFaceID);
1752           }
1753         }
1754       } // end nr<nl
1755
1756       StdMeshers_Array2OfNode NodesLast(1,nt,1,2);
1757       for(i=1; i<=nt; i++) {
1758         NodesLast.SetValue(i,2,uv_et[i-1].node);
1759       }
1760       int nnn=0;
1761       for(i=n1; i<drl+addv+1; i++) {
1762         nnn++;
1763         NodesLast.SetValue(nnn,1,NodesC.Value(1,i));
1764       }
1765       for(i=1; i<=nb; i++) {
1766         nnn++;
1767         NodesLast.SetValue(nnn,1,NodesC.Value(i,drl+addv+1));
1768       }
1769       for(i=drl+addv; i>=n2; i--) {
1770         nnn++;
1771         NodesLast.SetValue(nnn,1,NodesC.Value(nb,i));
1772       }
1773       for(i=1; i<nt; i++) {
1774         nbf++;
1775         if(WisF) {
1776           SMDS_MeshFace* F =
1777             myTool->AddFace(NodesLast.Value(i,1), NodesLast.Value(i+1,1),
1778                             NodesLast.Value(i+1,2), NodesLast.Value(i,2));
1779           meshDS->SetMeshElementOnShape(F, geomFaceID);
1780         }
1781         else {
1782           SMDS_MeshFace* F =
1783             myTool->AddFace(NodesLast.Value(i,1), NodesLast.Value(i,2),
1784                             NodesLast.Value(i+1,2), NodesLast.Value(i+1,2));
1785           meshDS->SetMeshElementOnShape(F, geomFaceID);
1786         }
1787       }
1788     } // if( (drl+addv) > 0 )
1789
1790   } // end new version implementation
1791
1792   bool isOk = true;
1793   return isOk;
1794 }
1795
1796
1797 //=======================================================================
1798 /*!
1799  * Evaluate only quandrangle faces
1800  */
1801 //=======================================================================
1802
1803 bool StdMeshers_Quadrangle_2D::EvaluateQuadPref(SMESH_Mesh &        aMesh,
1804                                                 const TopoDS_Shape& aShape,
1805                                                 std::vector<int>& aNbNodes,
1806                                                 MapShapeNbElems& aResMap,
1807                                                 bool IsQuadratic)
1808 {
1809   // Auxilary key in order to keep old variant
1810   // of meshing after implementation new variant
1811   // for bug 0016220 from Mantis.
1812   bool OldVersion = false;
1813
1814   const TopoDS_Face& F = TopoDS::Face(aShape);
1815   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
1816
1817   int nb = aNbNodes[0];
1818   int nr = aNbNodes[1];
1819   int nt = aNbNodes[2];
1820   int nl = aNbNodes[3];
1821   int dh = abs(nb-nt);
1822   int dv = abs(nr-nl);
1823
1824   if( dh>=dv ) {
1825     if( nt>nb ) {
1826       // it is a base case => not shift
1827     }
1828     else {
1829       // we have to shift on 2
1830       nb = aNbNodes[2];
1831       nr = aNbNodes[3];
1832       nt = aNbNodes[0];
1833       nl = aNbNodes[1];
1834     }
1835   }
1836   else {
1837     if( nr>nl ) {
1838       // we have to shift quad on 1
1839       nb = aNbNodes[3];
1840       nr = aNbNodes[0];
1841       nt = aNbNodes[1];
1842       nl = aNbNodes[2];
1843     }
1844     else {
1845       // we have to shift quad on 3
1846       nb = aNbNodes[1];
1847       nr = aNbNodes[2];
1848       nt = aNbNodes[3];
1849       nl = aNbNodes[0];
1850     }
1851   }
1852
1853   dh = abs(nb-nt);
1854   dv = abs(nr-nl);
1855   int nbh  = Max(nb,nt);
1856   int nbv = Max(nr,nl);
1857   int addh = 0;
1858   int addv = 0;
1859
1860   if(dh>dv) {
1861     addv = (dh-dv)/2;
1862     nbv = nbv + addv;
1863   }
1864   else { // dv>=dh
1865     addh = (dv-dh)/2;
1866     nbh = nbh + addh;
1867   }
1868
1869   int dl,dr;
1870   if(OldVersion) {
1871     // add some params to right and left after the first param
1872     // insert to right
1873     dr = nbv - nr;
1874     // insert to left
1875     dl = nbv - nl;
1876   }
1877
1878   int nnn = Min(nr,nl);
1879
1880   int nbNodes = 0;
1881   int nbFaces = 0;
1882   if(OldVersion) {
1883     // step1: create faces for left domain
1884     if(dl>0) {
1885       nbNodes += dl*(nl-1);
1886       nbFaces += dl*(nl-1);
1887     }
1888     // step2: create faces for right domain
1889     if(dr>0) {
1890       nbNodes += dr*(nr-1);
1891       nbFaces += dr*(nr-1);
1892     }
1893     // step3: create faces for central domain
1894     nbNodes += (nb-2)*(nnn-1) + (nbv-nnn-1)*(nb-2);
1895     nbFaces += (nb-1)*(nbv-1);
1896   }
1897   else { // New version (!OldVersion)
1898     nbNodes += (nnn-2)*(nb-2);
1899     nbFaces += (nnn-2)*(nb-1);
1900     int drl = abs(nr-nl);
1901     nbNodes += drl*(nb-1) + addv*nb;
1902     nbFaces += (drl+addv)*(nb-1) + (nt-1);
1903   } // end new version implementation
1904
1905   std::vector<int> aVec(17);
1906   for(int i=0; i<17; i++) aVec[i] = 0;
1907   if(IsQuadratic) {
1908     aVec[6] = nbFaces;
1909     aVec[0] = nbNodes + nbFaces*4;
1910   }
1911   else {
1912     aVec[0] = nbNodes;
1913     aVec[5] = nbFaces;
1914   }
1915   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
1916   aResMap.insert(std::make_pair(sm,aVec));
1917
1918   return true;
1919 }
1920
1921
1922 //=============================================================================
1923 /*! Split quadrangle in to 2 triangles by smallest diagonal
1924  *
1925  */
1926 //=============================================================================
1927 void StdMeshers_Quadrangle_2D::SplitQuad(SMESHDS_Mesh *theMeshDS,
1928                                     int theFaceID,
1929                                     const SMDS_MeshNode* theNode1,
1930                                     const SMDS_MeshNode* theNode2,
1931                                     const SMDS_MeshNode* theNode3,
1932                                     const SMDS_MeshNode* theNode4)
1933 {
1934   gp_Pnt a(theNode1->X(),theNode1->Y(),theNode1->Z());
1935   gp_Pnt b(theNode2->X(),theNode2->Y(),theNode2->Z());
1936   gp_Pnt c(theNode3->X(),theNode3->Y(),theNode3->Z());
1937   gp_Pnt d(theNode4->X(),theNode4->Y(),theNode4->Z());
1938   SMDS_MeshFace* face;
1939   if(a.Distance(c) > b.Distance(d)){
1940     face = myTool->AddFace(theNode2, theNode4 , theNode1);
1941     theMeshDS->SetMeshElementOnShape(face, theFaceID );
1942     face = myTool->AddFace(theNode2, theNode3, theNode4);
1943     theMeshDS->SetMeshElementOnShape(face, theFaceID );
1944
1945   }
1946   else{
1947     face = myTool->AddFace(theNode1, theNode2 ,theNode3);
1948     theMeshDS->SetMeshElementOnShape(face, theFaceID );
1949     face = myTool->AddFace(theNode1, theNode3, theNode4);
1950     theMeshDS->SetMeshElementOnShape(face, theFaceID );
1951   }
1952 }