src/StdMeshers/StdMeshers_Quadrangle_2D.cxx

   1 //  SMESH SMESH : implementaion of SMESH idl descriptions
   2 //
   3 //  Copyright (C) 2003  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.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  21 //
  22 //
  23 //
  24 //  File   : StdMeshers_Quadrangle_2D.cxx
  25 //           Moved here from SMESH_Quadrangle_2D.cxx
  26 //  Author : Paul RASCLE, EDF
  27 //  Module : SMESH
  28 //  $Header$
  29
  30 using namespace std;
  31 #include "StdMeshers_Quadrangle_2D.hxx"
  32 #include "SMESH_Gen.hxx"
  33 #include "SMESH_Mesh.hxx"
  34 #include "SMESH_subMesh.hxx"
  35
  36 #include "SMDS_MeshElement.hxx"
  37 #include "SMDS_MeshNode.hxx"
  38 #include "SMDS_EdgePosition.hxx"
  39 #include "SMDS_FacePosition.hxx"
  40
  41 #include <BRep_Tool.hxx>
  42 #include <BRepTools.hxx>
  43 #include <BRepTools_WireExplorer.hxx>
  44
  45 #include <Geom_Surface.hxx>
  46 #include <Geom_Curve.hxx>
  47 #include <Geom2d_Curve.hxx>
  48 #include <GeomAdaptor_Curve.hxx>
  49 #include <GCPnts_UniformAbscissa.hxx>
  50
  51 #include <Precision.hxx>
  52 #include <gp_Pnt2d.hxx>
  53 #include <TColStd_ListIteratorOfListOfInteger.hxx>
  54
  55 #include "utilities.h"
  56 #include "Utils_ExceptHandlers.hxx"
  57
  58
  59 //=============================================================================
  60 /*!
  61  *
  62  */
  63 //=============================================================================
  64
  65 StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, SMESH_Gen* gen)
  66      : SMESH_2D_Algo(hypId, studyId, gen)
  67 {
  68   MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D");
  69   _name = "Quadrangle_2D";
  70   //  _shapeType = TopAbs_FACE;
  71   _shapeType = (1 << TopAbs_FACE);
  72 }
  73
  74 //=============================================================================
  75 /*!
  76  *
  77  */
  78 //=============================================================================
  79
  80 StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D()
  81 {
  82   MESSAGE("StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D");
  83 }
  84
  85 //=============================================================================
  86 /*!
  87  *
  88  */
  89 //=============================================================================
  90
  91 bool StdMeshers_Quadrangle_2D::CheckHypothesis
  92                          (SMESH_Mesh& aMesh,
  93                           const TopoDS_Shape& aShape,
  94                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
  95 {
  96   //MESSAGE("StdMeshers_Quadrangle_2D::CheckHypothesis");
  97
  98   bool isOk = true;
  99   aStatus = SMESH_Hypothesis::HYP_OK;
 100
 101   // nothing to check
 102
 103   return isOk;
 104 }
 105
 106 //=============================================================================
 107 /*!
 108  *
 109  */
 110 //=============================================================================
 111
 112 bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh,
 113                                         const TopoDS_Shape& aShape) throw (SALOME_Exception)
 114 {
 115   Unexpect aCatch(SalomeException);
 116   //MESSAGE("StdMeshers_Quadrangle_2D::Compute");
 117   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 118   aMesh.GetSubMesh(aShape);
 119
 120   FaceQuadStruct *quad = CheckAnd2Dcompute(aMesh, aShape);
 121   if (!quad)
 122     return false;
 123
 124   // --- compute 3D values on points, store points & quadrangles
 125
 126   int nbdown  = quad->nbPts[0];
 127   int nbup    = quad->nbPts[2];
 128 //  bool isDownOut = (nbdown > nbup);
 129 //  bool isUpOut   = (nbdown < nbup);
 130
 131   int nbright = quad->nbPts[1];
 132   int nbleft  = quad->nbPts[3];
 133 //  bool isRightOut = (nbright > nbleft);
 134 //  bool isLeftOut  = (nbright < nbleft);
 135
 136   int nbhoriz  = Min(nbdown, nbup);
 137   int nbvertic = Min(nbright, nbleft);
 138
 139   int nbVertices = nbhoriz * nbvertic;
 140   int nbQuad = (nbhoriz - 1) * (nbvertic - 1);
 141   //SCRUTE(nbVertices);
 142   //SCRUTE(nbQuad);
 143
 144   //   const TopoDS_Face& FF = TopoDS::Face(aShape);
 145   //   bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
 146   //   TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
 147   const TopoDS_Face& F = TopoDS::Face(aShape);
 148   bool faceIsForward = (F.Orientation() == TopAbs_FORWARD);
 149   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
 150
 151   // internal mesh nodes
 152   int i, j;
 153   for (i = 1; i < nbhoriz - 1; i++) {
 154     for (j = 1; j < nbvertic - 1; j++) {
 155       int ij = j * nbhoriz + i;
 156       double u = quad->uv_grid[ij].u;
 157       double v = quad->uv_grid[ij].v;
 158       gp_Pnt P = S->Value(u, v);
 159       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 160       meshDS->SetNodeOnFace(node, F);
 161       quad->uv_grid[ij].node = node;
 162       SMDS_FacePosition* fpos =
 163         dynamic_cast<SMDS_FacePosition*>(node->GetPosition().get());
 164       fpos->SetUParameter(u);
 165       fpos->SetVParameter(v);
 166     }
 167   }
 168
 169   // mesh faces
 170
 171   //             [2]
 172   //      --.--.--.--.--.--  nbvertic
 173   //     |                 | ^
 174   //     |                 | ^
 175   // [3] |                 | ^ j  [1]
 176   //     |                 | ^
 177   //     |                 | ^
 178   //      ---.----.----.---  0
 179   //     0 > > > > > > > > nbhoriz
 180   //              i
 181   //             [0]
 182
 183   i = 0;
 184   int ilow = 0;
 185   int iup = nbhoriz - 1;
 186   if (quad->isEdgeOut[3]) { ilow++; } else { if (quad->isEdgeOut[1]) iup--; }
 187
 188   int jlow = 0;
 189   int jup = nbvertic - 1;
 190   if (quad->isEdgeOut[0]) { jlow++; } else { if (quad->isEdgeOut[2]) jup--; }
 191
 192   // regular quadrangles
 193   //   bool isQuadForward = ( faceIsForward == quad->isEdgeForward[0]);
 194   for (i = ilow; i < iup; i++) {
 195     for (j = jlow; j < jup; j++) {
 196       const SMDS_MeshNode *a, *b, *c, *d;
 197       a = quad->uv_grid[j * nbhoriz + i].node;
 198       b = quad->uv_grid[j * nbhoriz + i + 1].node;
 199       c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node;
 200       d = quad->uv_grid[(j + 1) * nbhoriz + i].node;
 201       //  if (isQuadForward) faceId = meshDS->AddFace(a,b,c,d);
 202       //  else faceId = meshDS->AddFace(a,d,c,b);
 203       SMDS_MeshFace * face = meshDS->AddFace(a, b, c, d);
 204       meshDS->SetMeshElementOnShape(face, F);
 205     }
 206   }
 207
 208   UVPtStruct *uv_e0 = quad->uv_edges[0];
 209   UVPtStruct *uv_e1 = quad->uv_edges[1];
 210   UVPtStruct *uv_e2 = quad->uv_edges[2];
 211   UVPtStruct *uv_e3 = quad->uv_edges[3];
 212
 213   double eps = Precision::Confusion();
 214
 215   // Boundary quadrangles
 216
 217   if (quad->isEdgeOut[0]) {
 218     // Down edge is out
 219     //
 220     // |___|___|___|___|___|___|
 221     // |   |   |   |   |   |   |
 222     // |___|___|___|___|___|___|
 223     // |   |   |   |   |   |   |
 224     // |___|___|___|___|___|___| __ first row of the regular grid
 225     // .  .  .  .  .  .  .  .  . __ down edge nodes
 226     //
 227     // >->->->->->->->->->->->-> -- direction of processing
 228
 229     int g = 0; // number of last processed node in the regular grid
 230
 231     // number of last node of the down edge to be processed
 232     int stop = nbdown - 1;
 233     // if right edge is out, we will stop at a node, previous to the last one
 234     if (quad->isEdgeOut[1]) stop--;
 235
 236     // for each node of the down edge find nearest node
 237     // in the first row of the regular grid and link them
 238     for (i = 0; i < stop; i++) {
 239       const SMDS_MeshNode *a, *b, *c, *d;
 240       a = uv_e0[i].node;
 241       b = uv_e0[i + 1].node;
 242       gp_Pnt pb (b->X(), b->Y(), b->Z());
 243
 244       // find node c in the regular grid, which will be linked with node b
 245       int near = g;
 246       if (i == stop - 1) {
 247         // right bound reached, link with the rightmost node
 248         near = iup;
 249         c = quad->uv_grid[nbhoriz + iup].node;
 250       } else {
 251         // find in the grid node c, nearest to the b
 252         double mind = RealLast();
 253         for (int k = g; k <= iup; k++) {
 254
 255           const SMDS_MeshNode *nk;
 256           if (k < ilow) // this can be, if left edge is out
 257             nk = uv_e3[1].node; // get node from the left edge
 258           else
 259             nk = quad->uv_grid[nbhoriz + k].node; // get one of middle nodes
 260
 261           gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 262           double dist = pb.Distance(pnk);
 263           if (dist < mind - eps) {
 264             c = nk;
 265             near = k;
 266             mind = dist;
 267           } else {
 268             break;
 269           }
 270         }
 271       }
 272
 273       if (near == g) { // make triangle
 274         SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
 275         meshDS->SetMeshElementOnShape(face, F);
 276       } else { // make quadrangle
 277         if (near - 1 < ilow)
 278           d = uv_e3[1].node;
 279         else
 280           d = quad->uv_grid[nbhoriz + near - 1].node;
 281         SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 282         meshDS->SetMeshElementOnShape(face, F);
 283
 284         // if node d is not at position g - make additional triangles
 285         if (near - 1 > g) {
 286           for (int k = near - 1; k > g; k--) {
 287             c = quad->uv_grid[nbhoriz + k].node;
 288             if (k - 1 < ilow)
 289               d = uv_e3[1].node;
 290             else
 291               d = quad->uv_grid[nbhoriz + k - 1].node;
 292             SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
 293             meshDS->SetMeshElementOnShape(face, F);
 294           }
 295         }
 296         g = near;
 297       }
 298     }
 299   } else {
 300     if (quad->isEdgeOut[2]) {
 301       // Up edge is out
 302       //
 303       // <-<-<-<-<-<-<-<-<-<-<-<-< -- direction of processing
 304       //
 305       // .  .  .  .  .  .  .  .  . __ up edge nodes
 306       //  ___ ___ ___ ___ ___ ___  __ first row of the regular grid
 307       // |   |   |   |   |   |   |
 308       // |___|___|___|___|___|___|
 309       // |   |   |   |   |   |   |
 310       // |___|___|___|___|___|___|
 311       // |   |   |   |   |   |   |
 312
 313       int g = nbhoriz - 1; // last processed node in the regular grid
 314
 315       int stop = 0;
 316       // if left edge is out, we will stop at a second node
 317       if (quad->isEdgeOut[3]) stop++;
 318
 319       // for each node of the up edge find nearest node
 320       // in the first row of the regular grid and link them
 321       for (i = nbup - 1; i > stop; i--) {
 322         const SMDS_MeshNode *a, *b, *c, *d;
 323         a = uv_e2[i].node;
 324         b = uv_e2[i - 1].node;
 325         gp_Pnt pb (b->X(), b->Y(), b->Z());
 326
 327         // find node c in the grid, which will be linked with node b
 328         int near = g;
 329         if (i == stop + 1) { // left bound reached, link with the leftmost node
 330           c = quad->uv_grid[nbhoriz*(nbvertic - 2) + ilow].node;
 331           near = ilow;
 332         } else {
 333           // find node c in the grid, nearest to the b
 334           double mind = RealLast();
 335           for (int k = g; k >= ilow; k--) {
 336             const SMDS_MeshNode *nk;
 337             if (k > iup)
 338               nk = uv_e1[nbright - 2].node;
 339             else
 340               nk = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
 341             gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 342             double dist = pb.Distance(pnk);
 343             if (dist < mind - eps) {
 344               c = nk;
 345               near = k;
 346               mind = dist;
 347             } else {
 348               break;
 349             }
 350           }
 351         }
 352
 353         if (near == g) { // make triangle
 354           SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
 355           meshDS->SetMeshElementOnShape(face, F);
 356         } else { // make quadrangle
 357           if (near + 1 > iup)
 358             d = uv_e1[nbright - 2].node;
 359           else
 360             d = quad->uv_grid[nbhoriz*(nbvertic - 2) + near + 1].node;
 361           SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 362           meshDS->SetMeshElementOnShape(face, F);
 363
 364           if (near + 1 < g) { // if d not is at g - make additional triangles
 365             for (int k = near + 1; k < g; k++) {
 366               c = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
 367               if (k + 1 > iup)
 368                 d = uv_e1[nbright - 2].node;
 369               else
 370                 d = quad->uv_grid[nbhoriz*(nbvertic - 2) + k + 1].node;
 371               SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
 372               meshDS->SetMeshElementOnShape(face, F);
 373             }
 374           }
 375           g = near;
 376         }
 377       }
 378     }
 379   }
 380
 381   // right or left boundary quadrangles
 382   if (quad->isEdgeOut[1]) {
 383 //    MESSAGE("right edge is out");
 384     int g = 0; // last processed node in the grid
 385     int stop = nbright - 1;
 386     if (quad->isEdgeOut[2]) stop--;
 387     for (i = 0; i < stop; i++) {
 388       const SMDS_MeshNode *a, *b, *c, *d;
 389       a = uv_e1[i].node;
 390       b = uv_e1[i + 1].node;
 391       gp_Pnt pb (b->X(), b->Y(), b->Z());
 392
 393       // find node c in the grid, nearest to the b
 394       int near = g;
 395       if (i == stop - 1) { // up bondary reached
 396         c = quad->uv_grid[nbhoriz*(jup + 1) - 2].node;
 397         near = jup;
 398       } else {
 399         double mind = RealLast();
 400         for (int k = g; k <= jup; k++) {
 401           const SMDS_MeshNode *nk;
 402           if (k < jlow)
 403             nk = uv_e0[nbdown - 2].node;
 404           else
 405             nk = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
 406           gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 407           double dist = pb.Distance(pnk);
 408           if (dist < mind - eps) {
 409             c = nk;
 410             near = k;
 411             mind = dist;
 412           } else {
 413             break;
 414           }
 415         }
 416       }
 417
 418       if (near == g) { // make triangle
 419         SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
 420         meshDS->SetMeshElementOnShape(face, F);
 421       } else { // make quadrangle
 422         if (near - 1 < jlow)
 423           d = uv_e0[nbdown - 2].node;
 424         else
 425           d = quad->uv_grid[nbhoriz*near - 2].node;
 426         SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 427         meshDS->SetMeshElementOnShape(face, F);
 428
 429         if (near - 1 > g) { // if d not is at g - make additional triangles
 430           for (int k = near - 1; k > g; k--) {
 431             c = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
 432             if (k - 1 < jlow)
 433               d = uv_e0[nbdown - 2].node;
 434             else
 435               d = quad->uv_grid[nbhoriz*k - 2].node;
 436             SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
 437             meshDS->SetMeshElementOnShape(face, F);
 438           }
 439         }
 440         g = near;
 441       }
 442     }
 443   } else {
 444     if (quad->isEdgeOut[3]) {
 445 //      MESSAGE("left edge is out");
 446       int g = nbvertic - 1; // last processed node in the grid
 447       int stop = 0;
 448       if (quad->isEdgeOut[0]) stop++;
 449       for (i = nbleft - 1; i > stop; i--) {
 450         const SMDS_MeshNode *a, *b, *c, *d;
 451         a = uv_e3[i].node;
 452         b = uv_e3[i - 1].node;
 453         gp_Pnt pb (b->X(), b->Y(), b->Z());
 454
 455         // find node c in the grid, nearest to the b
 456         int near = g;
 457         if (i == stop + 1) { // down bondary reached
 458           c = quad->uv_grid[nbhoriz*jlow + 1].node;
 459           near = jlow;
 460         } else {
 461           double mind = RealLast();
 462           for (int k = g; k >= jlow; k--) {
 463             const SMDS_MeshNode *nk;
 464             if (k > jup)
 465               nk = uv_e2[1].node;
 466             else
 467               nk = quad->uv_grid[nbhoriz*k + 1].node;
 468             gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
 469             double dist = pb.Distance(pnk);
 470             if (dist < mind - eps) {
 471               c = nk;
 472               near = k;
 473               mind = dist;
 474             } else {
 475               break;
 476             }
 477           }
 478         }
 479
 480         if (near == g) { // make triangle
 481           SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
 482           meshDS->SetMeshElementOnShape(face, F);
 483         } else { // make quadrangle
 484           if (near + 1 > jup)
 485             d = uv_e2[1].node;
 486           else
 487             d = quad->uv_grid[nbhoriz*(near + 1) + 1].node;
 488           SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
 489           meshDS->SetMeshElementOnShape(face, F);
 490
 491           if (near + 1 < g) { // if d not is at g - make additional triangles
 492             for (int k = near + 1; k < g; k++) {
 493               c = quad->uv_grid[nbhoriz*k + 1].node;
 494               if (k + 1 > jup)
 495                 d = uv_e2[1].node;
 496               else
 497                 d = quad->uv_grid[nbhoriz*(k + 1) + 1].node;
 498               SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
 499               meshDS->SetMeshElementOnShape(face, F);
 500             }
 501           }
 502           g = near;
 503         }
 504       }
 505     }
 506   }
 507
 508   QuadDelete(quad);
 509   bool isOk = true;
 510   return isOk;
 511 }
 512
 513 //=============================================================================
 514 /*!
 515  *
 516  */
 517 //=============================================================================
 518
 519 FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute
 520   (SMESH_Mesh & aMesh, const TopoDS_Shape & aShape) throw(SALOME_Exception)
 521 {
 522   Unexpect aCatch(SalomeException);
 523 //  MESSAGE("StdMeshers_Quadrangle_2D::CheckAnd2Dcompute");
 524
 525   SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape);
 526
 527   //   const TopoDS_Face& FF = TopoDS::Face(aShape);
 528   //   bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
 529   //   TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
 530   const TopoDS_Face & F = TopoDS::Face(aShape);
 531   bool faceIsForward = (F.Orientation() == TopAbs_FORWARD);
 532
 533   // verify 1 wire only, with 4 edges
 534
 535   if (NumberOfWires(F) != 1)
 536   {
 537     MESSAGE("only 1 wire by face (quadrangles)");
 538     return 0;
 539     //throw SALOME_Exception(LOCALIZED("only 1 wire by face (quadrangles)"));
 540   }
 541   //   const TopoDS_Wire WW = BRepTools::OuterWire(F);
 542   //   TopoDS_Wire W = TopoDS::Wire(WW.Oriented(TopAbs_FORWARD));
 543   const TopoDS_Wire& W = BRepTools::OuterWire(F);
 544   BRepTools_WireExplorer wexp (W, F);
 545
 546   FaceQuadStruct *quad = new FaceQuadStruct;
 547   for (int i = 0; i < 4; i++)
 548     quad->uv_edges[i] = 0;
 549   quad->uv_grid = 0;
 550
 551   int nbEdges = 0;
 552   for (wexp.Init(W, F); wexp.More(); wexp.Next())
 553   {
 554     //       const TopoDS_Edge& EE = wexp.Current();
 555     //       TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
 556     const TopoDS_Edge& E = wexp.Current();
 557     int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
 558     if (nbEdges < 4)
 559     {
 560       quad->edge[nbEdges] = E;
 561       quad->nbPts[nbEdges] = nb + 2; // internal points + 2 extrema
 562     }
 563     nbEdges++;
 564   }
 565
 566   if (nbEdges != 4)
 567   {
 568     MESSAGE("face must have 4 edges /quadrangles");
 569     QuadDelete(quad);
 570     return 0;
 571     //throw SALOME_Exception(LOCALIZED("face must have 4 edges /quadrangles"));
 572   }
 573
 574 //  if (quad->nbPts[0] != quad->nbPts[2]) {
 575 //    MESSAGE("different point number-opposed edge");
 576 //    QuadDelete(quad);
 577 //    return 0;
 578 //    //throw SALOME_Exception(LOCALIZED("different point number-opposed edge"));
 579 //  }
 580 //
 581 //  if (quad->nbPts[1] != quad->nbPts[3]) {
 582 //    MESSAGE("different point number-opposed edge");
 583 //    QuadDelete(quad);
 584 //    return 0;
 585 //    //throw SALOME_Exception(LOCALIZED("different point number-opposed edge"));
 586 //  }
 587
 588   // set normalized grid on unit square in parametric domain
 589
 590   SetNormalizedGrid(aMesh, F, quad);
 591
 592   return quad;
 593 }
 594
 595 //=============================================================================
 596 /*!
 597  *
 598  */
 599 //=============================================================================
 600
 601 void StdMeshers_Quadrangle_2D::QuadDelete (FaceQuadStruct * quad)
 602 {
 603   //MESSAGE("StdMeshers_Quadrangle_2D::QuadDelete");
 604   if (quad)
 605   {
 606     for (int i = 0; i < 4; i++)
 607     {
 608       if (quad->uv_edges[i])
 609         delete [] quad->uv_edges[i];
 610       quad->edge[i].Nullify();
 611     }
 612     if (quad->uv_grid)
 613       delete [] quad->uv_grid;
 614     delete quad;
 615   }
 616 }
 617
 618 //=============================================================================
 619 /*!
 620  *
 621  */
 622 //=============================================================================
 623
 624 void StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh,
 625                                                   const TopoDS_Shape& aShape,
 626                                                   FaceQuadStruct* quad) throw (SALOME_Exception)
 627 {
 628   Unexpect aCatch(SalomeException);
 629   // Algorithme décrit dans "Génération automatique de maillages"
 630   // P.L. GEORGE, MASSON, § 6.4.1 p. 84-85
 631   // traitement dans le domaine paramétrique 2d u,v
 632   // transport - projection sur le carré unité
 633
 634 //  MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid");
 635   const TopoDS_Face& F = TopoDS::Face(aShape);
 636
 637   // 1 --- find orientation of the 4 edges, by test on extrema
 638
 639   //      max             min                    0     x1     1
 640   //     |<----north-2-------^                a3 -------------> a2
 641   //     |                   |                   ^1          1^
 642   //    west-3            east-1 =right          |            |
 643   //     |                   |         ==>       |            |
 644   //  y0 |                   | y1                |            |
 645   //     |                   |                   |0          0|
 646   //     v----south-0-------->                a0 -------------> a1
 647   //      min             max                    0     x0     1
 648   //             =down
 649   //
 650
 651   Handle(Geom2d_Curve) c2d[4];
 652   gp_Pnt2d pf[4];
 653   gp_Pnt2d pl[4];
 654   for (int i = 0; i < 4; i++)
 655   {
 656     c2d[i] = BRep_Tool::CurveOnSurface(quad->edge[i], F,
 657                                        quad->first[i], quad->last[i]);
 658     pf[i] = c2d[i]->Value(quad->first[i]);
 659     pl[i] = c2d[i]->Value(quad->last[i]);
 660     quad->isEdgeForward[i] = false;
 661   }
 662
 663   double eps2d = 1.e-3; // *** utiliser plutot TopExp::CommonVertex, puis
 664   // distances si piece fausse
 665 //  int i = 0;
 666   if ((pf[1].Distance(pl[0]) < eps2d) || (pl[1].Distance(pl[0]) < eps2d))
 667   {
 668     quad->isEdgeForward[0] = true;
 669   } else {
 670     double tmp = quad->first[0];
 671     quad->first[0] = quad->last[0];
 672     quad->last[0] = tmp;
 673     pf[0] = c2d[0]->Value(quad->first[0]);
 674     pl[0] = c2d[0]->Value(quad->last[0]);
 675   }
 676
 677   for (int i = 1; i < 4; i++)
 678   {
 679     quad->isEdgeForward[i] = (pf[i].Distance(pl[i - 1]) < eps2d);
 680     if (!quad->isEdgeForward[i])
 681     {
 682       double tmp = quad->first[i];
 683       quad->first[i] = quad->last[i];
 684       quad->last[i] = tmp;
 685       pf[i] = c2d[i]->Value(quad->first[i]);
 686       pl[i] = c2d[i]->Value(quad->last[i]);
 687       //SCRUTE(pf[i].Distance(pl[i-1]));
 688       ASSERT(pf[i].Distance(pl[i - 1]) < eps2d);
 689     }
 690   }
 691   //SCRUTE(pf[0].Distance(pl[3]));
 692   ASSERT(pf[0].Distance(pl[3]) < eps2d);
 693
 694 //   for (int i=0; i<4; i++)
 695 //     {
 696 //       SCRUTE(quad->isEdgeForward[i]);
 697 //       MESSAGE(" -first "<<i<<" "<<pf[i].X()<<" "<<pf[i].Y());
 698 //       MESSAGE(" -last  "<<i<<" "<<pl[i].X()<<" "<<pl[i].Y());
 699 //     }
 700
 701   // 2 --- load 2d edge points (u,v) with orientation and value on unit square
 702
 703   bool loadOk = true;
 704   for (int i = 0; i < 2; i++)
 705   {
 706     quad->uv_edges[i] = LoadEdgePoints(aMesh, F, quad->edge[i],
 707                                        quad->first[i], quad->last[i]);
 708     if (!quad->uv_edges[i]) loadOk = false;
 709     //    quad->isEdgeForward[i]);
 710   }
 711
 712   for (int i = 2; i < 4; i++)
 713   {
 714     quad->uv_edges[i] = LoadEdgePoints(aMesh, F, quad->edge[i],
 715                                        quad->last[i], quad->first[i]);
 716     if (!quad->uv_edges[i]) loadOk = false;
 717     //    !quad->isEdgeForward[i]);
 718   }
 719
 720   if (!loadOk)
 721   {
 722 //    MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid - LoadEdgePoints failed");
 723     QuadDelete( quad );
 724     quad = 0;
 725     return;
 726   }
 727   // 3 --- 2D normalized values on unit square [0..1][0..1]
 728
 729 //  int nbdown = quad->nbPts[0];
 730 //  int nbright = quad->nbPts[1];
 731   int nbhoriz  = Min(quad->nbPts[0], quad->nbPts[2]);
 732   int nbvertic = Min(quad->nbPts[1], quad->nbPts[3]);
 733 //  MESSAGE("nbhoriz, nbvertic = " << nbhoriz << nbvertic);
 734
 735   quad->isEdgeOut[0] = (quad->nbPts[0] > quad->nbPts[2]);
 736   quad->isEdgeOut[1] = (quad->nbPts[1] > quad->nbPts[3]);
 737   quad->isEdgeOut[2] = (quad->nbPts[2] > quad->nbPts[0]);
 738   quad->isEdgeOut[3] = (quad->nbPts[3] > quad->nbPts[1]);
 739
 740   quad->uv_grid = new UVPtStruct[nbvertic * nbhoriz];
 741
 742   UVPtStruct *uv_grid = quad->uv_grid;
 743   UVPtStruct *uv_e0 = quad->uv_edges[0];
 744   UVPtStruct *uv_e1 = quad->uv_edges[1];
 745   UVPtStruct *uv_e2 = quad->uv_edges[2];
 746   UVPtStruct *uv_e3 = quad->uv_edges[3];
 747
 748   // nodes Id on "in" edges
 749   if (! quad->isEdgeOut[0]) {
 750     int j = 0;
 751     for (int i = 0; i < nbhoriz; i++) { // down
 752       int ij = j * nbhoriz + i;
 753       uv_grid[ij].node = uv_e0[i].node;
 754     }
 755   }
 756   if (! quad->isEdgeOut[1]) {
 757     int i = nbhoriz - 1;
 758     for (int j = 0; j < nbvertic; j++) { // right
 759       int ij = j * nbhoriz + i;
 760       uv_grid[ij].node = uv_e1[j].node;
 761     }
 762   }
 763   if (! quad->isEdgeOut[2]) {
 764     int j = nbvertic - 1;
 765     for (int i = 0; i < nbhoriz; i++) { // up
 766       int ij = j * nbhoriz + i;
 767       uv_grid[ij].node = uv_e2[i].node;
 768     }
 769   }
 770   if (! quad->isEdgeOut[3]) {
 771     int i = 0;
 772     for (int j = 0; j < nbvertic; j++) { // left
 773       int ij = j * nbhoriz + i;
 774       uv_grid[ij].node = uv_e3[j].node;
 775     }
 776   }
 777
 778   // falsificate "out" edges
 779   if (quad->isEdgeOut[0]) // down
 780     uv_e0 = MakeEdgePoints
 781       (aMesh, F, quad->edge[0], quad->first[0], quad->last[0], nbhoriz - 1);
 782   else if (quad->isEdgeOut[2]) // up
 783     uv_e2 = MakeEdgePoints
 784       (aMesh, F, quad->edge[2], quad->last[2], quad->first[2], nbhoriz - 1);
 785
 786   if (quad->isEdgeOut[1]) // right
 787     uv_e1 = MakeEdgePoints
 788       (aMesh, F, quad->edge[1], quad->first[1], quad->last[1], nbvertic - 1);
 789   else if (quad->isEdgeOut[3]) // left
 790     uv_e3 = MakeEdgePoints
 791       (aMesh, F, quad->edge[3], quad->last[3], quad->first[3], nbvertic - 1);
 792
 793   // normalized 2d values on grid
 794   for (int i = 0; i < nbhoriz; i++)
 795   {
 796     for (int j = 0; j < nbvertic; j++)
 797     {
 798       int ij = j * nbhoriz + i;
 799       // --- droite i cste : x = x0 + y(x1-x0)
 800       double x0 = uv_e0[i].normParam;   // bas - sud
 801       double x1 = uv_e2[i].normParam;   // haut - nord
 802       // --- droite j cste : y = y0 + x(y1-y0)
 803       double y0 = uv_e3[j].normParam;   // gauche-ouest
 804       double y1 = uv_e1[j].normParam;   // droite - est
 805       // --- intersection : x=x0+(y0+x(y1-y0))(x1-x0)
 806       double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
 807       double y = y0 + x * (y1 - y0);
 808       uv_grid[ij].x = x;
 809       uv_grid[ij].y = y;
 810       //MESSAGE("-xy-01 "<<x0<<" "<<x1<<" "<<y0<<" "<<y1);
 811       //MESSAGE("-xy-norm "<<i<<" "<<j<<" "<<x<<" "<<y);
 812     }
 813   }
 814
 815   // 4 --- projection on 2d domain (u,v)
 816   gp_Pnt2d a0 = pf[0];
 817   gp_Pnt2d a1 = pf[1];
 818   gp_Pnt2d a2 = pf[2];
 819   gp_Pnt2d a3 = pf[3];
 820
 821   for (int i = 0; i < nbhoriz; i++)
 822   {
 823     for (int j = 0; j < nbvertic; j++)
 824     {
 825       int ij = j * nbhoriz + i;
 826       double x = uv_grid[ij].x;
 827       double y = uv_grid[ij].y;
 828       double param_0 = uv_e0[0].param + x * (uv_e0[nbhoriz - 1].param - uv_e0[0].param); // sud
 829       double param_2 = uv_e2[0].param + x * (uv_e2[nbhoriz - 1].param - uv_e2[0].param); // nord
 830       double param_1 = uv_e1[0].param + y * (uv_e1[nbvertic - 1].param - uv_e1[0].param); // est
 831       double param_3 = uv_e3[0].param + y * (uv_e3[nbvertic - 1].param - uv_e3[0].param); // ouest
 832
 833       //MESSAGE("params "<<param_0<<" "<<param_1<<" "<<param_2<<" "<<param_3);
 834       gp_Pnt2d p0 = c2d[0]->Value(param_0);
 835       gp_Pnt2d p1 = c2d[1]->Value(param_1);
 836       gp_Pnt2d p2 = c2d[2]->Value(param_2);
 837       gp_Pnt2d p3 = c2d[3]->Value(param_3);
 838
 839       double u = (1 - y) * p0.X() + x * p1.X() + y * p2.X() + (1 - x) * p3.X();
 840       double v = (1 - y) * p0.Y() + x * p1.Y() + y * p2.Y() + (1 - x) * p3.Y();
 841
 842       u -= (1 - x) * (1 - y) * a0.X() + x * (1 - y) * a1.X() +
 843         x * y * a2.X() + (1 - x) * y * a3.X();
 844       v -= (1 - x) * (1 - y) * a0.Y() + x * (1 - y) * a1.Y() +
 845         x * y * a2.Y() + (1 - x) * y * a3.Y();
 846
 847       uv_grid[ij].u = u;
 848       uv_grid[ij].v = v;
 849
 850       //MESSAGE("-uv- "<<i<<" "<<j<<" "<<uv_grid[ij].u<<" "<<uv_grid[ij].v);
 851     }
 852   }
 853 }
 854
 855 //=============================================================================
 856 /*!
 857  *  LoadEdgePoints
 858  */
 859 //=============================================================================
 860 UVPtStruct* StdMeshers_Quadrangle_2D::LoadEdgePoints (SMESH_Mesh & aMesh,
 861                                                       const TopoDS_Face& F,
 862                                                       const TopoDS_Edge& E,
 863                                                       double first, double last)
 864 //                        bool isForward)
 865 {
 866   //MESSAGE("StdMeshers_Quadrangle_2D::LoadEdgePoints");
 867
 868   SMDS_Mesh* meshDS = aMesh.GetMeshDS();
 869
 870   // --- IDNodes of first and last Vertex
 871
 872   TopoDS_Vertex VFirst, VLast;
 873   TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
 874
 875   ASSERT(!VFirst.IsNull());
 876   SMDS_NodeIteratorPtr lid = aMesh.GetSubMesh(VFirst)->GetSubMeshDS()->GetNodes();
 877   if (!lid->more())
 878   {
 879     MESSAGE ( "NO NODE BUILT ON VERTEX" );
 880     return 0;
 881   }
 882   const SMDS_MeshNode* idFirst = lid->next();
 883
 884   ASSERT(!VLast.IsNull());
 885   lid = aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes();
 886   if (!lid->more())
 887   {
 888     MESSAGE ( "NO NODE BUILT ON VERTEX" );
 889     return 0;
 890   }
 891   const SMDS_MeshNode* idLast = lid->next();
 892
 893   // --- edge internal IDNodes (relies on good order storage, not checked)
 894
 895   map<double, const SMDS_MeshNode *> params;
 896   SMDS_NodeIteratorPtr ite = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes();
 897
 898   while(ite->more())
 899   {
 900     const SMDS_MeshNode* node = ite->next();
 901     const SMDS_EdgePosition* epos =
 902       static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
 903     double param = epos->GetUParameter();
 904     params[param] = node;
 905   }
 906
 907   int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
 908   if (nbPoints != params.size())
 909   {
 910     MESSAGE( "BAD NODE ON EDGE POSITIONS" );
 911     return 0;
 912   }
 913   UVPtStruct* uvslf = new UVPtStruct[nbPoints + 2];
 914
 915   double f, l;
 916   Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
 917
 918   bool isForward = (((l - f) * (last - first)) > 0);
 919   double paramin = 0;
 920   double paramax = 0;
 921   if (isForward)
 922   {
 923     paramin = f;
 924     paramax = l;
 925     gp_Pnt2d p = C2d->Value(f); // first point = Vertex Forward
 926     uvslf[0].x = p.X();
 927     uvslf[0].y = p.Y();
 928     uvslf[0].param = f;
 929     uvslf[0].node = idFirst;
 930     //MESSAGE("__ f "<<f<<" "<<uvslf[0].x <<" "<<uvslf[0].y);
 931     map < double, const SMDS_MeshNode* >::iterator itp = params.begin();
 932     for (int i = 1; i <= nbPoints; i++) // nbPoints internal
 933     {
 934       double param = (*itp).first;
 935       gp_Pnt2d p = C2d->Value(param);
 936       uvslf[i].x = p.X();
 937       uvslf[i].y = p.Y();
 938       uvslf[i].param = param;
 939       uvslf[i].node = (*itp).second;
 940       //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[i].x <<" "<<uvslf[i].y);
 941       itp++;
 942     }
 943     p = C2d->Value(l);          // last point = Vertex Reversed
 944     uvslf[nbPoints + 1].x = p.X();
 945     uvslf[nbPoints + 1].y = p.Y();
 946     uvslf[nbPoints + 1].param = l;
 947     uvslf[nbPoints + 1].node = idLast;
 948     //MESSAGE("__ l "<<l<<" "<<uvslf[nbPoints+1].x <<" "<<uvslf[nbPoints+1].y);
 949   } else
 950   {
 951     paramin = l;
 952     paramax = f;
 953     gp_Pnt2d p = C2d->Value(l); // first point = Vertex Reversed
 954     uvslf[0].x = p.X();
 955     uvslf[0].y = p.Y();
 956     uvslf[0].param = l;
 957     uvslf[0].node = idLast;
 958     //MESSAGE("__ l "<<l<<" "<<uvslf[0].x <<" "<<uvslf[0].y);
 959     map < double, const SMDS_MeshNode* >::reverse_iterator itp = params.rbegin();
 960
 961     for (int j = nbPoints; j >= 1; j--) // nbPoints internal
 962     {
 963       double param = (*itp).first;
 964       int i = nbPoints + 1 - j;
 965       gp_Pnt2d p = C2d->Value(param);
 966       uvslf[i].x = p.X();
 967       uvslf[i].y = p.Y();
 968       uvslf[i].param = param;
 969       uvslf[i].node = (*itp).second;
 970       //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[i].x <<" "<<uvslf[i].y);
 971       itp++;
 972     }
 973     p = C2d->Value(f);          // last point = Vertex Forward
 974     uvslf[nbPoints + 1].x = p.X();
 975     uvslf[nbPoints + 1].y = p.Y();
 976     uvslf[nbPoints + 1].param = f;
 977     uvslf[nbPoints + 1].node = idFirst;
 978     //MESSAGE("__ f "<<f<<" "<<uvslf[nbPoints+1].x <<" "<<uvslf[nbPoints+1].y);
 979   }
 980
 981   ASSERT(paramin != paramax);
 982   for (int i = 0; i < nbPoints + 2; i++)
 983   {
 984     uvslf[i].normParam = (uvslf[i].param - paramin) / (paramax - paramin);
 985     //SCRUTE(uvslf[i].normParam);
 986   }
 987
 988   return uvslf;
 989 }
 990
 991 //=============================================================================
 992 /*!
 993  *  MakeEdgePoints
 994  */
 995 //=============================================================================
 996 UVPtStruct* StdMeshers_Quadrangle_2D::MakeEdgePoints (SMESH_Mesh & aMesh,
 997                                                       const TopoDS_Face& F,
 998                                                       const TopoDS_Edge& E,
 999                                                       double first, double last,
1000                                                       int nb_segm)
1001 {
1002 //  MESSAGE("StdMeshers_Quadrangle_2D::MakeEdgePoints");
1003
1004   UVPtStruct* uvslf = new UVPtStruct[nb_segm + 1];
1005   list<double> params;
1006
1007   // --- edge internal points
1008   double fi, li;
1009   Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, fi, li);
1010   if (!Curve.IsNull()) {
1011     try {
1012       GeomAdaptor_Curve C3d (Curve);
1013       double length = EdgeLength(E);
1014       double eltSize = length / nb_segm;
1015       GCPnts_UniformAbscissa Discret (C3d, eltSize, fi, li);
1016       if (!Discret.IsDone()) return false;
1017       int NbPoints = Discret.NbPoints();
1018       for (int i = 1; i <= NbPoints; i++) {
1019         double param = Discret.Parameter(i);
1020         params.push_back(param);
1021       }
1022     }
1023     catch (Standard_Failure) {
1024       return 0;
1025     }
1026   }
1027   else
1028   {
1029     // Edge is a degenerated Edge
1030     BRep_Tool::Range(E, fi, li);
1031     double du = (li - fi) / nb_segm;
1032     for (int i = 1; i <= nb_segm + 1; i++)
1033     {
1034       double param = fi + (i - 1) * du;
1035       params.push_back(param);
1036     }
1037   }
1038
1039   double f, l;
1040   Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
1041   ASSERT(f != l);
1042
1043   bool isForward = (((l - f) * (last - first)) > 0);
1044   if (isForward) {
1045     list<double>::iterator itU = params.begin();
1046     for (int i = 0; i <= nb_segm; i++) // nbPoints internal
1047     {
1048       double param = *itU;
1049       gp_Pnt2d p = C2d->Value(param);
1050       uvslf[i].x = p.X();
1051       uvslf[i].y = p.Y();
1052       uvslf[i].param = param;
1053       uvslf[i].normParam = (param - f) / (l - f);
1054       itU++;
1055     }
1056   } else {
1057     list<double>::reverse_iterator itU = params.rbegin();
1058     for (int j = nb_segm; j >= 0; j--) // nbPoints internal
1059     {
1060       double param = *itU;
1061       int i = nb_segm - j;
1062       gp_Pnt2d p = C2d->Value(param);
1063       uvslf[i].x = p.X();
1064       uvslf[i].y = p.Y();
1065       uvslf[i].param = param;
1066       uvslf[i].normParam = (param - l) / (f - l);
1067       itU++;
1068     }
1069   }
1070
1071   return uvslf;
1072 }
1073
1074
1075 //=============================================================================
1076 /*!
1077  *
1078  */
1079 //=============================================================================
1080
1081 ostream & StdMeshers_Quadrangle_2D::SaveTo(ostream & save)
1082 {
1083   return save;
1084 }
1085
1086 //=============================================================================
1087 /*!
1088  *
1089  */
1090 //=============================================================================
1091
1092 istream & StdMeshers_Quadrangle_2D::LoadFrom(istream & load)
1093 {
1094   return load;
1095 }
1096
1097 //=============================================================================
1098 /*!
1099  *
1100  */
1101 //=============================================================================
1102
1103 ostream & operator <<(ostream & save, StdMeshers_Quadrangle_2D & hyp)
1104 {
1105   return hyp.SaveTo( save );
1106 }
1107
1108 //=============================================================================
1109 /*!
1110  *
1111  */
1112 //=============================================================================
1113
1114 istream & operator >>(istream & load, StdMeshers_Quadrangle_2D & hyp)
1115 {
1116   return hyp.LoadFrom( load );
1117 }