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