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