src/StdMeshers/StdMeshers_Quadrangle_2D.cxx

   1 //  Copyright (C) 2007-2010  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 //  Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 //  SMESH SMESH : implementaion of SMESH idl descriptions
  24 //  File   : StdMeshers_Quadrangle_2D.cxx
  25 //           Moved here from SMESH_Quadrangle_2D.cxx
  26 //  Author : Paul RASCLE, EDF
  27 //  Module : SMESH
  28 //
  29 #include "StdMeshers_Quadrangle_2D.hxx"
  30
  31 #include "StdMeshers_FaceSide.hxx"
  32
  33 #include "StdMeshers_QuadrangleParams.hxx"
  34
  35 #include "SMESH_Gen.hxx"
  36 #include "SMESH_Mesh.hxx"
  37 #include "SMESH_subMesh.hxx"
  38 #include "SMESH_MesherHelper.hxx"
  39 #include "SMESH_Block.hxx"
  40 #include "SMESH_Comment.hxx"
  41
  42 #include "SMDS_MeshElement.hxx"
  43 #include "SMDS_MeshNode.hxx"
  44 #include "SMDS_EdgePosition.hxx"
  45 #include "SMDS_FacePosition.hxx"
  46
  47 #include <BRep_Tool.hxx>
  48 #include <Geom_Surface.hxx>
  49 #include <NCollection_DefineArray2.hxx>
  50 #include <Precision.hxx>
  51 #include <TColStd_SequenceOfReal.hxx>
  52 #include <TColgp_SequenceOfXY.hxx>
  53 #include <TopExp.hxx>
  54 #include <TopExp_Explorer.hxx>
  55 #include <TopTools_ListIteratorOfListOfShape.hxx>
  56 #include <TopTools_MapOfShape.hxx>
  57 #include <TopoDS.hxx>
  58
  59 #include "utilities.h"
  60 #include "Utils_ExceptHandlers.hxx"
  61
  62 #ifndef StdMeshers_Array2OfNode_HeaderFile
  63 #define StdMeshers_Array2OfNode_HeaderFile
  64 typedef const SMDS_MeshNode* SMDS_MeshNodePtr;
  65 DEFINE_BASECOLLECTION (StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr)
  66 DEFINE_ARRAY2(StdMeshers_Array2OfNode,
  67               StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr)
  68 #endif
  69
  70 using namespace std;
  71
  72 typedef gp_XY gp_UV;
  73 typedef SMESH_Comment TComm;
  74
  75 //=============================================================================
  76 /*!
  77  *
  78  */
  79 //=============================================================================
  80
  81 StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId,
  82                                                     SMESH_Gen* gen)
  83      : SMESH_2D_Algo(hypId, studyId, gen)
  84 {
  85   MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D");
  86   _name = "Quadrangle_2D";
  87   _shapeType = (1 << TopAbs_FACE);
  88   _compatibleHypothesis.push_back("QuadrangleParams");
  89   _compatibleHypothesis.push_back("QuadranglePreference");
  90   _compatibleHypothesis.push_back("TrianglePreference");
  91   myTool = 0;
  92 }
  93
  94 //=============================================================================
  95 /*!
  96  *
  97  */
  98 //=============================================================================
  99
 100 StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D()
 101 {
 102   MESSAGE("StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D");
 103 }
 104
 105 //=============================================================================
 106 /*!
 107  *
 108  */
 109 //=============================================================================
 110
 111 bool StdMeshers_Quadrangle_2D::CheckHypothesis
 112                          (SMESH_Mesh&                          aMesh,
 113                           const TopoDS_Shape&                  aShape,
 114                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 115 {
 116   bool isOk = true;
 117   aStatus = SMESH_Hypothesis::HYP_OK;
 118
 119   const list <const SMESHDS_Hypothesis * >&hyps =
 120     GetUsedHypothesis(aMesh, aShape, false);
 121   const SMESHDS_Hypothesis *theHyp = 0;
 122
 123   if( hyps.size() == 1 ) {
 124     myTriaVertexID = -1;
 125     theHyp = hyps.front();
 126     if(strcmp("QuadrangleParams", theHyp->GetName()) == 0) {
 127       const StdMeshers_QuadrangleParams* theHyp1 =
 128         (const StdMeshers_QuadrangleParams*)theHyp;
 129       myTriaVertexID = theHyp1->GetTriaVertex();
 130       myQuadranglePreference= false;
 131       myTrianglePreference= false;
 132     }
 133     if(strcmp("QuadranglePreference", theHyp->GetName()) == 0) {
 134       myQuadranglePreference= true;
 135       myTrianglePreference= false;
 136       myTriaVertexID = -1;
 137     }
 138     else if(strcmp("TrianglePreference", theHyp->GetName()) == 0){
 139       myQuadranglePreference= false;
 140       myTrianglePreference= true;
 141       myTriaVertexID = -1;
 142     }
 143   }
 144
 145   else if( hyps.size() > 1 ) {
 146     theHyp = hyps.front();
 147     if(strcmp("QuadrangleParams", theHyp->GetName()) == 0) {
 148       const StdMeshers_QuadrangleParams* theHyp1 =
 149         (const StdMeshers_QuadrangleParams*)theHyp;
 150       myTriaVertexID = theHyp1->GetTriaVertex();
 151       theHyp = hyps.back();
 152       if(strcmp("QuadranglePreference", theHyp->GetName()) == 0) {
 153         myQuadranglePreference= true;
 154         myTrianglePreference= false;
 155       }
 156       else if(strcmp("TrianglePreference", theHyp->GetName()) == 0){
 157         myQuadranglePreference= false;
 158         myTrianglePreference= true;
 159       }
 160     }
 161     else {
 162       if(strcmp("QuadranglePreference", theHyp->GetName()) == 0) {
 163         myQuadranglePreference= true;
 164         myTrianglePreference= false;
 165       }
 166       else if(strcmp("TrianglePreference", theHyp->GetName()) == 0){
 167         myQuadranglePreference= false;
 168         myTrianglePreference= true;
 169       }
 170       const StdMeshers_QuadrangleParams* theHyp2 =
 171         (const StdMeshers_QuadrangleParams*)hyps.back();
 172       myTriaVertexID = theHyp2->GetTriaVertex();
 173     }
 174   }
 175
 176   else {
 177     myQuadranglePreference = false;
 178     myTrianglePreference = false;
 179     myTriaVertexID = -1;
 180   }
 181
 182   return isOk;
 183 }
 184
 185 //=============================================================================
 186 /*!
 187  *
 188  */
 189 //=============================================================================
 190
 191 bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh,
 192                                         const TopoDS_Shape& aShape)// throw (SALOME_Exception)
 193 {
 194   // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside
 195   //Unexpect aCatchSalomeException);
 196
 197   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 198   aMesh.GetSubMesh(aShape);
 199
 200   SMESH_MesherHelper helper(aMesh);
 201   myTool = &helper;
 202
 203   _quadraticMesh = myTool->IsQuadraticSubMesh(aShape);
 204
 205   FaceQuadStruct *quad = CheckNbEdges( aMesh, aShape );
 206   std::auto_ptr<FaceQuadStruct> quadDeleter( quad ); // to delete quad at exit from Compute()
 207   if (!quad)
 208     return false;
 209
 210   if(myQuadranglePreference) {
 211     int n1 = quad->side[0]->NbPoints();
 212     int n2 = quad->side[1]->NbPoints();
 213     int n3 = quad->side[2]->NbPoints();
 214     int n4 = quad->side[3]->NbPoints();
 215     int nfull = n1+n2+n3+n4;
 216     int ntmp = nfull/2;
 217     ntmp = ntmp*2;
 218     if( nfull==ntmp && ( (n1!=n3) || (n2!=n4) ) ) {
 219       // special path for using only quandrangle faces
 220       bool ok = ComputeQuadPref(aMesh, aShape, quad);
 221       return ok;
 222     }
 223   }
 224
 225   // set normalized grid on unit square in parametric domain
 226
 227   if (!SetNormalizedGrid(aMesh, aShape, quad))
 228     return false;
 229
 230   // --- compute 3D values on points, store points & quadrangles
 231
 232   int nbdown  = quad->side[0]->NbPoints();
 233   int nbup    = quad->side[2]->NbPoints();
 234
 235   int nbright = quad->side[1]->NbPoints();
 236   int nbleft  = quad->side[3]->NbPoints();
 237
 238   int nbhoriz  = Min(nbdown, nbup);
 239   int nbvertic = Min(nbright, nbleft);
 240
 241   const TopoDS_Face& F = TopoDS::Face(aShape);
 242   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
 243
 244   // internal mesh nodes
 245   int i, j, geomFaceID = meshDS->ShapeToIndex( F );
 246   for (i = 1; i < nbhoriz - 1; i++) {
 247     for (j = 1; j < nbvertic - 1; j++) {
 248       int ij = j * nbhoriz + i;
 249       double u = quad->uv_grid[ij].u;
 250       double v = quad->uv_grid[ij].v;
 251       gp_Pnt P = S->Value(u, v);
 252       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 253       meshDS->SetNodeOnFace(node, geomFaceID, u, v);
 254       quad->uv_grid[ij].node = node;
 255     }
 256   }
 257
 258   // mesh faces
 259
 260   //             [2]
 261   //      --.--.--.--.--.--  nbvertic
 262   //     |                 | ^
 263   //     |                 | ^
 264   // [3] |                 | ^ j  [1]
 265   //     |                 | ^
 266   //     |                 | ^
 267   //      ---.----.----.---  0
 268   //     0 > > > > > > > > nbhoriz
 269   //              i
 270   //             [0]
 271
 272   i = 0;
 273   int ilow = 0;
 274   int iup = nbhoriz - 1;
 275   if (quad->isEdgeOut[3]) { ilow++; } else { if (quad->isEdgeOut[1]) iup--; }
 276
 277   int jlow = 0;
 278   int jup = nbvertic - 1;
 279   if (quad->isEdgeOut[0]) { jlow++; } else { if (quad->isEdgeOut[2]) jup--; }
 280
 281   // regular quadrangles
 282   for (i = ilow; i < iup; i++) {
 283     for (j = jlow; j < jup; j++) {
 284       const SMDS_MeshNode *a, *b, *c, *d;
 285       a = quad->uv_grid[j * nbhoriz + i].node;
 286       b = quad->uv_grid[j * nbhoriz + i + 1].node;
 287       c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node;
 288       d = quad->uv_grid[(j + 1) * nbhoriz + i].node;
 289       SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
 290       if(face) {
 291         meshDS->SetMeshElementOnShape(face, geomFaceID);
 292       }
 293     }
 294   }
 295
 296   const vector<UVPtStruct>& uv_e0 = quad->side[0]->GetUVPtStruct(true,0 );
 297   const vector<UVPtStruct>& uv_e1 = quad->side[1]->GetUVPtStruct(false,1);
 298   const vector<UVPtStruct>& uv_e2 = quad->side[2]->GetUVPtStruct(true,1 );
 299   const vector<UVPtStruct>& uv_e3 = quad->side[3]->GetUVPtStruct(false,0);
 300
 301   if ( uv_e0.empty() || uv_e1.empty() || uv_e2.empty() || uv_e3.empty() )
 302     return error( COMPERR_BAD_INPUT_MESH );
 303
 304   double eps = Precision::Confusion();
 305
 306   // Boundary quadrangles
 307
 308   if (quad->isEdgeOut[0]) {
 309     // Down edge is out
 310     //
 311     // |___|___|___|___|___|___|
 312     // |   |   |   |   |   |   |
 313     // |___|___|___|___|___|___|
 314     // |   |   |   |   |   |   |
 315     // |___|___|___|___|___|___| __ first row of the regular grid
 316     // .  .  .  .  .  .  .  .  . __ down edge nodes
 317     //
 318     // >->->->->->->->->->->->-> -- direction of processing
 319
 320     int g = 0; // number of last processed node in the regular grid
 321
 322     // number of last node of the down edge to be processed
 323     int stop = nbdown - 1;
 324     // if right edge is out, we will stop at a node, previous to the last one
 325     if (quad->isEdgeOut[1]) stop--;
 326
 327     // for each node of the down edge find nearest node
 328     // in the first row of the regular grid and link them
 329     for (i = 0; i < stop; i++) {
 330       const SMDS_MeshNode *a, *b, *c, *d;
 331       a = uv_e0[i].node;
 332       b = uv_e0[i + 1].node;
 333       gp_Pnt pb (b->X(), b->Y(), b->Z());
 334
 335       // find node c in the regular grid, which will be linked with node b
 336       int near = g;
 337       if (i == stop - 1) {
 338         // right bound reached, link with the rightmost node
 339         near = iup;
 340         c = quad->uv_grid[nbhoriz + iup].node;
 341       }
 342       else {
 343         // find in the grid node c, nearest to the b
 344         double mind = RealLast();
 345         for (int k = g; k <= iup; k++) {
 346
 347           const SMDS_MeshNode *nk;
 348           if (k < ilow) // this can be, if left edge is out
 349             nk = uv_e3[1].node; // get node from the left edge
 350           else
 351             nk = quad->uv_grid[nbhoriz + k].node; // get one of middle nodes
 352
 353           gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 354           double dist = pb.Distance(pnk);
 355           if (dist < mind - eps) {
 356             c = nk;
 357             near = k;
 358             mind = dist;
 359           } else {
 360             break;
 361           }
 362         }
 363       }
 364
 365       if (near == g) { // make triangle
 366         SMDS_MeshFace* face = myTool->AddFace(a, b, c);
 367         if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 368       }
 369       else { // make quadrangle
 370         if (near - 1 < ilow)
 371           d = uv_e3[1].node;
 372         else
 373           d = quad->uv_grid[nbhoriz + near - 1].node;
 374         //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 375
 376         if(!myTrianglePreference){
 377           SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
 378           if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 379         }
 380         else {
 381           SplitQuad(meshDS, geomFaceID, a, b, c, d);
 382         }
 383
 384         // if node d is not at position g - make additional triangles
 385         if (near - 1 > g) {
 386           for (int k = near - 1; k > g; k--) {
 387             c = quad->uv_grid[nbhoriz + k].node;
 388             if (k - 1 < ilow)
 389               d = uv_e3[1].node;
 390             else
 391               d = quad->uv_grid[nbhoriz + k - 1].node;
 392             SMDS_MeshFace* face = myTool->AddFace(a, c, d);
 393             if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 394           }
 395         }
 396         g = near;
 397       }
 398     }
 399   } else {
 400     if (quad->isEdgeOut[2]) {
 401       // Up edge is out
 402       //
 403       // <-<-<-<-<-<-<-<-<-<-<-<-< -- direction of processing
 404       //
 405       // .  .  .  .  .  .  .  .  . __ up edge nodes
 406       //  ___ ___ ___ ___ ___ ___  __ first row of the regular grid
 407       // |   |   |   |   |   |   |
 408       // |___|___|___|___|___|___|
 409       // |   |   |   |   |   |   |
 410       // |___|___|___|___|___|___|
 411       // |   |   |   |   |   |   |
 412
 413       int g = nbhoriz - 1; // last processed node in the regular grid
 414
 415       int stop = 0;
 416       // if left edge is out, we will stop at a second node
 417       if (quad->isEdgeOut[3]) stop++;
 418
 419       // for each node of the up edge find nearest node
 420       // in the first row of the regular grid and link them
 421       for (i = nbup - 1; i > stop; i--) {
 422         const SMDS_MeshNode *a, *b, *c, *d;
 423         a = uv_e2[i].node;
 424         b = uv_e2[i - 1].node;
 425         gp_Pnt pb (b->X(), b->Y(), b->Z());
 426
 427         // find node c in the grid, which will be linked with node b
 428         int near = g;
 429         if (i == stop + 1) { // left bound reached, link with the leftmost node
 430           c = quad->uv_grid[nbhoriz*(nbvertic - 2) + ilow].node;
 431           near = ilow;
 432         } else {
 433           // find node c in the grid, nearest to the b
 434           double mind = RealLast();
 435           for (int k = g; k >= ilow; k--) {
 436             const SMDS_MeshNode *nk;
 437             if (k > iup)
 438               nk = uv_e1[nbright - 2].node;
 439             else
 440               nk = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
 441             gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 442             double dist = pb.Distance(pnk);
 443             if (dist < mind - eps) {
 444               c = nk;
 445               near = k;
 446               mind = dist;
 447             } else {
 448               break;
 449             }
 450           }
 451         }
 452
 453         if (near == g) { // make triangle
 454           SMDS_MeshFace* face = myTool->AddFace(a, b, c);
 455           if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 456         }
 457         else { // make quadrangle
 458           if (near + 1 > iup)
 459             d = uv_e1[nbright - 2].node;
 460           else
 461             d = quad->uv_grid[nbhoriz*(nbvertic - 2) + near + 1].node;
 462           //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 463           if(!myTrianglePreference){
 464             SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
 465             if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 466           }
 467           else {
 468             SplitQuad(meshDS, geomFaceID, a, b, c, d);
 469           }
 470
 471           if (near + 1 < g) { // if d not is at g - make additional triangles
 472             for (int k = near + 1; k < g; k++) {
 473               c = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
 474               if (k + 1 > iup)
 475                 d = uv_e1[nbright - 2].node;
 476               else
 477                 d = quad->uv_grid[nbhoriz*(nbvertic - 2) + k + 1].node;
 478               SMDS_MeshFace* face = myTool->AddFace(a, c, d);
 479               if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 480             }
 481           }
 482           g = near;
 483         }
 484       }
 485     }
 486   }
 487
 488   // right or left boundary quadrangles
 489   if (quad->isEdgeOut[1]) {
 490 //    MESSAGE("right edge is out");
 491     int g = 0; // last processed node in the grid
 492     int stop = nbright - 1;
 493     if (quad->isEdgeOut[2]) stop--;
 494     for (i = 0; i < stop; i++) {
 495       const SMDS_MeshNode *a, *b, *c, *d;
 496       a = uv_e1[i].node;
 497       b = uv_e1[i + 1].node;
 498       gp_Pnt pb (b->X(), b->Y(), b->Z());
 499
 500       // find node c in the grid, nearest to the b
 501       int near = g;
 502       if (i == stop - 1) { // up bondary reached
 503         c = quad->uv_grid[nbhoriz*(jup + 1) - 2].node;
 504         near = jup;
 505       } else {
 506         double mind = RealLast();
 507         for (int k = g; k <= jup; k++) {
 508           const SMDS_MeshNode *nk;
 509           if (k < jlow)
 510             nk = uv_e0[nbdown - 2].node;
 511           else
 512             nk = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
 513           gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 514           double dist = pb.Distance(pnk);
 515           if (dist < mind - eps) {
 516             c = nk;
 517             near = k;
 518             mind = dist;
 519           } else {
 520             break;
 521           }
 522         }
 523       }
 524
 525       if (near == g) { // make triangle
 526         SMDS_MeshFace* face = myTool->AddFace(a, b, c);
 527         if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 528       }
 529       else { // make quadrangle
 530         if (near - 1 < jlow)
 531           d = uv_e0[nbdown - 2].node;
 532         else
 533           d = quad->uv_grid[nbhoriz*near - 2].node;
 534         //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 535
 536         if(!myTrianglePreference){
 537           SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
 538           if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 539         }
 540         else {
 541           SplitQuad(meshDS, geomFaceID, a, b, c, d);
 542         }
 543
 544         if (near - 1 > g) { // if d not is at g - make additional triangles
 545           for (int k = near - 1; k > g; k--) {
 546             c = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
 547             if (k - 1 < jlow)
 548               d = uv_e0[nbdown - 2].node;
 549             else
 550               d = quad->uv_grid[nbhoriz*k - 2].node;
 551             SMDS_MeshFace* face = myTool->AddFace(a, c, d);
 552             if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 553           }
 554         }
 555         g = near;
 556       }
 557     }
 558   } else {
 559     if (quad->isEdgeOut[3]) {
 560 //      MESSAGE("left edge is out");
 561       int g = nbvertic - 1; // last processed node in the grid
 562       int stop = 0;
 563       if (quad->isEdgeOut[0]) stop++;
 564       for (i = nbleft - 1; i > stop; i--) {
 565         const SMDS_MeshNode *a, *b, *c, *d;
 566         a = uv_e3[i].node;
 567         b = uv_e3[i - 1].node;
 568         gp_Pnt pb (b->X(), b->Y(), b->Z());
 569
 570         // find node c in the grid, nearest to the b
 571         int near = g;
 572         if (i == stop + 1) { // down bondary reached
 573           c = quad->uv_grid[nbhoriz*jlow + 1].node;
 574           near = jlow;
 575         } else {
 576           double mind = RealLast();
 577           for (int k = g; k >= jlow; k--) {
 578             const SMDS_MeshNode *nk;
 579             if (k > jup)
 580               nk = uv_e2[1].node;
 581             else
 582               nk = quad->uv_grid[nbhoriz*k + 1].node;
 583             gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 584             double dist = pb.Distance(pnk);
 585             if (dist < mind - eps) {
 586               c = nk;
 587               near = k;
 588               mind = dist;
 589             } else {
 590               break;
 591             }
 592           }
 593         }
 594
 595         if (near == g) { // make triangle
 596           SMDS_MeshFace* face = myTool->AddFace(a, b, c);
 597           if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 598         }
 599         else { // make quadrangle
 600           if (near + 1 > jup)
 601             d = uv_e2[1].node;
 602           else
 603             d = quad->uv_grid[nbhoriz*(near + 1) + 1].node;
 604           //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 605           if(!myTrianglePreference){
 606             SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
 607             if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 608           }
 609           else {
 610             SplitQuad(meshDS, geomFaceID, a, b, c, d);
 611           }
 612
 613           if (near + 1 < g) { // if d not is at g - make additional triangles
 614             for (int k = near + 1; k < g; k++) {
 615               c = quad->uv_grid[nbhoriz*k + 1].node;
 616               if (k + 1 > jup)
 617                 d = uv_e2[1].node;
 618               else
 619                 d = quad->uv_grid[nbhoriz*(k + 1) + 1].node;
 620               SMDS_MeshFace* face = myTool->AddFace(a, c, d);
 621               if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
 622             }
 623           }
 624           g = near;
 625         }
 626       }
 627     }
 628   }
 629
 630   bool isOk = true;
 631   return isOk;
 632 }
 633
 634
 635 //=============================================================================
 636 /*!
 637  *  Evaluate
 638  */
 639 //=============================================================================
 640
 641 bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh& aMesh,
 642                                         const TopoDS_Shape& aShape,
 643                                         MapShapeNbElems& aResMap)
 644
 645 {
 646   aMesh.GetSubMesh(aShape);
 647
 648   std::vector<int> aNbNodes(4);
 649   bool IsQuadratic = false;
 650   if( !CheckNbEdgesForEvaluate( aMesh, aShape, aResMap, aNbNodes, IsQuadratic ) ) {
 651     std::vector<int> aResVec(SMDSEntity_Last);
 652     for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
 653     SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 654     aResMap.insert(std::make_pair(sm,aResVec));
 655     SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 656     smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
 657     return false;
 658   }
 659
 660   if(myQuadranglePreference) {
 661     int n1 = aNbNodes[0];
 662     int n2 = aNbNodes[1];
 663     int n3 = aNbNodes[2];
 664     int n4 = aNbNodes[3];
 665     int nfull = n1+n2+n3+n4;
 666     int ntmp = nfull/2;
 667     ntmp = ntmp*2;
 668     if( nfull==ntmp && ( (n1!=n3) || (n2!=n4) ) ) {
 669       // special path for using only quandrangle faces
 670       return EvaluateQuadPref(aMesh, aShape, aNbNodes, aResMap, IsQuadratic);
 671       //return true;
 672     }
 673   }
 674
 675   int nbdown  = aNbNodes[0];
 676   int nbup    = aNbNodes[2];
 677
 678   int nbright = aNbNodes[1];
 679   int nbleft  = aNbNodes[3];
 680
 681   int nbhoriz  = Min(nbdown, nbup);
 682   int nbvertic = Min(nbright, nbleft);
 683
 684   int dh = Max(nbdown, nbup) - nbhoriz;
 685   int dv = Max(nbright, nbleft) - nbvertic;
 686
 687   //int kdh = 0;
 688   //if(dh>0) kdh = 1;
 689   //int kdv = 0;
 690   //if(dv>0) kdv = 1;
 691
 692   int nbNodes = (nbhoriz-2)*(nbvertic-2);
 693   //int nbFaces3 = dh + dv + kdh*(nbvertic-1)*2 + kdv*(nbhoriz-1)*2;
 694   int nbFaces3 = dh + dv;
 695   //if( kdh==1 && kdv==1 ) nbFaces3 -= 2;
 696   //if( dh>0 && dv>0 ) nbFaces3 -= 2;
 697   //int nbFaces4 = (nbhoriz-1-kdh)*(nbvertic-1-kdv);
 698   int nbFaces4 = (nbhoriz-1)*(nbvertic-1);
 699
 700   std::vector<int> aVec(SMDSEntity_Last);
 701   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i] = 0;
 702   if(IsQuadratic) {
 703     aVec[SMDSEntity_Quad_Triangle] = nbFaces3;
 704     aVec[SMDSEntity_Quad_Quadrangle] = nbFaces4;
 705     int nbbndedges = nbdown + nbup + nbright + nbleft -4;
 706     int nbintedges = ( nbFaces4*4 + nbFaces3*3 - nbbndedges ) / 2;
 707     aVec[SMDSEntity_Node] = nbNodes + nbintedges;
 708     if( aNbNodes.size()==5 ) {
 709       aVec[SMDSEntity_Quad_Triangle] = nbFaces3 + aNbNodes[3] -1;
 710       aVec[SMDSEntity_Quad_Quadrangle] = nbFaces4 - aNbNodes[3] +1;
 711     }
 712   }
 713   else {
 714     aVec[SMDSEntity_Node] = nbNodes;
 715     aVec[SMDSEntity_Triangle] = nbFaces3;
 716     aVec[SMDSEntity_Quadrangle] = nbFaces4;
 717     if( aNbNodes.size()==5 ) {
 718       aVec[SMDSEntity_Triangle] = nbFaces3 + aNbNodes[3] - 1;
 719       aVec[SMDSEntity_Quadrangle] = nbFaces4 - aNbNodes[3] + 1;
 720     }
 721   }
 722   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 723   aResMap.insert(std::make_pair(sm,aVec));
 724
 725   return true;
 726 }
 727
 728
 729 //================================================================================
 730 /*!
 731  * \brief Return true if only two given edges meat at their common vertex
 732  */
 733 //================================================================================
 734
 735 static bool twoEdgesMeatAtVertex(const TopoDS_Edge& e1,
 736                                  const TopoDS_Edge& e2,
 737                                  SMESH_Mesh &       mesh)
 738 {
 739   TopoDS_Vertex v;
 740   if ( !TopExp::CommonVertex( e1, e2, v ))
 741     return false;
 742   TopTools_ListIteratorOfListOfShape ancestIt( mesh.GetAncestors( v ));
 743   for ( ; ancestIt.More() ; ancestIt.Next() )
 744     if ( ancestIt.Value().ShapeType() == TopAbs_EDGE )
 745       if ( !e1.IsSame( ancestIt.Value() ) && !e2.IsSame( ancestIt.Value() ))
 746         return false;
 747   return true;
 748 }
 749
 750 //=============================================================================
 751 /*!
 752  *
 753  */
 754 //=============================================================================
 755
 756 FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh &         aMesh,
 757                                                        const TopoDS_Shape & aShape)
 758   //throw(SALOME_Exception)
 759 {
 760   TopoDS_Face F = TopoDS::Face(aShape);
 761   if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD );
 762   const bool ignoreMediumNodes = _quadraticMesh;
 763
 764   // verify 1 wire only, with 4 edges
 765   TopoDS_Vertex V;
 766   list< TopoDS_Edge > edges;
 767   list< int > nbEdgesInWire;
 768   int nbWire = SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
 769   if (nbWire != 1) {
 770     error(COMPERR_BAD_SHAPE, TComm("Wrong number of wires: ") << nbWire);
 771     return 0;
 772   }
 773   FaceQuadStruct* quad = new FaceQuadStruct;
 774   quad->uv_grid = 0;
 775   quad->side.reserve(nbEdgesInWire.front());
 776
 777   int nbSides = 0;
 778   list< TopoDS_Edge >::iterator edgeIt = edges.begin();
 779   if ( nbEdgesInWire.front() == 3 ) // exactly 3 edges
 780   {
 781     SMESH_Comment comment;
 782     SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 783     if ( myTriaVertexID == -1)
 784     {
 785       comment << "No Base vertex parameter provided for a trilateral geometrical face";
 786     }
 787     else
 788     {
 789       TopoDS_Vertex V = TopoDS::Vertex(meshDS->IndexToShape(myTriaVertexID));
 790       if ( !V.IsNull() ) {
 791         TopoDS_Edge E1,E2,E3;
 792         for(; edgeIt != edges.end(); ++edgeIt) {
 793           TopoDS_Edge E =  *edgeIt;
 794           TopoDS_Vertex VF, VL;
 795           TopExp::Vertices(E, VF, VL, true);
 796           if( VF.IsSame(V) )
 797             E1 = E;
 798           else if( VL.IsSame(V) )
 799             E3 = E;
 800           else
 801             E2 = E;
 802         }
 803         if ( !E1.IsNull() && !E2.IsNull() && !E3.IsNull() )
 804         {
 805           quad->side.push_back( new StdMeshers_FaceSide(F, E1, &aMesh, true, ignoreMediumNodes));
 806           quad->side.push_back( new StdMeshers_FaceSide(F, E2, &aMesh, true, ignoreMediumNodes));
 807           quad->side.push_back( new StdMeshers_FaceSide(F, E3, &aMesh, false,ignoreMediumNodes));
 808           const vector<UVPtStruct>& UVPSleft  = quad->side[0]->GetUVPtStruct(true,0);
 809           /*  vector<UVPtStruct>& UVPStop   = */quad->side[1]->GetUVPtStruct(false,1);
 810           /*  vector<UVPtStruct>& UVPSright = */quad->side[2]->GetUVPtStruct(true,1);
 811           const SMDS_MeshNode* aNode = UVPSleft[0].node;
 812           gp_Pnt2d aPnt2d( UVPSleft[0].u, UVPSleft[0].v );
 813           quad->side.push_back( new StdMeshers_FaceSide(aNode, aPnt2d, quad->side[1]));
 814           return quad;
 815         }
 816       }
 817       comment << "Invalid Base vertex parameter: " << myTriaVertexID << " is not among [";
 818       TopTools_MapOfShape vMap;
 819       for ( TopExp_Explorer v( aShape, TopAbs_VERTEX ); v.More(); v.Next())
 820         if ( vMap.Add( v.Current() ))
 821           comment << meshDS->ShapeToIndex( v.Current() ) << ( vMap.Extent()==3 ? "]" : ", ");
 822     }
 823     error( comment );
 824     delete quad;
 825     return quad = 0;
 826   }
 827   else if ( nbEdgesInWire.front() == 4 ) { // exactly 4 edges
 828     for ( ; edgeIt != edges.end(); ++edgeIt, nbSides++ )
 829       quad->side.push_back( new StdMeshers_FaceSide(F, *edgeIt, &aMesh,
 830                                                     nbSides<TOP_SIDE, ignoreMediumNodes));
 831   }
 832   else if ( nbEdgesInWire.front() > 4 ) { // more than 4 edges - try to unite some
 833     list< TopoDS_Edge > sideEdges;
 834     while ( !edges.empty()) {
 835       sideEdges.clear();
 836       sideEdges.splice( sideEdges.end(), edges, edges.begin()); // edges.front() -> sideEdges.end()
 837       bool sameSide = true;
 838       while ( !edges.empty() && sameSide ) {
 839         sameSide = SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() );
 840         if ( sameSide )
 841           sideEdges.splice( sideEdges.end(), edges, edges.begin());
 842       }
 843       if ( nbSides == 0 ) { // go backward from the first edge
 844         sameSide = true;
 845         while ( !edges.empty() && sameSide ) {
 846           sameSide = SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() );
 847           if ( sameSide )
 848             sideEdges.splice( sideEdges.begin(), edges, --edges.end());
 849         }
 850       }
 851       quad->side.push_back( new StdMeshers_FaceSide(F, sideEdges, &aMesh,
 852                                                     nbSides<TOP_SIDE, ignoreMediumNodes));
 853       ++nbSides;
 854     }
 855     // issue 20222. Try to unite only edges shared by two same faces
 856     if (nbSides < 4) {
 857       // delete found sides
 858       { FaceQuadStruct cleaner( *quad ); }
 859       quad->side.clear();
 860       quad->side.reserve(nbEdgesInWire.front());
 861       nbSides = 0;
 862
 863       SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
 864       while ( !edges.empty()) {
 865         sideEdges.clear();
 866         sideEdges.splice( sideEdges.end(), edges, edges.begin());
 867         bool sameSide = true;
 868         while ( !edges.empty() && sameSide ) {
 869           sameSide =
 870             SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() ) &&
 871             twoEdgesMeatAtVertex( sideEdges.back(), edges.front(), aMesh );
 872           if ( sameSide )
 873             sideEdges.splice( sideEdges.end(), edges, edges.begin());
 874         }
 875         if ( nbSides == 0 ) { // go backward from the first edge
 876           sameSide = true;
 877           while ( !edges.empty() && sameSide ) {
 878             sameSide =
 879               SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() ) &&
 880               twoEdgesMeatAtVertex( sideEdges.front(), edges.back(), aMesh );
 881             if ( sameSide )
 882               sideEdges.splice( sideEdges.begin(), edges, --edges.end());
 883           }
 884         }
 885         quad->side.push_back( new StdMeshers_FaceSide(F, sideEdges, &aMesh,
 886                                                       nbSides<TOP_SIDE, ignoreMediumNodes));
 887         ++nbSides;
 888       }
 889     }
 890   }
 891   if (nbSides != 4) {
 892 #ifdef _DEBUG_
 893     MESSAGE ( "StdMeshers_Quadrangle_2D. Edge IDs of " << nbSides << " sides:\n" );
 894     for ( int i = 0; i < nbSides; ++i ) {
 895       MESSAGE ( " ( " );
 896       for ( int e = 0; e < quad->side[i]->NbEdges(); ++e )
 897         MESSAGE ( myTool->GetMeshDS()->ShapeToIndex( quad->side[i]->Edge( e )) << " " );
 898       MESSAGE ( ")\n" );
 899     }
 900     //cout << endl;
 901 #endif
 902     if ( !nbSides )
 903       nbSides = nbEdgesInWire.front();
 904     error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides);
 905     delete quad;
 906     quad = 0;
 907   }
 908
 909   return quad;
 910 }
 911
 912
 913 //=============================================================================
 914 /*!
 915  *
 916  */
 917 //=============================================================================
 918
 919 bool StdMeshers_Quadrangle_2D::CheckNbEdgesForEvaluate(SMESH_Mesh& aMesh,
 920                                                        const TopoDS_Shape & aShape,
 921                                                        MapShapeNbElems& aResMap,
 922                                                        std::vector<int>& aNbNodes,
 923                                                        bool& IsQuadratic)
 924
 925 {
 926   const TopoDS_Face & F = TopoDS::Face(aShape);
 927
 928   // verify 1 wire only, with 4 edges
 929   TopoDS_Vertex V;
 930   list< TopoDS_Edge > edges;
 931   list< int > nbEdgesInWire;
 932   int nbWire = SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
 933   if (nbWire != 1) {
 934     return false;
 935   }
 936
 937   aNbNodes.resize(4);
 938
 939   int nbSides = 0;
 940   list< TopoDS_Edge >::iterator edgeIt = edges.begin();
 941   SMESH_subMesh * sm = aMesh.GetSubMesh( *edgeIt );
 942   MapShapeNbElemsItr anIt = aResMap.find(sm);
 943   if(anIt==aResMap.end()) {
 944     return false;
 945   }
 946   std::vector<int> aVec = (*anIt).second;
 947   IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
 948   if ( nbEdgesInWire.front() == 3 ) { // exactly 3 edges
 949     if(myTriaVertexID>0) {
 950       SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 951       TopoDS_Vertex V = TopoDS::Vertex(meshDS->IndexToShape(myTriaVertexID));
 952       if(!V.IsNull()) {
 953         TopoDS_Edge E1,E2,E3;
 954         for(; edgeIt != edges.end(); ++edgeIt) {
 955           TopoDS_Edge E =  TopoDS::Edge(*edgeIt);
 956           TopoDS_Vertex VF, VL;
 957           TopExp::Vertices(E, VF, VL, true);
 958           if( VF.IsSame(V) )
 959             E1 = E;
 960           else if( VL.IsSame(V) )
 961             E3 = E;
 962           else
 963             E2 = E;
 964         }
 965         SMESH_subMesh * sm = aMesh.GetSubMesh(E1);
 966         MapShapeNbElemsItr anIt = aResMap.find(sm);
 967         if(anIt==aResMap.end()) return false;
 968         std::vector<int> aVec = (*anIt).second;
 969         if(IsQuadratic)
 970           aNbNodes[0] = (aVec[SMDSEntity_Node]-1)/2 + 2;
 971         else
 972           aNbNodes[0] = aVec[SMDSEntity_Node] + 2;
 973         sm = aMesh.GetSubMesh(E2);
 974         anIt = aResMap.find(sm);
 975         if(anIt==aResMap.end()) return false;
 976         aVec = (*anIt).second;
 977         if(IsQuadratic)
 978           aNbNodes[1] = (aVec[SMDSEntity_Node]-1)/2 + 2;
 979         else
 980           aNbNodes[1] = aVec[SMDSEntity_Node] + 2;
 981         sm = aMesh.GetSubMesh(E3);
 982         anIt = aResMap.find(sm);
 983         if(anIt==aResMap.end()) return false;
 984         aVec = (*anIt).second;
 985         if(IsQuadratic)
 986           aNbNodes[2] = (aVec[SMDSEntity_Node]-1)/2 + 2;
 987         else
 988           aNbNodes[2] = aVec[SMDSEntity_Node] + 2;
 989         aNbNodes[3] = aNbNodes[1];
 990         aNbNodes.resize(5);
 991         nbSides = 4;
 992       }
 993     }
 994   }
 995   if ( nbEdgesInWire.front() == 4 ) { // exactly 4 edges
 996     for(; edgeIt != edges.end(); edgeIt++) {
 997       SMESH_subMesh * sm = aMesh.GetSubMesh( *edgeIt );
 998       MapShapeNbElemsItr anIt = aResMap.find(sm);
 999       if(anIt==aResMap.end()) {
1000         return false;
1001       }
1002       std::vector<int> aVec = (*anIt).second;
1003       if(IsQuadratic)
1004         aNbNodes[nbSides] = (aVec[SMDSEntity_Node]-1)/2 + 2;
1005       else
1006         aNbNodes[nbSides] = aVec[SMDSEntity_Node] + 2;
1007       nbSides++;
1008     }
1009   }
1010   else if ( nbEdgesInWire.front() > 4 ) { // more than 4 edges - try to unite some
1011     list< TopoDS_Edge > sideEdges;
1012     while ( !edges.empty()) {
1013       sideEdges.clear();
1014       sideEdges.splice( sideEdges.end(), edges, edges.begin()); // edges.front() -> sideEdges.end()
1015       bool sameSide = true;
1016       while ( !edges.empty() && sameSide ) {
1017         sameSide = SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() );
1018         if ( sameSide )
1019           sideEdges.splice( sideEdges.end(), edges, edges.begin());
1020       }
1021       if ( nbSides == 0 ) { // go backward from the first edge
1022         sameSide = true;
1023         while ( !edges.empty() && sameSide ) {
1024           sameSide = SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() );
1025           if ( sameSide )
1026             sideEdges.splice( sideEdges.begin(), edges, --edges.end());
1027         }
1028       }
1029       list<TopoDS_Edge>::iterator ite = sideEdges.begin();
1030       aNbNodes[nbSides] = 1;
1031       for(; ite!=sideEdges.end(); ite++) {
1032         SMESH_subMesh * sm = aMesh.GetSubMesh( *ite );
1033         MapShapeNbElemsItr anIt = aResMap.find(sm);
1034         if(anIt==aResMap.end()) {
1035           return false;
1036         }
1037         std::vector<int> aVec = (*anIt).second;
1038         if(IsQuadratic)
1039           aNbNodes[nbSides] += (aVec[SMDSEntity_Node]-1)/2 + 1;
1040         else
1041           aNbNodes[nbSides] += aVec[SMDSEntity_Node] + 1;
1042       }
1043       ++nbSides;
1044     }
1045     // issue 20222. Try to unite only edges shared by two same faces
1046     if (nbSides < 4) {
1047       nbSides = 0;
1048       SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
1049       while ( !edges.empty()) {
1050         sideEdges.clear();
1051         sideEdges.splice( sideEdges.end(), edges, edges.begin());
1052         bool sameSide = true;
1053         while ( !edges.empty() && sameSide ) {
1054           sameSide =
1055             SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() ) &&
1056             twoEdgesMeatAtVertex( sideEdges.back(), edges.front(), aMesh );
1057           if ( sameSide )
1058             sideEdges.splice( sideEdges.end(), edges, edges.begin());
1059         }
1060         if ( nbSides == 0 ) { // go backward from the first edge
1061           sameSide = true;
1062           while ( !edges.empty() && sameSide ) {
1063             sameSide =
1064               SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() ) &&
1065               twoEdgesMeatAtVertex( sideEdges.front(), edges.back(), aMesh );
1066             if ( sameSide )
1067               sideEdges.splice( sideEdges.begin(), edges, --edges.end());
1068           }
1069         }
1070         list<TopoDS_Edge>::iterator ite = sideEdges.begin();
1071         aNbNodes[nbSides] = 1;
1072         for(; ite!=sideEdges.end(); ite++) {
1073           SMESH_subMesh * sm = aMesh.GetSubMesh( *ite );
1074           MapShapeNbElemsItr anIt = aResMap.find(sm);
1075           if(anIt==aResMap.end()) {
1076             return false;
1077           }
1078           std::vector<int> aVec = (*anIt).second;
1079           if(IsQuadratic)
1080             aNbNodes[nbSides] += (aVec[SMDSEntity_Node]-1)/2 + 1;
1081           else
1082             aNbNodes[nbSides] += aVec[SMDSEntity_Node] + 1;
1083         }
1084         ++nbSides;
1085       }
1086     }
1087   }
1088   if (nbSides != 4) {
1089     if ( !nbSides )
1090       nbSides = nbEdgesInWire.front();
1091     error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides);
1092     return false;
1093   }
1094
1095   return true;
1096 }
1097
1098
1099 //=============================================================================
1100 /*!
1101  *  CheckAnd2Dcompute
1102  */
1103 //=============================================================================
1104
1105 FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute
1106                            (SMESH_Mesh &         aMesh,
1107                             const TopoDS_Shape & aShape,
1108                             const bool           CreateQuadratic) //throw(SALOME_Exception)
1109 {
1110   _quadraticMesh = CreateQuadratic;
1111
1112   FaceQuadStruct *quad = CheckNbEdges(aMesh, aShape);
1113
1114   if(!quad) return 0;
1115
1116   // set normalized grid on unit square in parametric domain
1117   bool stat = SetNormalizedGrid(aMesh, aShape, quad);
1118   if(!stat) {
1119     if(!quad)
1120       delete quad;
1121     quad = 0;
1122   }
1123
1124   return quad;
1125 }
1126
1127 //=============================================================================
1128 /*!
1129  *
1130  */
1131 //=============================================================================
1132
1133 faceQuadStruct::~faceQuadStruct()
1134 {
1135   for (int i = 0; i < side.size(); i++) {
1136     if (side[i])     delete side[i];
1137   }
1138   if (uv_grid)       delete [] uv_grid;
1139 }
1140
1141 namespace {
1142   inline const vector<UVPtStruct>& GetUVPtStructIn(FaceQuadStruct* quad, int i, int nbSeg)
1143   {
1144     bool   isXConst   = ( i == BOTTOM_SIDE || i == TOP_SIDE );
1145     double constValue = ( i == BOTTOM_SIDE || i == LEFT_SIDE ) ? 0 : 1;
1146     return
1147       quad->isEdgeOut[i] ?
1148       quad->side[i]->SimulateUVPtStruct(nbSeg,isXConst,constValue) :
1149       quad->side[i]->GetUVPtStruct(isXConst,constValue);
1150   }
1151 }
1152
1153 //=============================================================================
1154 /*!
1155  *
1156  */
1157 //=============================================================================
1158
1159 bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh,
1160                                                   const TopoDS_Shape& aShape,
1161                                                   FaceQuadStruct* & quad) //throw (SALOME_Exception)
1162 {
1163   // Algorithme décrit dans "Génération automatique de maillages"
1164   // P.L. GEORGE, MASSON, § 6.4.1 p. 84-85
1165   // traitement dans le domaine paramétrique 2d u,v
1166   // transport - projection sur le carré unité
1167
1168 //  MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid");
1169 //  const TopoDS_Face& F = TopoDS::Face(aShape);
1170
1171   // 1 --- find orientation of the 4 edges, by test on extrema
1172
1173   //      max             min                    0     x1     1
1174   //     |<----north-2-------^                a3 -------------> a2
1175   //     |                   |                   ^1          1^
1176   //    west-3            east-1 =right          |            |
1177   //     |                   |         ==>       |            |
1178   //  y0 |                   | y1                |            |
1179   //     |                   |                   |0          0|
1180   //     v----south-0-------->                a0 -------------> a1
1181   //      min             max                    0     x0     1
1182   //             =down
1183   //
1184
1185   // 3 --- 2D normalized values on unit square [0..1][0..1]
1186
1187   int nbhoriz  = Min(quad->side[0]->NbPoints(), quad->side[2]->NbPoints());
1188   int nbvertic = Min(quad->side[1]->NbPoints(), quad->side[3]->NbPoints());
1189
1190   quad->isEdgeOut[0] = (quad->side[0]->NbPoints() > quad->side[2]->NbPoints());
1191   quad->isEdgeOut[1] = (quad->side[1]->NbPoints() > quad->side[3]->NbPoints());
1192   quad->isEdgeOut[2] = (quad->side[2]->NbPoints() > quad->side[0]->NbPoints());
1193   quad->isEdgeOut[3] = (quad->side[3]->NbPoints() > quad->side[1]->NbPoints());
1194
1195   UVPtStruct *uv_grid = quad->uv_grid = new UVPtStruct[nbvertic * nbhoriz];
1196
1197   const vector<UVPtStruct>& uv_e0 = GetUVPtStructIn( quad, 0, nbhoriz - 1 );
1198   const vector<UVPtStruct>& uv_e1 = GetUVPtStructIn( quad, 1, nbvertic - 1 );
1199   const vector<UVPtStruct>& uv_e2 = GetUVPtStructIn( quad, 2, nbhoriz - 1 );
1200   const vector<UVPtStruct>& uv_e3 = GetUVPtStructIn( quad, 3, nbvertic - 1 );
1201
1202   if ( uv_e0.empty() || uv_e1.empty() || uv_e2.empty() || uv_e3.empty() )
1203     //return error( "Can't find nodes on sides");
1204     return error( COMPERR_BAD_INPUT_MESH );
1205
1206   // nodes Id on "in" edges
1207   if (! quad->isEdgeOut[0]) {
1208     int j = 0;
1209     for (int i = 0; i < nbhoriz; i++) { // down
1210       int ij = j * nbhoriz + i;
1211       uv_grid[ij].node = uv_e0[i].node;
1212     }
1213   }
1214   if (! quad->isEdgeOut[1]) {
1215     int i = nbhoriz - 1;
1216     for (int j = 0; j < nbvertic; j++) { // right
1217       int ij = j * nbhoriz + i;
1218       uv_grid[ij].node = uv_e1[j].node;
1219     }
1220   }
1221   if (! quad->isEdgeOut[2]) {
1222     int j = nbvertic - 1;
1223     for (int i = 0; i < nbhoriz; i++) { // up
1224       int ij = j * nbhoriz + i;
1225       uv_grid[ij].node = uv_e2[i].node;
1226     }
1227   }
1228   if (! quad->isEdgeOut[3]) {
1229     int i = 0;
1230     for (int j = 0; j < nbvertic; j++) { // left
1231       int ij = j * nbhoriz + i;
1232       uv_grid[ij].node = uv_e3[j].node;
1233     }
1234   }
1235
1236   // normalized 2d values on grid
1237   for (int i = 0; i < nbhoriz; i++) {
1238     for (int j = 0; j < nbvertic; j++) {
1239       int ij = j * nbhoriz + i;
1240       // --- droite i cste : x = x0 + y(x1-x0)
1241       double x0 = uv_e0[i].normParam;   // bas - sud
1242       double x1 = uv_e2[i].normParam;   // haut - nord
1243       // --- droite j cste : y = y0 + x(y1-y0)
1244       double y0 = uv_e3[j].normParam;   // gauche-ouest
1245       double y1 = uv_e1[j].normParam;   // droite - est
1246       // --- intersection : x=x0+(y0+x(y1-y0))(x1-x0)
1247       double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
1248       double y = y0 + x * (y1 - y0);
1249       uv_grid[ij].x = x;
1250       uv_grid[ij].y = y;
1251       //MESSAGE("-xy-01 "<<x0<<" "<<x1<<" "<<y0<<" "<<y1);
1252       //MESSAGE("-xy-norm "<<i<<" "<<j<<" "<<x<<" "<<y);
1253     }
1254   }
1255
1256   // 4 --- projection on 2d domain (u,v)
1257   gp_UV a0( uv_e0.front().u, uv_e0.front().v );
1258   gp_UV a1( uv_e0.back().u,  uv_e0.back().v );
1259   gp_UV a2( uv_e2.back().u,  uv_e2.back().v );
1260   gp_UV a3( uv_e2.front().u, uv_e2.front().v );
1261
1262   for (int i = 0; i < nbhoriz; i++) {
1263     for (int j = 0; j < nbvertic; j++) {
1264       int ij = j * nbhoriz + i;
1265       double x = uv_grid[ij].x;
1266       double y = uv_grid[ij].y;
1267       double param_0 = uv_e0[0].normParam + x * (uv_e0.back().normParam - uv_e0[0].normParam); // sud
1268       double param_2 = uv_e2[0].normParam + x * (uv_e2.back().normParam - uv_e2[0].normParam); // nord
1269       double param_1 = uv_e1[0].normParam + y * (uv_e1.back().normParam - uv_e1[0].normParam); // est
1270       double param_3 = uv_e3[0].normParam + y * (uv_e3.back().normParam - uv_e3[0].normParam); // ouest
1271
1272       //MESSAGE("params "<<param_0<<" "<<param_1<<" "<<param_2<<" "<<param_3);
1273       gp_UV p0 = quad->side[0]->Value2d(param_0).XY();
1274       gp_UV p1 = quad->side[1]->Value2d(param_1).XY();
1275       gp_UV p2 = quad->side[2]->Value2d(param_2).XY();
1276       gp_UV p3 = quad->side[3]->Value2d(param_3).XY();
1277
1278       gp_UV uv = (1 - y) * p0 + x * p1 + y * p2 + (1 - x) * p3;
1279       uv -= (1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3;
1280
1281       uv_grid[ij].u = uv.X();
1282       uv_grid[ij].v = uv.Y();
1283     }
1284   }
1285   return true;
1286 }
1287
1288 //=======================================================================
1289 //function : ShiftQuad
1290 //purpose  : auxilary function for ComputeQuadPref
1291 //=======================================================================
1292
1293 static void ShiftQuad(FaceQuadStruct* quad, const int num, bool)
1294 {
1295   StdMeshers_FaceSide* side[4] = { quad->side[0], quad->side[1], quad->side[2], quad->side[3] };
1296   for (int i = BOTTOM_SIDE; i < NB_SIDES; ++i ) {
1297     int id = ( i + num ) % NB_SIDES;
1298     bool wasForward = ( i < TOP_SIDE );
1299     bool newForward = ( id < TOP_SIDE );
1300     if ( wasForward != newForward )
1301       side[ i ]->Reverse();
1302     quad->side[ id ] = side[ i ];
1303   }
1304 }
1305
1306 //=======================================================================
1307 //function : CalcUV
1308 //purpose  : auxilary function for ComputeQuadPref
1309 //=======================================================================
1310
1311 static gp_UV CalcUV(double x0, double x1, double y0, double y1,
1312                     FaceQuadStruct* quad,
1313                     const gp_UV& a0, const gp_UV& a1,
1314                     const gp_UV& a2, const gp_UV& a3)
1315 {
1316   const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0 );
1317   const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
1318   const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1 );
1319   const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
1320
1321   double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
1322   double y = y0 + x * (y1 - y0);
1323
1324   double param_b = uv_eb[0].normParam + x * (uv_eb.back().normParam - uv_eb[0].normParam);
1325   double param_t = uv_et[0].normParam + x * (uv_et.back().normParam - uv_et[0].normParam);
1326   double param_r = uv_er[0].normParam + y * (uv_er.back().normParam - uv_er[0].normParam);
1327   double param_l = uv_el[0].normParam + y * (uv_el.back().normParam - uv_el[0].normParam);
1328
1329   gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(param_b).XY();
1330   gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(param_r).XY();
1331   gp_UV p2 = quad->side[TOP_SIDE   ]->Value2d(param_t).XY();
1332   gp_UV p3 = quad->side[LEFT_SIDE  ]->Value2d(param_l).XY();
1333
1334   gp_UV uv = p0 * (1 - y) + p1 * x + p2 * y + p3 * (1 - x);
1335
1336   uv -= (1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3;
1337
1338   return uv;
1339 }
1340
1341 //=======================================================================
1342 //function : CalcUV2
1343 //purpose  : auxilary function for ComputeQuadPref
1344 //=======================================================================
1345
1346 static gp_UV CalcUV2(double x, double y,
1347                      FaceQuadStruct* quad,
1348                      const gp_UV& a0, const gp_UV& a1,
1349                      const gp_UV& a2, const gp_UV& a3)
1350 {
1351   const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0 );
1352   const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
1353   const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1 );
1354   const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
1355
1356   //double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
1357   //double y = y0 + x * (y1 - y0);
1358
1359   double param_b = uv_eb[0].normParam + x * (uv_eb.back().normParam - uv_eb[0].normParam);
1360   double param_t = uv_et[0].normParam + x * (uv_et.back().normParam - uv_et[0].normParam);
1361   double param_r = uv_er[0].normParam + y * (uv_er.back().normParam - uv_er[0].normParam);
1362   double param_l = uv_el[0].normParam + y * (uv_el.back().normParam - uv_el[0].normParam);
1363
1364   gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(param_b).XY();
1365   gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(param_r).XY();
1366   gp_UV p2 = quad->side[TOP_SIDE   ]->Value2d(param_t).XY();
1367   gp_UV p3 = quad->side[LEFT_SIDE  ]->Value2d(param_l).XY();
1368
1369   gp_UV uv = p0 * (1 - y) + p1 * x + p2 * y + p3 * (1 - x);
1370
1371   uv -= (1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3;
1372
1373   return uv;
1374 }
1375
1376
1377 //=======================================================================
1378 /*!
1379  * Create only quandrangle faces
1380  */
1381 //=======================================================================
1382
1383 bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh &        aMesh,
1384                                                 const TopoDS_Shape& aShape,
1385                                                 FaceQuadStruct*     quad)
1386 {
1387   // Auxilary key in order to keep old variant
1388   // of meshing after implementation new variant
1389   // for bug 0016220 from Mantis.
1390   bool OldVersion = false;
1391
1392   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
1393   const TopoDS_Face& F = TopoDS::Face(aShape);
1394   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
1395 //  const TopoDS_Wire& W = BRepTools::OuterWire(F);
1396   bool WisF = true;
1397 //   if(W.Orientation()==TopAbs_FORWARD)
1398 //     WisF = true;
1399   //if(WisF) cout<<"W is FORWARD"<<endl;
1400   //else cout<<"W is REVERSED"<<endl;
1401 //   bool FisF = (F.Orientation()==TopAbs_FORWARD);
1402 //   if(!FisF) WisF = !WisF;
1403 //  WisF = FisF;
1404   int i,j,geomFaceID = meshDS->ShapeToIndex( F );
1405
1406   int nb = quad->side[0]->NbPoints();
1407   int nr = quad->side[1]->NbPoints();
1408   int nt = quad->side[2]->NbPoints();
1409   int nl = quad->side[3]->NbPoints();
1410   int dh = abs(nb-nt);
1411   int dv = abs(nr-nl);
1412
1413   if( dh>=dv ) {
1414     if( nt>nb ) {
1415       // it is a base case => not shift quad but me be replacement is need
1416       ShiftQuad(quad,0,WisF);
1417     }
1418     else {
1419       // we have to shift quad on 2
1420       ShiftQuad(quad,2,WisF);
1421     }
1422   }
1423   else {
1424     if( nr>nl ) {
1425       // we have to shift quad on 1
1426       ShiftQuad(quad,1,WisF);
1427     }
1428     else {
1429       // we have to shift quad on 3
1430       ShiftQuad(quad,3,WisF);
1431     }
1432   }
1433
1434   nb = quad->side[0]->NbPoints();
1435   nr = quad->side[1]->NbPoints();
1436   nt = quad->side[2]->NbPoints();
1437   nl = quad->side[3]->NbPoints();
1438   dh = abs(nb-nt);
1439   dv = abs(nr-nl);
1440   int nbh  = Max(nb,nt);
1441   int nbv = Max(nr,nl);
1442   int addh = 0;
1443   int addv = 0;
1444
1445   // ----------- Old version ---------------
1446   // orientation of face and 3 main domain for future faces
1447   //       0   top    1
1448   //      1------------1
1449   //       |   |  |   |
1450   //       |   |  |   |
1451   //       | L |  | R |
1452   //  left |   |  |   | rigth
1453   //       |  /    \  |
1454   //       | /  C   \ |
1455   //       |/        \|
1456   //      0------------0
1457   //       0  bottom  1
1458
1459   // ----------- New version ---------------
1460   // orientation of face and 3 main domain for future faces
1461   //       0   top    1
1462   //      1------------1
1463   //       |  |____|  |
1464   //       |  /    \  |
1465   //       | /  C   \ |
1466   //  left |/________\| rigth
1467   //       |          |
1468   //       |          |
1469   //       |          |
1470   //      0------------0
1471   //       0  bottom  1
1472
1473   if(dh>dv) {
1474     addv = (dh-dv)/2;
1475     nbv = nbv + addv;
1476   }
1477   else { // dv>=dh
1478     addh = (dv-dh)/2;
1479     nbh = nbh + addh;
1480   }
1481
1482   const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0 );
1483   const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
1484   const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1 );
1485   const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
1486
1487   if ( uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl )
1488     return error( COMPERR_BAD_INPUT_MESH );
1489
1490   // arrays for normalized params
1491   //cout<<"Dump B:"<<endl;
1492   TColStd_SequenceOfReal npb, npr, npt, npl;
1493   for(i=0; i<nb; i++) {
1494     npb.Append(uv_eb[i].normParam);
1495     //cout<<"i="<<i<<" par="<<uv_eb[i].normParam<<" npar="<<uv_eb[i].normParam;
1496     //const SMDS_MeshNode* N = uv_eb[i].node;
1497     //cout<<" node("<<N->X()<<","<<N->Y()<<","<<N->Z()<<")"<<endl;
1498   }
1499   for(i=0; i<nr; i++) {
1500     npr.Append(uv_er[i].normParam);
1501   }
1502   for(i=0; i<nt; i++) {
1503     npt.Append(uv_et[i].normParam);
1504   }
1505   for(i=0; i<nl; i++) {
1506     npl.Append(uv_el[i].normParam);
1507   }
1508
1509   int dl,dr;
1510   if(OldVersion) {
1511     // add some params to right and left after the first param
1512     // insert to right
1513     dr = nbv - nr;
1514     double dpr = (npr.Value(2) - npr.Value(1))/(dr+1);
1515     for(i=1; i<=dr; i++) {
1516       npr.InsertAfter(1,npr.Value(2)-dpr);
1517     }
1518     // insert to left
1519     dl = nbv - nl;
1520     dpr = (npl.Value(2) - npl.Value(1))/(dl+1);
1521     for(i=1; i<=dl; i++) {
1522       npl.InsertAfter(1,npl.Value(2)-dpr);
1523     }
1524   }
1525   //cout<<"npb:";
1526   //for(i=1; i<=npb.Length(); i++) {
1527   //  cout<<" "<<npb.Value(i);
1528   //}
1529   //cout<<endl;
1530
1531   gp_XY a0( uv_eb.front().u, uv_eb.front().v );
1532   gp_XY a1( uv_eb.back().u,  uv_eb.back().v );
1533   gp_XY a2( uv_et.back().u,  uv_et.back().v );
1534   gp_XY a3( uv_et.front().u, uv_et.front().v );
1535   //cout<<" a0("<<a0.X()<<","<<a0.Y()<<")"<<" a1("<<a1.X()<<","<<a1.Y()<<")"
1536   //    <<" a2("<<a2.X()<<","<<a2.Y()<<")"<<" a3("<<a3.X()<<","<<a3.Y()<<")"<<endl;
1537
1538   int nnn = Min(nr,nl);
1539   // auxilary sequence of XY for creation nodes
1540   // in the bottom part of central domain
1541   // it's length must be == nbv-nnn-1
1542   TColgp_SequenceOfXY UVL;
1543   TColgp_SequenceOfXY UVR;
1544
1545   if(OldVersion) {
1546     // step1: create faces for left domain
1547     StdMeshers_Array2OfNode NodesL(1,dl+1,1,nl);
1548     // add left nodes
1549     for(j=1; j<=nl; j++)
1550       NodesL.SetValue(1,j,uv_el[j-1].node);
1551     if(dl>0) {
1552       // add top nodes
1553       for(i=1; i<=dl; i++)
1554         NodesL.SetValue(i+1,nl,uv_et[i].node);
1555       // create and add needed nodes
1556       TColgp_SequenceOfXY UVtmp;
1557       for(i=1; i<=dl; i++) {
1558         double x0 = npt.Value(i+1);
1559         double x1 = x0;
1560         // diagonal node
1561         double y0 = npl.Value(i+1);
1562         double y1 = npr.Value(i+1);
1563         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1564         gp_Pnt P = S->Value(UV.X(),UV.Y());
1565         SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1566         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1567         NodesL.SetValue(i+1,1,N);
1568         if(UVL.Length()<nbv-nnn-1) UVL.Append(UV);
1569         // internal nodes
1570         for(j=2; j<nl; j++) {
1571           double y0 = npl.Value(dl+j);
1572           double y1 = npr.Value(dl+j);
1573           gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1574           gp_Pnt P = S->Value(UV.X(),UV.Y());
1575           SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1576           meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1577           NodesL.SetValue(i+1,j,N);
1578           if( i==dl ) UVtmp.Append(UV);
1579         }
1580       }
1581       for(i=1; i<=UVtmp.Length() && UVL.Length()<nbv-nnn-1; i++) {
1582         UVL.Append(UVtmp.Value(i));
1583       }
1584       //cout<<"Dump NodesL:"<<endl;
1585       //for(i=1; i<=dl+1; i++) {
1586       //  cout<<"i="<<i;
1587       //  for(j=1; j<=nl; j++) {
1588       //    cout<<" ("<<NodesL.Value(i,j)->X()<<","<<NodesL.Value(i,j)->Y()<<","<<NodesL.Value(i,j)->Z()<<")";
1589       //  }
1590       //  cout<<endl;
1591       //}
1592       // create faces
1593       for(i=1; i<=dl; i++) {
1594         for(j=1; j<nl; j++) {
1595           if(WisF) {
1596             SMDS_MeshFace* F =
1597               myTool->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j),
1598                               NodesL.Value(i+1,j+1), NodesL.Value(i,j+1));
1599             if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1600           }
1601           else {
1602             SMDS_MeshFace* F =
1603               myTool->AddFace(NodesL.Value(i,j), NodesL.Value(i,j+1),
1604                               NodesL.Value(i+1,j+1), NodesL.Value(i+1,j));
1605             if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1606           }
1607         }
1608       }
1609     }
1610     else {
1611       // fill UVL using c2d
1612       for(i=1; i<npl.Length() && UVL.Length()<nbv-nnn-1; i++) {
1613         UVL.Append( gp_UV ( uv_el[i].u, uv_el[i].v ));
1614       }
1615     }
1616
1617     // step2: create faces for right domain
1618     StdMeshers_Array2OfNode NodesR(1,dr+1,1,nr);
1619     // add right nodes
1620     for(j=1; j<=nr; j++)
1621       NodesR.SetValue(1,j,uv_er[nr-j].node);
1622     if(dr>0) {
1623       // add top nodes
1624       for(i=1; i<=dr; i++)
1625         NodesR.SetValue(i+1,1,uv_et[nt-1-i].node);
1626       // create and add needed nodes
1627       TColgp_SequenceOfXY UVtmp;
1628       for(i=1; i<=dr; i++) {
1629         double x0 = npt.Value(nt-i);
1630         double x1 = x0;
1631         // diagonal node
1632         double y0 = npl.Value(i+1);
1633         double y1 = npr.Value(i+1);
1634         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1635         gp_Pnt P = S->Value(UV.X(),UV.Y());
1636         SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1637         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1638         NodesR.SetValue(i+1,nr,N);
1639         if(UVR.Length()<nbv-nnn-1) UVR.Append(UV);
1640         // internal nodes
1641         for(j=2; j<nr; j++) {
1642           double y0 = npl.Value(nbv-j+1);
1643           double y1 = npr.Value(nbv-j+1);
1644           gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1645           gp_Pnt P = S->Value(UV.X(),UV.Y());
1646           SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1647           meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1648           NodesR.SetValue(i+1,j,N);
1649           if( i==dr ) UVtmp.Prepend(UV);
1650         }
1651       }
1652       for(i=1; i<=UVtmp.Length() && UVR.Length()<nbv-nnn-1; i++) {
1653         UVR.Append(UVtmp.Value(i));
1654       }
1655       // create faces
1656       for(i=1; i<=dr; i++) {
1657         for(j=1; j<nr; j++) {
1658           if(WisF) {
1659             SMDS_MeshFace* F =
1660               myTool->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j),
1661                               NodesR.Value(i+1,j+1), NodesR.Value(i,j+1));
1662             if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1663           }
1664           else {
1665             SMDS_MeshFace* F =
1666               myTool->AddFace(NodesR.Value(i,j), NodesR.Value(i,j+1),
1667                               NodesR.Value(i+1,j+1), NodesR.Value(i+1,j));
1668             if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1669           }
1670         }
1671       }
1672     }
1673     else {
1674       // fill UVR using c2d
1675       for(i=1; i<npr.Length() && UVR.Length()<nbv-nnn-1; i++) {
1676         UVR.Append( gp_UV( uv_er[i].u, uv_er[i].v ));
1677       }
1678     }
1679
1680     // step3: create faces for central domain
1681     StdMeshers_Array2OfNode NodesC(1,nb,1,nbv);
1682     // add first string using NodesL
1683     for(i=1; i<=dl+1; i++)
1684       NodesC.SetValue(1,i,NodesL(i,1));
1685     for(i=2; i<=nl; i++)
1686       NodesC.SetValue(1,dl+i,NodesL(dl+1,i));
1687     // add last string using NodesR
1688     for(i=1; i<=dr+1; i++)
1689       NodesC.SetValue(nb,i,NodesR(i,nr));
1690     for(i=1; i<nr; i++)
1691       NodesC.SetValue(nb,dr+i+1,NodesR(dr+1,nr-i));
1692     // add top nodes (last columns)
1693     for(i=dl+2; i<nbh-dr; i++)
1694       NodesC.SetValue(i-dl,nbv,uv_et[i-1].node);
1695     // add bottom nodes (first columns)
1696     for(i=2; i<nb; i++)
1697       NodesC.SetValue(i,1,uv_eb[i-1].node);
1698
1699     // create and add needed nodes
1700     // add linear layers
1701     for(i=2; i<nb; i++) {
1702       double x0 = npt.Value(dl+i);
1703       double x1 = x0;
1704       for(j=1; j<nnn; j++) {
1705         double y0 = npl.Value(nbv-nnn+j);
1706         double y1 = npr.Value(nbv-nnn+j);
1707         gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1708         gp_Pnt P = S->Value(UV.X(),UV.Y());
1709         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1710         meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1711         NodesC.SetValue(i,nbv-nnn+j,N);
1712       }
1713     }
1714     // add diagonal layers
1715     //cout<<"UVL.Length()="<<UVL.Length()<<" UVR.Length()="<<UVR.Length()<<endl;
1716     //cout<<"Dump UVL:"<<endl;
1717     //for(i=1; i<=UVL.Length(); i++) {
1718     //  cout<<" ("<<UVL.Value(i).X()<<","<<UVL.Value(i).Y()<<")";
1719     //}
1720     //cout<<endl;
1721     for(i=1; i<nbv-nnn; i++) {
1722       double du = UVR.Value(i).X() - UVL.Value(i).X();
1723       double dv = UVR.Value(i).Y() - UVL.Value(i).Y();
1724       for(j=2; j<nb; j++) {
1725         double u = UVL.Value(i).X() + du*npb.Value(j);
1726         double v = UVL.Value(i).Y() + dv*npb.Value(j);
1727         gp_Pnt P = S->Value(u,v);
1728         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1729         meshDS->SetNodeOnFace(N, geomFaceID, u, v);
1730         NodesC.SetValue(j,i+1,N);
1731       }
1732     }
1733     // create faces
1734     for(i=1; i<nb; i++) {
1735       for(j=1; j<nbv; j++) {
1736         if(WisF) {
1737           SMDS_MeshFace* F =
1738             myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
1739                             NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
1740           if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1741         }
1742         else {
1743           SMDS_MeshFace* F =
1744             myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1),
1745                             NodesC.Value(i+1,j+1), NodesC.Value(i+1,j));
1746           if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1747         }
1748       }
1749     }
1750   }
1751
1752   else { // New version (!OldVersion)
1753     // step1: create faces for bottom rectangle domain
1754     StdMeshers_Array2OfNode NodesBRD(1,nb,1,nnn-1);
1755     // fill UVL and UVR using c2d
1756     for(j=0; j<nb; j++) {
1757       NodesBRD.SetValue(j+1,1,uv_eb[j].node);
1758     }
1759     for(i=1; i<nnn-1; i++) {
1760       NodesBRD.SetValue(1,i+1,uv_el[i].node);
1761       NodesBRD.SetValue(nb,i+1,uv_er[i].node);
1762       double du = uv_er[i].u - uv_el[i].u;
1763       double dv = uv_er[i].v - uv_el[i].v;
1764       for(j=2; j<nb; j++) {
1765         double u = uv_el[i].u + du*npb.Value(j);
1766         double v = uv_el[i].v + dv*npb.Value(j);
1767         gp_Pnt P = S->Value(u,v);
1768         SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1769         meshDS->SetNodeOnFace(N, geomFaceID, u, v);
1770         NodesBRD.SetValue(j,i+1,N);
1771
1772       }
1773     }
1774     int nbf=0;
1775     for(j=1; j<nnn-1; j++) {
1776       for(i=1; i<nb; i++) {
1777         nbf++;
1778         if(WisF) {
1779           SMDS_MeshFace* F =
1780             myTool->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i+1,j),
1781                             NodesBRD.Value(i+1,j+1), NodesBRD.Value(i,j+1));
1782           if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1783         }
1784         else {
1785           SMDS_MeshFace* F =
1786             myTool->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i,j+1),
1787                             NodesBRD.Value(i+1,j+1), NodesBRD.Value(i+1,j));
1788           if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1789         }
1790       }
1791     }
1792     int drl = abs(nr-nl);
1793     // create faces for region C
1794     StdMeshers_Array2OfNode NodesC(1,nb,1,drl+1+addv);
1795     // add nodes from previous region
1796     for(j=1; j<=nb; j++) {
1797       NodesC.SetValue(j,1,NodesBRD.Value(j,nnn-1));
1798     }
1799     if( (drl+addv) > 0 ) {
1800       int n1,n2;
1801       if(nr>nl) {
1802         n1 = 1;
1803         n2 = drl + 1;
1804         TColgp_SequenceOfXY UVtmp;
1805         double drparam = npr.Value(nr) - npr.Value(nnn-1);
1806         double dlparam = npl.Value(nnn) - npl.Value(nnn-1);
1807         double y0,y1;
1808         for(i=1; i<=drl; i++) {
1809           // add existed nodes from right edge
1810           NodesC.SetValue(nb,i+1,uv_er[nnn+i-2].node);
1811           //double dtparam = npt.Value(i+1);
1812           y1 = npr.Value(nnn+i-1); // param on right edge
1813           double dpar = (y1 - npr.Value(nnn-1))/drparam;
1814           y0 = npl.Value(nnn-1) + dpar*dlparam; // param on left edge
1815           double dy = y1 - y0;
1816           for(j=1; j<nb; j++) {
1817             double x = npt.Value(i+1) + npb.Value(j)*(1-npt.Value(i+1));
1818             double y = y0 + dy*x;
1819             gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1820             gp_Pnt P = S->Value(UV.X(),UV.Y());
1821             SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1822             meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1823             NodesC.SetValue(j,i+1,N);
1824           }
1825         }
1826         double dy0 = (1-y0)/(addv+1);
1827         double dy1 = (1-y1)/(addv+1);
1828         for(i=1; i<=addv; i++) {
1829           double yy0 = y0 + dy0*i;
1830           double yy1 = y1 + dy1*i;
1831           double dyy = yy1 - yy0;
1832           for(j=1; j<=nb; j++) {
1833             double x = npt.Value(i+1+drl) +
1834               npb.Value(j) * ( npt.Value(nt-i) - npt.Value(i+1+drl) );
1835             double y = yy0 + dyy*x;
1836             gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1837             gp_Pnt P = S->Value(UV.X(),UV.Y());
1838             SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1839             meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1840             NodesC.SetValue(j,i+drl+1,N);
1841           }
1842         }
1843       }
1844       else { // nr<nl
1845         n2 = 1;
1846         n1 = drl + 1;
1847         TColgp_SequenceOfXY UVtmp;
1848         double dlparam = npl.Value(nl) - npl.Value(nnn-1);
1849         double drparam = npr.Value(nnn) - npr.Value(nnn-1);
1850         double y0 = npl.Value(nnn-1);
1851         double y1 = npr.Value(nnn-1);
1852         for(i=1; i<=drl; i++) {
1853           // add existed nodes from right edge
1854           NodesC.SetValue(1,i+1,uv_el[nnn+i-2].node);
1855           y0 = npl.Value(nnn+i-1); // param on left edge
1856           double dpar = (y0 - npl.Value(nnn-1))/dlparam;
1857           y1 = npr.Value(nnn-1) + dpar*drparam; // param on right edge
1858           double dy = y1 - y0;
1859           for(j=2; j<=nb; j++) {
1860             double x = npb.Value(j)*npt.Value(nt-i);
1861             double y = y0 + dy*x;
1862             gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1863             gp_Pnt P = S->Value(UV.X(),UV.Y());
1864             SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1865             meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1866             NodesC.SetValue(j,i+1,N);
1867           }
1868         }
1869         double dy0 = (1-y0)/(addv+1);
1870         double dy1 = (1-y1)/(addv+1);
1871         for(i=1; i<=addv; i++) {
1872           double yy0 = y0 + dy0*i;
1873           double yy1 = y1 + dy1*i;
1874           double dyy = yy1 - yy0;
1875           for(j=1; j<=nb; j++) {
1876             double x = npt.Value(i+1) +
1877               npb.Value(j) * ( npt.Value(nt-i-drl) - npt.Value(i+1) );
1878             double y = yy0 + dyy*x;
1879             gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
1880             gp_Pnt P = S->Value(UV.X(),UV.Y());
1881             SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1882             meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1883             NodesC.SetValue(j,i+drl+1,N);
1884           }
1885         }
1886       }
1887       // create faces
1888       for(j=1; j<=drl+addv; j++) {
1889         for(i=1; i<nb; i++) {
1890           nbf++;
1891           if(WisF) {
1892             SMDS_MeshFace* F =
1893               myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
1894                               NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
1895             if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1896           }
1897           else {
1898             SMDS_MeshFace* F =
1899               myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1),
1900                               NodesC.Value(i+1,j+1), NodesC.Value(i+1,j));
1901             if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1902           }
1903         }
1904       } // end nr<nl
1905
1906       StdMeshers_Array2OfNode NodesLast(1,nt,1,2);
1907       for(i=1; i<=nt; i++) {
1908         NodesLast.SetValue(i,2,uv_et[i-1].node);
1909       }
1910       int nnn=0;
1911       for(i=n1; i<drl+addv+1; i++) {
1912         nnn++;
1913         NodesLast.SetValue(nnn,1,NodesC.Value(1,i));
1914       }
1915       for(i=1; i<=nb; i++) {
1916         nnn++;
1917         NodesLast.SetValue(nnn,1,NodesC.Value(i,drl+addv+1));
1918       }
1919       for(i=drl+addv; i>=n2; i--) {
1920         nnn++;
1921         NodesLast.SetValue(nnn,1,NodesC.Value(nb,i));
1922       }
1923       for(i=1; i<nt; i++) {
1924         nbf++;
1925         if(WisF) {
1926           SMDS_MeshFace* F =
1927             myTool->AddFace(NodesLast.Value(i,1), NodesLast.Value(i+1,1),
1928                             NodesLast.Value(i+1,2), NodesLast.Value(i,2));
1929           if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1930         }
1931         else {
1932           SMDS_MeshFace* F =
1933             myTool->AddFace(NodesLast.Value(i,1), NodesLast.Value(i,2),
1934                             NodesLast.Value(i+1,2), NodesLast.Value(i+1,2));
1935           if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1936         }
1937       }
1938     } // if( (drl+addv) > 0 )
1939
1940   } // end new version implementation
1941
1942   bool isOk = true;
1943   return isOk;
1944 }
1945
1946
1947 //=======================================================================
1948 /*!
1949  * Evaluate only quandrangle faces
1950  */
1951 //=======================================================================
1952
1953 bool StdMeshers_Quadrangle_2D::EvaluateQuadPref(SMESH_Mesh &        aMesh,
1954                                                 const TopoDS_Shape& aShape,
1955                                                 std::vector<int>& aNbNodes,
1956                                                 MapShapeNbElems& aResMap,
1957                                                 bool IsQuadratic)
1958 {
1959   // Auxilary key in order to keep old variant
1960   // of meshing after implementation new variant
1961   // for bug 0016220 from Mantis.
1962   bool OldVersion = false;
1963
1964   const TopoDS_Face& F = TopoDS::Face(aShape);
1965   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
1966
1967   int nb = aNbNodes[0];
1968   int nr = aNbNodes[1];
1969   int nt = aNbNodes[2];
1970   int nl = aNbNodes[3];
1971   int dh = abs(nb-nt);
1972   int dv = abs(nr-nl);
1973
1974   if( dh>=dv ) {
1975     if( nt>nb ) {
1976       // it is a base case => not shift
1977     }
1978     else {
1979       // we have to shift on 2
1980       nb = aNbNodes[2];
1981       nr = aNbNodes[3];
1982       nt = aNbNodes[0];
1983       nl = aNbNodes[1];
1984     }
1985   }
1986   else {
1987     if( nr>nl ) {
1988       // we have to shift quad on 1
1989       nb = aNbNodes[3];
1990       nr = aNbNodes[0];
1991       nt = aNbNodes[1];
1992       nl = aNbNodes[2];
1993     }
1994     else {
1995       // we have to shift quad on 3
1996       nb = aNbNodes[1];
1997       nr = aNbNodes[2];
1998       nt = aNbNodes[3];
1999       nl = aNbNodes[0];
2000     }
2001   }
2002
2003   dh = abs(nb-nt);
2004   dv = abs(nr-nl);
2005   int nbh  = Max(nb,nt);
2006   int nbv = Max(nr,nl);
2007   int addh = 0;
2008   int addv = 0;
2009
2010   if(dh>dv) {
2011     addv = (dh-dv)/2;
2012     nbv = nbv + addv;
2013   }
2014   else { // dv>=dh
2015     addh = (dv-dh)/2;
2016     nbh = nbh + addh;
2017   }
2018
2019   int dl,dr;
2020   if(OldVersion) {
2021     // add some params to right and left after the first param
2022     // insert to right
2023     dr = nbv - nr;
2024     // insert to left
2025     dl = nbv - nl;
2026   }
2027
2028   int nnn = Min(nr,nl);
2029
2030   int nbNodes = 0;
2031   int nbFaces = 0;
2032   if(OldVersion) {
2033     // step1: create faces for left domain
2034     if(dl>0) {
2035       nbNodes += dl*(nl-1);
2036       nbFaces += dl*(nl-1);
2037     }
2038     // step2: create faces for right domain
2039     if(dr>0) {
2040       nbNodes += dr*(nr-1);
2041       nbFaces += dr*(nr-1);
2042     }
2043     // step3: create faces for central domain
2044     nbNodes += (nb-2)*(nnn-1) + (nbv-nnn-1)*(nb-2);
2045     nbFaces += (nb-1)*(nbv-1);
2046   }
2047   else { // New version (!OldVersion)
2048     nbNodes += (nnn-2)*(nb-2);
2049     nbFaces += (nnn-2)*(nb-1);
2050     int drl = abs(nr-nl);
2051     nbNodes += drl*(nb-1) + addv*nb;
2052     nbFaces += (drl+addv)*(nb-1) + (nt-1);
2053   } // end new version implementation
2054
2055   std::vector<int> aVec(SMDSEntity_Last);
2056   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i] = 0;
2057   if(IsQuadratic) {
2058     aVec[SMDSEntity_Quad_Quadrangle] = nbFaces;
2059     aVec[SMDSEntity_Node] = nbNodes + nbFaces*4;
2060     if( aNbNodes.size()==5 ) {
2061       aVec[SMDSEntity_Quad_Triangle] = aNbNodes[3] - 1;
2062       aVec[SMDSEntity_Quad_Quadrangle] = nbFaces - aNbNodes[3] + 1;
2063     }
2064   }
2065   else {
2066     aVec[SMDSEntity_Node] = nbNodes;
2067     aVec[SMDSEntity_Quadrangle] = nbFaces;
2068     if( aNbNodes.size()==5 ) {
2069       aVec[SMDSEntity_Triangle] = aNbNodes[3] - 1;
2070       aVec[SMDSEntity_Quadrangle] = nbFaces - aNbNodes[3] + 1;
2071     }
2072   }
2073   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
2074   aResMap.insert(std::make_pair(sm,aVec));
2075
2076   return true;
2077 }
2078
2079
2080 //=============================================================================
2081 /*! Split quadrangle in to 2 triangles by smallest diagonal
2082  *
2083  */
2084 //=============================================================================
2085 void StdMeshers_Quadrangle_2D::SplitQuad(SMESHDS_Mesh *theMeshDS,
2086                                     int theFaceID,
2087                                     const SMDS_MeshNode* theNode1,
2088                                     const SMDS_MeshNode* theNode2,
2089                                     const SMDS_MeshNode* theNode3,
2090                                     const SMDS_MeshNode* theNode4)
2091 {
2092   gp_Pnt a(theNode1->X(),theNode1->Y(),theNode1->Z());
2093   gp_Pnt b(theNode2->X(),theNode2->Y(),theNode2->Z());
2094   gp_Pnt c(theNode3->X(),theNode3->Y(),theNode3->Z());
2095   gp_Pnt d(theNode4->X(),theNode4->Y(),theNode4->Z());
2096   SMDS_MeshFace* face;
2097   if(a.Distance(c) > b.Distance(d)){
2098     face = myTool->AddFace(theNode2, theNode4 , theNode1);
2099     if(face) theMeshDS->SetMeshElementOnShape(face, theFaceID );
2100     face = myTool->AddFace(theNode2, theNode3, theNode4);
2101     if(face) theMeshDS->SetMeshElementOnShape(face, theFaceID );
2102
2103   }
2104   else{
2105     face = myTool->AddFace(theNode1, theNode2 ,theNode3);
2106     if(face) theMeshDS->SetMeshElementOnShape(face, theFaceID );
2107     face = myTool->AddFace(theNode1, theNode3, theNode4);
2108     if(face) theMeshDS->SetMeshElementOnShape(face, theFaceID );
2109   }
2110 }
2111
2112