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