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