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