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