src/HEXABLOCK/HexElements_piq.cxx

   1
   2 // C++ : Table d'hexaedres (Evol Versions 3)
   3
   4 // Copyright (C) 2009-2012  CEA/DEN, EDF R&D
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include "HexElements.hxx"
  24
  25 #include "HexVector.hxx"
  26 #include "HexVertex.hxx"
  27 #include "HexEdge.hxx"
  28 #include "HexDiagnostics.hxx"
  29
  30 #include <cmath>
  31 #include <map>
  32
  33 BEGIN_NAMESPACE_HEXA
  34
  35 static bool db = false;
  36
  37 // --------------------------------------------------------
  38 struct PatQuad
  39 {
  40    int   q_edge [QUAD4];
  41    Quad* refer;
  42 };
  43 // --------------------------------------------------------
  44 struct PatVertex
  45 {
  46    double  v_x, v_y;
  47    Vertex* refer;
  48 };
  49 // --------------------------------------------------------
  50 struct PatEdge
  51 {
  52    int   v_amont, v_aval;
  53    Edge* refer;
  54    int   nbr_refer;
  55 };
  56
  57 // --------------------------------------------------------
  58 class Pattern
  59 {
  60 friend class Elements;
  61 public  :
  62     int initialize (Vertex* v1, Vertex* v2, Vertex* v3);
  63     int addQuad    (Quad* quad);
  64     int verify     (int &nbed, int &nbver);
  65
  66 private :
  67     int addEdge   (Edge*   edge);
  68     int addVertex (Vertex* vertex);
  69
  70 private :
  71    enum EnumProj { ProjXY, ProjYZ, ProjZX };
  72
  73    vector <PatVertex> pat_vertex;
  74    vector <PatEdge  > pat_edge;
  75    vector <PatQuad  > pat_quad;
  76
  77    int      nbr_vertex, nbr_edges, nbr_quads;
  78    double   determinant;
  79    EnumProj projection;
  80    Real3    base_i, base_j, origine;
  81
  82    int pos_vertex4;
  83    int pos_edge3, pos_edge4;
  84 };
  85
  86 // ====================================================== initialize
  87 int Pattern::initialize (Vertex* v1, Vertex* v2, Vertex* v3)
  88 {
  89    nbr_vertex  = nbr_edges = nbr_quads = 0;
  90    projection  = ProjXY;
  91    determinant = 1;
  92    base_i[0]   = base_i[1] =  base_i[2] = 0;
  93    base_j[0]   = base_j[1] =  base_j[2] = 0;
  94
  95    if (v1==NULL || v2==NULL || v3==NULL)
  96       return HERR;
  97
  98    addVertex (v1);
  99    addVertex (v2);
 100    addVertex (v3);
 101    if (nbr_vertex!=3)
 102       return HERR;
 103
 104    Document* doc = v1->dad ();
 105    Edge* edc1  = doc->findEdge (v1, v2);
 106    Edge* edc2  = doc->findEdge (v2, v3);
 107    if (edc1==NULL || edc2==NULL)
 108       return HERR;
 109
 110    addEdge (edc1);
 111    addEdge (edc2);
 112    if (nbr_edges!=2)
 113       return HERR;
 114
 115    pat_edge[0].nbr_refer = pat_edge[1].nbr_refer = 0;
 116
 117    pat_vertex[0].v_x = pat_vertex[1].v_x = 0;
 118    pat_vertex[2].v_x = 1;
 119    pat_vertex[0].v_y = 1;
 120    pat_vertex[1].v_y = pat_vertex[2].v_y = 0;
 121
 122    Real3 pb, pc;
 123    v1->getPoint (pc);
 124    v2->getPoint (origine);
 125    v3->getPoint (pb);
 126
 127    calc_vecteur (origine, pb, base_i);
 128    calc_vecteur (origine, pc, base_j);
 129
 130    if (db)
 131       {
 132       PutCoord (origine);
 133       PutCoord (base_i);
 134       PutCoord (base_j);
 135       PutData (determinant);
 136       }
 137
 138 /* ******************************
 139  * AB = pu.vI + pv.vJ
 140  *
 141  * vx = pu.vI[0] + pv.vJ[0]   (1) *+vJ[1] ) };
 142  *
 143  * (1 & 2 ) : pu =  (vx.vJ[1] - vy.vJ[0]) / detxy
 144  *            pv = -(vx.vI[1] - vy.vI[0]) / detxy
 145  *x
 146  * (2 & 3 ) : pu =  (vy.vJ[2] - vz.vJ[1]) / detyz
 147  *            pv = -(vy.vI[2] - vz.vI[1]) / detyz
 148  *
 149  * (3 & 1 ) : pu =  (vz.vJ[0] - vx.vJ[2]) / detzx
 150  *            pv = -(vz.vI[0] - vx.vI[2]) / detzx
 151  *
 152  *  Les 3 systemes d'equations sont valides.
 153  *  On va choisir celui dont la valeur absolue du determinant est maximale
 154    ****************************** */
 155
 156    double detxy = base_i[dir_x]*base_j[dir_y] - base_i[dir_y]*base_j[dir_x];
 157    double detyz = base_i[dir_y]*base_j[dir_z] - base_i[dir_z]*base_j[dir_y];
 158    double detzx = base_i[dir_z]*base_j[dir_x] - base_i[dir_x]*base_j[dir_z];
 159
 160    determinant = detxy;
 161    projection  = ProjXY;
 162
 163    if (fabs (detyz) > fabs (determinant))
 164       {
 165       determinant = detyz;
 166       projection  = ProjYZ;
 167       }
 168    if (fabs (detzx) > fabs (determinant))
 169       {
 170       determinant = detzx;
 171       projection  = ProjZX;
 172       }
 173
 174    return HOK;
 175 }
 176 // ====================================================== verify
 177 int Pattern::verify (int &nbed, int &nbver)
 178 {
 179    nbed  = nbr_edges;
 180    nbver = nbr_vertex;
 181    pos_edge3 = pos_edge4 = pos_vertex4 = NOTHING;
 182
 183    if (pat_edge[0].nbr_refer!=1 || pat_edge[1].nbr_refer!=1)
 184       return HERR;
 185
 186    for (int nro=2 ; nro<nbr_edges; nro++)
 187        {
 188        if (pat_edge[nro].nbr_refer==1)
 189           {
 190           int pv1 = pat_edge[nro].v_amont;
 191           int pv2 = pat_edge[nro].v_aval;
 192
 193           if (pv1==2)
 194              {
 195              pos_edge3   = nro;
 196              if (pos_vertex4 == NOTHING)
 197                  pos_vertex4 =  pv2;
 198              else if (pos_vertex4 != pv2)
 199                  {
 200                  return HERR;
 201                  }
 202              }
 203           else if (pv2==2)
 204              {
 205              pos_edge3   = nro;
 206              if (pos_vertex4 == NOTHING)
 207                  pos_vertex4 =  pv1;
 208              else if (pos_vertex4 != pv1)
 209                  {
 210                  return HERR;
 211                  }
 212              }
 213           else if (pv1==0)
 214              {
 215              pos_edge4   = nro;
 216              if (pos_vertex4 == NOTHING)
 217                  pos_vertex4 =  pv2;
 218              else if (pos_vertex4 != pv2)
 219                  {
 220                  return HERR;
 221                  }
 222              }
 223           else if (pv2==0)
 224              {
 225              pos_edge4   = nro;
 226              if (pos_vertex4 == NOTHING)
 227                  pos_vertex4 =  pv1;
 228              else if (pos_vertex4 != pv1)
 229                  {
 230                  return HERR;
 231                  }
 232              }
 233           }
 234        }
 235
 236    if (pos_edge3==NOTHING || pos_edge4==NOTHING)
 237       return HERR;
 238
 239    return HOK;
 240 }
 241 // ====================================================== addQuad
 242 int Pattern::addQuad (Quad* elt)
 243 {
 244    if (elt==NULL)
 245       return HERR;
 246
 247    PatQuad quad;
 248    quad.refer = elt;
 249    for (int nro=0; nro<QUAD4 ; nro++)
 250        {
 251        Edge* edge = elt->getEdge (nro);
 252        quad.q_edge [nro] = addEdge (edge);
 253        }
 254
 255    pat_quad.push_back (quad);
 256    nbr_quads++;
 257    return HOK;
 258 }
 259 // ====================================================== addEdge
 260 int Pattern::addEdge (Edge* elt)
 261 {
 262    for (int nro=0; nro<nbr_edges ; nro++)
 263        {
 264        if (elt==pat_edge [nro].refer)
 265           {
 266           pat_edge [nro].nbr_refer++;
 267           return nro;
 268           }
 269        }
 270
 271    PatEdge edge;
 272    edge.nbr_refer = 1;
 273    edge.refer    = elt;
 274    edge.v_amont  = addVertex (elt->getVertex (V_AMONT));
 275    edge.v_aval   = addVertex (elt->getVertex (V_AVAL));
 276
 277    pat_edge.push_back (edge);
 278    nbr_edges++;
 279    return nbr_edges-1;
 280 }
 281 // ====================================================== addVertex
 282 int Pattern::addVertex (Vertex* elt)
 283 {
 284    for (int nro=0; nro<nbr_vertex ; nro++)
 285        {
 286        if (elt==pat_vertex [nro].refer)
 287           return nro;
 288        }
 289
 290    PatVertex vertex;
 291    vertex.refer = elt;
 292    double vx = elt->getX() - origine [dir_x];
 293    double vy = elt->getY() - origine [dir_y];
 294    double vz = elt->getZ() - origine [dir_z];
 295    switch (projection)
 296       {
 297       case ProjXY :  default :
 298            vertex.v_x =  (vx*base_j[dir_y] - vy*base_j[dir_x]) / determinant;
 299            vertex.v_y = -(vx*base_i[dir_y] - vy*base_i[dir_x]) / determinant;
 300            break;
 301
 302       case ProjYZ :
 303            vertex.v_y =  (vy*base_j[dir_z] - vz*base_j[dir_y]) / determinant;
 304            vertex.v_x = -(vy*base_i[dir_z] - vz*base_i[dir_y]) / determinant;
 305            break;
 306
 307       case ProjZX :
 308            vertex.v_x =  (vz*base_j[dir_x] - vx*base_j[dir_z]) / determinant;
 309            vertex.v_y = -(vz*base_i[dir_x] - vx*base_i[dir_z]) / determinant;
 310            break;
 311       }
 312
 313    if (db)
 314       printf (" Vertex nro %d : (%g, %g, %g)\t -> (%g, %g)\n",
 315                         nbr_vertex, vx, vy, vz, vertex.v_x, vertex.v_y);
 316    pat_vertex.push_back (vertex);
 317    nbr_vertex++;
 318    return nbr_vertex-1;
 319 }
 320 // -------------------------------------------------------------------
 321 // -------------------------------------------------------------------
 322 // -------------------------------------------------------------------
 323 // -------------------------------------------------------------------
 324 // ====================================================== replaceHexas
 325 int Elements::replaceHexas (Quads& liste, Vertex* p1, Vertex* c1,
 326                             Vertex* p2, Vertex* c2, Vertex* p3, Vertex* c3)
 327 {
 328     Edge* edp1  = el_root->findEdge (p1, p2);
 329     Edge* edp2  = el_root->findEdge (p2, p3);
 330     Edge* edc1  = el_root->findEdge (c1, c2);
 331     Edge* edc2  = el_root->findEdge (c2, c3);
 332
 333     Quad* quadc = el_root->findQuad (c1, c3);
 334
 335     if (edp1==NULL || edp2==NULL || edc1==NULL || edc2==NULL || quadc==NULL)
 336        {
 337        printf ("... Error in HexaBlock function \n");
 338        printf ("... doc.replace (lquads, p1,c1, p2,c2, p3,c3)\n");
 339        if (edp1==NULL)
 340            printf ("Vertices p1 and p2 don't define an edge\n");
 341        else if (edp2==NULL)
 342            printf ("Vertices p2 and p3 don't define an edge\n");
 343        else if (edc1==NULL)
 344            printf ("Vertices c1 and c2 don't define an edge\n");
 345        else if (edc2==NULL)
 346            printf ("Vertices c2 and c3 don't define an edge\n");
 347        else if (quadc==NULL)
 348            printf ("Vertices c1 and c3 don't define a quad\n");
 349        return HERR;
 350        }
 351
 352     int np = quadc->getNbrParents ();
 353     Hexa* hexac = quadc->getParent (0);
 354
 355     if (np!=1 || hexac==NULL)
 356        {
 357        printf ("... Error in HexaBlock function \n");
 358        printf ("... doc.replace (lquads, p1,c1, p2,c2, p3,c3)\n");
 359        printf ("Quad (c1,c2,c3) is not an external quad\n");
 360        return HERR;
 361        }
 362                               // Analyse du pattern
 363     int nbquads = liste.size();
 364     Pattern   pattern;
 365     int ier = pattern.initialize (p1, p2, p3);
 366     if (ier!=HOK)
 367        {
 368        printf ("... Error in HexaBlock function \n");
 369        printf ("... doc.replace (lquads, p1,c1, p2,c2, p3,c3)\n");
 370        printf ("Vertices (p1,p2,p3) don't define a virtual quad\n");
 371        return HERR;
 372        }
 373
 374     for (int nq=0 ; nq<nbquads ; nq++)
 375         {
 376         ier = pattern.addQuad (liste[nq]);
 377         if (ier!=HOK)
 378            {
 379            printf ("... Error in HexaBlock function \n");
 380            printf ("... doc.replace (lquads, p1,c1, p2,c2, p3,c3)\n");
 381            printf ("Quad nr %d of the list is NULL\n", nq+1);
 382            return HERR;
 383            }
 384         }
 385     int nbv  = 0;
 386     int nbed = 0;
 387     int nbh  = 0;
 388     ier = pattern.verify (nbed, nbv);
 389     if (ier!=HOK)
 390        {
 391        printf ("... Error in HexaBlock function \n");
 392        printf ("... doc.replace (lquads, p1,c1, p2,c2, p3,c3)\n");
 393        printf ("Vertices (p1,p2,p3) don't define a virtual quad\n");
 394        return HERR;
 395        }
 396                               // Analyse de la cible
 397     Quads pil_quad;
 398     Hexas pil_hexa;
 399
 400     pil_hexa.push_back (hexac);
 401     pil_quad.push_back (quadc);
 402     bool more = true;
 403                               // Constitution de la pile des hexaedres
 404     while (more)
 405           {
 406           nbh ++;
 407           Quad* oppos = hexac->getOpposedQuad (quadc);
 408           pil_quad.push_back (oppos);
 409           more = oppos->getNbrParents() == 2;
 410
 411           if (more)
 412              {
 413              if (oppos->getParent (0)==hexac)
 414                  hexac = oppos->getParent(1);
 415              else
 416                  hexac = oppos->getParent(0);
 417              pil_hexa.push_back (hexac);
 418              }
 419           quadc = oppos;
 420           }
 421
 422     resize (GR_REPLACE, nbquads, nbh, nbv, nbed);
 423
 424                 // 1) Constitution des 4 coins de la cible
 425                 // 2) Quadriller la face du haut
 426                 // 3) Trouver les 4 coins de la face du bas
 427                 // 4) Quadriller la face du bas
 428                 // 6) Decouper l'hexaedre defini par ces 2 faces
 429                 // 7) S'il y a un hexaedre dessous goto 3)
 430
 431     Vertex* tvert[QUAD4] = { c1, c2, c3, pil_quad[0]->getOpposVertex(c2) };
 432     replaceQuad (1, &pattern, pil_quad[0], tvert);
 433     for (int nh=1 ; nh<gr_hauteur ; nh++)
 434         {
 435         for (int nro=0 ; nro<QUAD4 ; nro++)
 436             {
 437             Edge*  edv  = pil_hexa[nh-1]->getPerpendicularEdge (pil_quad[nh-1],
 438                                                                 tvert[nro]);
 439             tvert [nro] = edv ->opposedVertex(tvert[nro]);
 440             }
 441         replaceQuad (nh+1, &pattern, pil_quad[nh], tvert);
 442         replaceHexa (nh,   &pattern, pil_hexa[nh-1]);
 443         }
 444
 445     for (int nh=1 ; nh<=nbh ; nh++)
 446         pil_quad[nh]->remove ();
 447
 448     extrudeQuad (&pattern);
 449     replaceHexa (0, &pattern, NULL);
 450     return HOK;
 451 }
 452 // ====================================================== repVertex
 453 void Elements::repVertex (int nh, int nro, Vertex* elt)
 454 {
 455    int addr = nh * pat_nbvertex + nro;
 456    if (tab_vertex[addr] == NULL)
 457        tab_vertex[addr]  = elt;
 458 }
 459 // ====================================================== repVertex
 460 Vertex* Elements::repVertex (int nh, int nro, double px, double py, double pz)
 461 {
 462    int addr = nh * pat_nbvertex + nro;
 463
 464    if (tab_vertex[addr] == NULL)
 465        tab_vertex[addr] = el_root->addVertex (px, py, pz);
 466
 467    return tab_vertex[addr];
 468 }
 469 // ====================================================== repVertex
 470 Vertex* Elements::repVertex (int nh, int nro)
 471 {
 472    int addr = nh * pat_nbvertex + nro;
 473    return tab_vertex[addr];
 474 }
 475 // ====================================================== repEdgeH
 476 void Elements::repEdgeH (int nh, int nro, Edge* elt)
 477 {
 478    int addr = nh * (pat_nbedges + pat_nbvertex) + nro;
 479    if (tab_edge[addr] == NULL)
 480        tab_edge[addr]  = elt;
 481 }
 482 // ====================================================== repEdgeH
 483 Edge* Elements::repEdgeH (int nh, int nro, Vertex* v1, Vertex* v2)
 484 {
 485    int addr = nh * (pat_nbedges + pat_nbvertex) + nro;
 486    if (tab_edge[addr] == NULL)
 487        tab_edge[addr]  = el_root->addEdge (v1, v2);
 488    return tab_edge[addr];
 489 }
 490 // ====================================================== repEdgeH
 491 Edge* Elements::repEdgeH (int nh, int nro)
 492 {
 493    int addr = nh * (pat_nbedges + pat_nbvertex) + nro;
 494    return tab_edge[addr];
 495 }
 496 // ====================================================== repEdgeV
 497 void Elements::repEdgeV (int nh, int nro, Edge* elt)
 498 {
 499    int addr = nh * (pat_nbedges + pat_nbvertex) + pat_nbedges + nro;
 500    if (tab_edge[addr] == NULL)
 501        tab_edge[addr]  = elt;
 502 }
 503 // ====================================================== repEdgeV
 504 Edge* Elements::repEdgeV (int nh, int nro, Vertex* v1, Vertex* v2)
 505 {
 506    int addr = nh * (pat_nbedges + pat_nbvertex) + pat_nbedges + nro;
 507    if (tab_edge[addr] == NULL)
 508        tab_edge[addr]  = el_root->addEdge (v1, v2);
 509    return tab_edge[addr];
 510 }
 511 // ====================================================== repEdgeV
 512 Edge* Elements::repEdgeV (int nh, int nro)
 513 {
 514    int addr = nh * (pat_nbedges + pat_nbvertex) + pat_nbedges + nro;
 515    return tab_edge[addr];
 516 }
 517 // ====================================================== repQuadH
 518 void Elements::repQuadH (int nh, int nro, Quad* elt)
 519 {
 520    int addr = nh * (nbr_orig + pat_nbedges) + nro;
 521    if (tab_quad[addr] == NULL)
 522        tab_quad[addr]  = elt;
 523 }
 524 // ====================================================== repQuadH
 525 Quad* Elements::repQuadH (int nh, int nro)
 526 {
 527    int addr = nh * (nbr_orig + pat_nbedges) + nro;
 528    return tab_quad [addr];
 529 }
 530 // ====================================================== repQuadH
 531 Quad* Elements::repQuadH (int nh, int nro, Edge* e1, Edge* e2, Edge* e3,
 532                                            Edge* e4)
 533 {
 534    int addr = nh * (nbr_orig + pat_nbedges) + nro;
 535    if (tab_quad[addr] == NULL)
 536        tab_quad[addr]  = el_root->addQuad (e1, e2, e3, e4);
 537    return tab_quad [addr];
 538 }
 539 // ====================================================== repQuadV
 540 void Elements::repQuadV (int nh, int nro, Quad* elt)
 541 {
 542    int addr = nh * (nbr_orig + pat_nbedges) + nbr_orig + nro;
 543    if (tab_quad[addr] == NULL)
 544        tab_quad[addr]  = elt;
 545 }
 546 // ====================================================== repQuadV
 547 Quad* Elements::repQuadV (int nh, int nro)
 548 {
 549    int addr = nh * (nbr_orig + pat_nbedges) + nbr_orig + nro;
 550    return tab_quad [addr];
 551 }
 552 // ====================================================== repQuadV
 553 Quad* Elements::repQuadV (int nh, int nro, Edge* e1, Edge* e2, Edge* e3,
 554                                            Edge* e4)
 555 {
 556    int addr = nh * (nbr_orig + pat_nbedges) + nbr_orig + nro;
 557    if (tab_quad[addr] == NULL)
 558        tab_quad[addr]  = el_root->addQuad (e1, e2, e3, e4);
 559    return tab_quad [addr];
 560 }
 561 // ====================================================== repHexa
 562 Hexa* Elements::repHexa (int nh, int nro, Quad* qa, Quad* qb, Quad* qc,
 563                                           Quad* qd, Quad* qe, Quad* qf)
 564 {
 565    int addr = nh * nbr_orig  + nro;
 566    if (tab_hexa[addr] == NULL)
 567        tab_hexa[addr]  = el_root->addHexa (qa, qb, qc, qd, qe, qf);
 568    return tab_hexa [addr];
 569 }
 570 // ====================================================== replaceQuad
 571 // ==== Creation des quads horizontaux
 572 int Elements::replaceQuad (int nh, Pattern* pat, Quad* quad, Vertex* tvert[])
 573 {
 574    Vertex* pvert[QUAD4] = { tvert[0],  tvert[1],  tvert[2],  tvert[3] };
 575
 576    int vnro [4] = { 0, 1, 2, pat->pos_vertex4};
 577    int ednro[4] = { 0, 1, pat->pos_edge3, pat->pos_edge4 };
 578
 579                               // Enregistrement des vertex & edges existant
 580    for (int nro=0 ; nro<QUAD4 ; nro++)
 581        {
 582        Vertex* vh    = pvert [nro];
 583        Vertex* vh1   = pvert [(nro+1) MODULO QUAD4];
 584        Edge*   edh   = quad->findEdge (vh, vh1);
 585
 586        repVertex (nh,   vnro  [nro], vh);
 587        repEdgeH  (nh,   ednro [nro], edh);
 588        }
 589
 590    Real3 orig, ib, jb, pb, pc;
 591                               // Creation des vertex
 592    tvert[0]->getPoint (pc);
 593    tvert[1]->getPoint (orig);
 594    tvert[2]->getPoint (pb);
 595    calc_vecteur (orig, pb, ib);
 596    calc_vecteur (orig, pc, jb);
 597
 598    for (int nro=0 ; nro<pat->nbr_vertex ; nro++)
 599        {
 600        double lambda = pat->pat_vertex [nro].v_x;
 601        double mu     = pat->pat_vertex [nro].v_y;
 602        double px = orig[dir_x] + lambda*ib[dir_x] + mu*jb[dir_x];
 603        double py = orig[dir_y] + lambda*ib[dir_y] + mu*jb[dir_y];
 604        double pz = orig[dir_z] + lambda*ib[dir_z] + mu*jb[dir_z];
 605        repVertex (nh, nro, px, py, pz);
 606        }
 607                               // Creation des edges horizontaux
 608    for (int nro=0 ; nro<pat->nbr_edges ; nro++)
 609        {
 610        int nv1 = pat->pat_edge [nro].v_amont;
 611        int nv2 = pat->pat_edge [nro].v_aval;
 612        Vertex* v1 = repVertex (nh, nv1);
 613        Vertex* v2 = repVertex (nh, nv2);
 614        repEdgeH (nh, nro, v1, v2);
 615        }
 616                               // Creation des quads horizontaux
 617    for (int nro=0 ; nro<pat->nbr_quads ; nro++)
 618        {
 619        Edge* eda = repEdgeH (nh, pat->pat_quad [nro].q_edge [0]);
 620        Edge* edb = repEdgeH (nh, pat->pat_quad [nro].q_edge [1]);
 621        Edge* edc = repEdgeH (nh, pat->pat_quad [nro].q_edge [2]);
 622        Edge* edd = repEdgeH (nh, pat->pat_quad [nro].q_edge [3]);
 623
 624        repQuadH (nh,  nro, eda, edb, edc, edd);
 625        }
 626    return HOK;
 627 }
 628 // ====================================================== extrudeQuad
 629 // ==== Creation des quads horizontaux
 630 int Elements::extrudeQuad (Pattern* pat)
 631 {
 632                               // Creation des vertex de niveau 0
 633    for (int nro=0 ; nro<pat->nbr_vertex ; nro++)
 634        {
 635        Vertex* v1 = repVertex (1, nro);
 636        Vertex* v2 = repVertex (2, nro);
 637        double px = 2*v1->getX() - v2->getX();
 638        double py = 2*v1->getY() - v2->getY();
 639        double pz = 2*v1->getZ() - v2->getZ();
 640        repVertex (0, nro, px, py, pz);
 641        }
 642                               // Creation des edges horizontaux
 643    for (int nro=0 ; nro<pat->nbr_edges ; nro++)
 644        {
 645        int nv1 = pat->pat_edge [nro].v_amont;
 646        int nv2 = pat->pat_edge [nro].v_aval;
 647        Vertex* v1 = repVertex (0, nv1);
 648        Vertex* v2 = repVertex (0, nv2);
 649        repEdgeH (0, nro, v1, v2);
 650        }
 651                               // Creation des quads horizontaux
 652    for (int nro=0 ; nro<pat->nbr_quads ; nro++)
 653        {
 654        Edge* eda = repEdgeH (0, pat->pat_quad [nro].q_edge [0]);
 655        Edge* edb = repEdgeH (0, pat->pat_quad [nro].q_edge [1]);
 656        Edge* edc = repEdgeH (0, pat->pat_quad [nro].q_edge [2]);
 657        Edge* edd = repEdgeH (0, pat->pat_quad [nro].q_edge [3]);
 658
 659        repQuadH (0,  nro, eda, edb, edc, edd);
 660        }
 661    return HOK;
 662 }
 663 // ====================================================== replaceHexa
 664 int Elements::replaceHexa (int nh, Pattern* pat, Hexa* hexa)
 665 {
 666    int vnro [4] = { 0, 1, 2, pat->pos_vertex4};
 667    int ednro[4] = { 0, 1, pat->pos_edge3, pat->pos_edge4 };
 668
 669                               // Enregistrement des edges & quads existants
 670    if (hexa!=NULL)
 671        {
 672        for (int nro=0 ; nro<QUAD4 ; nro++)
 673            {
 674            Vertex* vh  = repVertex (nh,   vnro [nro]);
 675            Vertex* vh1 = repVertex (nh+1, vnro [nro]);
 676            Edge*   edv = hexa->findEdge (vh, vh1);
 677            repEdgeV (nh, vnro [nro], edv);
 678
 679            Edge*   edh   = repEdgeH (nh,   ednro [nro]);
 680            Edge*   edh1  = repEdgeH (nh+1, ednro [nro]);
 681            Quad*   quadv = hexa->findQuad (edh, edh1);
 682            repQuadV (nh, ednro [nro], quadv);
 683            }
 684        }
 685                               // Creation des edges verticaux
 686    for (int nro=0 ; nro<pat->nbr_vertex ; nro++)
 687        {
 688        Vertex* v1 = repVertex (nh,   nro);
 689        Vertex* v2 = repVertex (nh+1, nro);
 690        repEdgeV (nh, nro, v1, v2);
 691        }
 692                               // Creation des quads verticaux
 693    for (int nro=0 ; nro<pat->nbr_edges ; nro++)
 694        {
 695        int nv1 = pat->pat_edge [nro].v_amont;
 696        int nv2 = pat->pat_edge [nro].v_aval;
 697
 698        Edge* eda = repEdgeH (nh,   nro);
 699        Edge* edb = repEdgeV (nh,   nv1);
 700        Edge* edc = repEdgeH (nh+1, nro);
 701        Edge* edd = repEdgeV (nh,   nv2);
 702        repQuadV (nh,  nro, eda, edb, edc, edd);
 703        }
 704                               // Creation des hexaedres
 705    for (int nro=0 ; nro<pat->nbr_quads ; nro++)
 706        {
 707        Quad* qa = repQuadH (nh,   nro);
 708        Quad* qb = repQuadH (nh+1, nro);
 709        Quad* qc = repQuadV (nh, pat->pat_quad [nro].q_edge [0]);
 710        Quad* qd = repQuadV (nh, pat->pat_quad [nro].q_edge [2]);
 711        Quad* qe = repQuadV (nh, pat->pat_quad [nro].q_edge [1]);
 712        Quad* qf = repQuadV (nh, pat->pat_quad [nro].q_edge [3]);
 713        repHexa (nh, nro, qa, qb, qc, qd, qe, qf);
 714        }
 715
 716    return HOK;
 717 }
 718 END_NAMESPACE_HEXA