src/HEXABLOCK/HexElements.cxx

   1
   2 // C++ : Grilles
   3
   4 // Copyright (C) 2009-2024  CEA, EDF
   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, or (at your option) any later version.
  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 #include "HexDocument.hxx"
  25 #include "HexVector.hxx"
  26 #include "HexVertex.hxx"
  27 #include "HexHexa.hxx"
  28 #include "HexEdge.hxx"
  29 #include "HexGlobale.hxx"
  30
  31 #include <cmath>
  32 #include <map>
  33
  34 BEGIN_NAMESPACE_HEXA
  35
  36 static bool db=on_debug();
  37
  38 // ====================================================== Constructeur
  39 Elements::Elements (Document* doc) : EltBase (doc, EL_GRID)
  40 {
  41    glob  = Globale::getInstance ();
  42
  43    grid_type  = GR_NONE;
  44    size_qx = size_ex = size_vx = size_hx = 0;
  45    size_qy = size_ey = size_vy = size_hy = 0;
  46    size_qz = size_ez = size_vz = size_hz = 0;
  47    size_qvplus = size_qhplus = size_ehplus = size_evplus = size_hplus = 0;
  48
  49    nbr_hexas  = nbr_quads  = nbr_edges  = nbr_vertex = 0;
  50    ker_vertex = 0;
  51    cyl_closed = false;
  52    cyl_fill   = false;
  53    grid_nocart = true;
  54    cyl_dispo  = CYL_NOFILL;
  55    revo_lution = false;
  56    prism_vec   = false;
  57    val_absolues = false;       // Version Hexa6
  58    cyl_phimax   = M_PI/2;
  59    under_v6     = true;
  60 }
  61 // ====================================================== Constructeur
  62 Elements::Elements (Document* doc, int nx, int ny, int nz)
  63         : EltBase (doc, EL_GRID)
  64 {
  65    glob  = Globale::getInstance ();
  66
  67    grid_type  = GR_NONE;
  68    size_qx = size_ex = size_vx = size_hx = 0;
  69    size_qy = size_ey = size_vy = size_hy = 0;
  70    size_qz = size_ez = size_vz = size_hz = 0;
  71    size_qvplus = size_qhplus = size_ehplus = size_evplus = size_hplus = 0;
  72
  73    nbr_hexas  = nbr_quads  = nbr_edges  = nbr_vertex = 0;
  74    cyl_closed = true;
  75    cyl_fill   = false;
  76    cyl_dispo  = CYL_NOFILL;
  77
  78    resize (GR_CYLINDRIC, nx, ny, nz);
  79    cyl_closed = true;
  80    grid_geom  = NULL;
  81    cyl_phimax   = M_PI/2;
  82 }
  83 // ====================================================== Constructeur (clonage)
  84 Elements::Elements (Elements* orig) : EltBase (orig->el_root)
  85 {
  86    glob = Globale::getInstance ();
  87
  88    grid_type  = orig->grid_type;
  89    cyl_closed = orig->cyl_closed;
  90    cyl_fill   = orig->cyl_fill;
  91    cyl_dispo  = orig->cyl_dispo;
  92
  93    size_qx = size_ex = size_vx = size_hx = 0;
  94    size_qy = size_ey = size_vy = size_hy = 0;
  95    size_qz = size_ez = size_vz = size_hz = 0;
  96    size_qvplus = size_qhplus = size_ehplus = size_evplus = size_hplus = 0;
  97    nbr_hexas  = nbr_quads  = nbr_edges  = nbr_vertex = 0;
  98
  99    resize (orig->grid_type, orig->size_hx, orig->size_hy, orig->size_hz);
 100    cyl_closed = orig->cyl_closed;
 101    cyl_phimax = orig->cyl_phimax;
 102 }
 103 // ====================================================== resize
 104 void Elements::resize (EnumGrid type, int nx, int ny, int nz, int nplus)
 105 {
 106    grid_type   = type;
 107    grid_nocart = true;
 108
 109    switch (grid_type)
 110       {
 111       case GR_CARTESIAN :
 112       case GR_CYLINDRIC :
 113       default :
 114            grid_nocart = false;
 115            size_hx = std::max (nx, 1);
 116            size_hy = std::max (ny, 1);
 117            size_hz = std::max (nz, 1);
 118
 119            size_qx = size_ex = size_vx = size_hx + 1;
 120            size_qy = size_ey = size_vy = size_hy + 1;
 121            size_qz = size_ez = size_vz = size_hz + 1;
 122
 123            nbr_hexas  = size_hx * size_hy * size_hz;
 124            nbr_quads  = size_qx * size_qy * size_qz * DIM3;
 125            nbr_edges  = size_ex * size_ey * size_ez * DIM3;
 126            nbr_vertex = size_vx * size_vy * size_vz;
 127            break;
 128
 129       case GR_SPHERIC :
 130            size_qx = size_ex = size_vx = size_hx = 0;
 131            size_qy = size_ey = size_vy = size_hy = 0;
 132            size_qz = size_ez = size_vz = size_hz = 0;
 133            nbr_quads  = nbr_edges  = nbr_vertex = 0;
 134            gr_rayon = std::max (nx, 1);
 135
 136            nbr_hexas = 1 +  gr_rayon*HQ_MAXI;
 137            tab_hexa.clear ();
 138            ker_vertex = nbr_vertex;
 139            return;
 140
 141       case GR_JOINT :
 142            nbr_orig = std::max (nx, 1);
 143            gr_hauteur  = std::max (ny, 1);
 144            size_hx  = nbr_orig;
 145            size_hy  = 1;
 146            size_hz  = gr_hauteur;
 147
 148            nbr_hexas  = nbr_orig * gr_hauteur;
 149            nbr_vertex = nbr_hexas * QUAD4;
 150            nbr_vertex = nbr_orig * (gr_hauteur+1)*QUAD4;
 151            nbr_quads  = nbr_vertex;
 152            nbr_edges  = 2*nbr_vertex;
 153            break;
 154
 155       case GR_REPLACE :
 156            nbr_orig    = std::max (nx, 1);    // nb quads du pattern
 157            gr_hauteur  = ny + 1;              // Hauteur des hexas
 158            pat_nbvertex  = std::max (nz, 1);    // nb vertex du pattern
 159            pat_nbedges   = std::max (nplus, 1); // nb edges du pattern
 160            size_hx  = nbr_orig;
 161            size_hy  = 1;
 162            size_hz  = gr_hauteur;
 163
 164            nbr_hexas  = nbr_orig  * gr_hauteur;
 165            nbr_vertex = pat_nbvertex * (gr_hauteur+1);
 166            nbr_edges  = pat_nbedges * (gr_hauteur+1) + pat_nbvertex*gr_hauteur;
 167            nbr_quads  = nbr_orig    * (gr_hauteur+1) + pat_nbedges *gr_hauteur;
 168            break;
 169
 170       case GR_BICYL :
 171            cyl_closed = true;
 172            size_hx = 1;
 173            size_hy = 8;
 174            size_hz = 4;
 175
 176            size_qx = size_ex = size_vx = size_hx + 1;
 177            size_qy = size_ey = size_vy = size_hy + 1;
 178            size_qz = size_ez = size_vz = size_hz + 1;
 179
 180            nbr_hexas  = size_hx * size_hy * size_hz;
 181            nbr_quads  = size_qx * size_qy * size_qz * DIM3;
 182            nbr_edges  = size_ex * size_ey * size_ez * DIM3;
 183            nbr_vertex = size_vx * size_vy * size_vz;
 184            break;
 185       }
 186
 187    tab_hexa  .resize (nbr_hexas,  NULL);
 188    tab_quad  .resize (nbr_quads,  NULL);
 189    tab_edge  .resize (nbr_edges,  NULL);
 190    tab_vertex.resize (nbr_vertex, NULL);
 191
 192    ker_vertex = nbr_vertex;
 193 }
 194 // ====================================================== saveVtk
 195 int Elements::saveVtk (cpchar nomfic)
 196 {
 197    DumpStart ("saveVtk", nomfic);
 198
 199    pfile    vtk = fopen (nomfic, "w");
 200    if (vtk==NULL)
 201       {
 202       DumpReturn (HERR);
 203       return HERR;
 204       }
 205                                            // -- 1) Raz numerotation precedente
 206    for (int nro=0 ; nro<nbr_hexas ; nro++)
 207        {
 208        Hexa* cell = tab_hexa[nro];
 209        if (cell!=NULL && cell->isHere())
 210            cell->razNodes ();
 211        }
 212    for (int nro=0 ; nro<size_hplus ; nro++)
 213        {
 214        Hexa* cell = ker_hexa[nro];
 215        if (cell!=NULL && cell->isHere())
 216            cell->razNodes ();
 217        }
 218
 219
 220    int nbnodes = 0;
 221    int nbcells = 0;
 222                                            // -- 2) Comptage
 223    for (int nro=0 ; nro<nbr_hexas ; nro++)
 224        {
 225        Hexa* cell = tab_hexa[nro];
 226        if (cell!=NULL && cell->isHere())
 227           {
 228           nbcells ++;
 229           nbnodes += cell->countNodes ();
 230           }
 231        }
 232
 233    for (int nro=0 ; nro<size_hplus ; nro++)
 234        {
 235        Hexa* cell = ker_hexa[nro];
 236        if (cell!=NULL && cell->isHere())
 237           {
 238           nbcells ++;
 239           nbnodes += cell->countNodes ();
 240           }
 241        }
 242
 243    fprintf (vtk, "# vtk DataFile Version 3.1\n");
 244    fprintf (vtk, "%s \n", nomfic);
 245    fprintf (vtk, "ASCII\n");
 246    fprintf (vtk, "DATASET UNSTRUCTURED_GRID\n");
 247    fprintf (vtk, "POINTS %d float\n", nbnodes);
 248
 249                                            // -- 2) Les noeuds
 250    int nronode = 0;
 251    for (int nro=0 ; nro<nbr_hexas ; nro++)
 252        {
 253        Hexa* cell = tab_hexa[nro];
 254        if (cell!=NULL && cell->isHere())
 255           cell->printNodes (vtk, nronode);
 256        }
 257    for (int nro=0 ; nro<size_hplus ; nro++)
 258        {
 259        Hexa* cell = ker_hexa[nro];
 260        if (cell!=NULL && cell->isHere())
 261           cell->printNodes (vtk, nronode);
 262        }
 263                                            // -- 2) Les hexas
 264
 265    fprintf (vtk, "CELLS %d %d\n", nbcells, nbcells*(HV_MAXI+1));
 266
 267    for (int nro=0 ; nro<nbr_hexas ; nro++)
 268        {
 269        Hexa* cell = tab_hexa[nro];
 270        if (cell!=NULL && cell->isHere())
 271           cell->printHexa (vtk);
 272        }
 273    for (int nro=0 ; nro<size_hplus ; nro++)
 274        {
 275        Hexa* cell = ker_hexa[nro];
 276        if (cell!=NULL && cell->isHere())
 277           cell->printHexa (vtk);
 278        }
 279
 280    fprintf (vtk, "CELL_TYPES %d\n", nbcells);
 281    for (int nro=0 ; nro<nbcells ; nro++)
 282        fprintf (vtk, "%d\n", HE_MAXI);
 283
 284    fclose (vtk);
 285    DumpReturn (HOK);
 286    return HOK;
 287 }
 288 // ============ = = = = = = = = = Vertices ..........
 289 // ====================================================== newEdge
 290 Edge* Elements::newEdge (Vertex* v1, Vertex* v2)
 291 {
 292    if (v1==NULL || v2==NULL)
 293       return NULL;
 294
 295    Edge* elt = new Edge (v1, v2);
 296    return elt;
 297 }
 298 // ====================================================== newQuad
 299 Quad* Elements::newQuad (Edge* e1, Edge* e2, Edge* e3, Edge* e4)
 300 {
 301    if (e1==NULL || e2==NULL || e3==NULL|| e4==NULL)
 302       return NULL;
 303
 304    Quad* elt = new Quad (e1, e2, e3, e4);
 305    return elt;
 306 }
 307 // ====================================================== newHexa
 308 Hexa* Elements::newHexa (Quad* qa, Quad* qb, Quad* qc, Quad* qd,
 309                                              Quad* qe, Quad* qf)
 310 {
 311    if (qa==NULL || qb==NULL || qc==NULL|| qd==NULL || qe==NULL|| qf==NULL)
 312       return NULL;
 313
 314    Hexa*  elt = new Hexa (qa, qb, qc, qd, qe, qf);
 315    return elt;
 316 }
 317 // ======================================================== coupler
 318 int Elements::coupler (int nquad, Quad* dest, StrOrient* orient)
 319 {
 320    Quad* orig = tab_orig [nquad];
 321
 322    setQuad (orig, dir_z, 0, 0, 0);
 323    setQuad (dest, dir_z, nquad, 0, 0);
 324
 325    int n11 = orig->indexVertex (orient->v11);
 326    int n12 = orig->indexVertex (orient->v12);
 327    int n21 = dest->indexVertex (orient->v21);
 328    int n22 = dest->indexVertex (orient->v22);
 329
 330                // ---------------- Les 4 sommets initiaux
 331    Vertex* vorig[QUAD4] = { orient->v11, orient->v12,
 332                             orig->getVertex((n11+2) MODULO QUAD4),
 333                             orig->getVertex((n12+2) MODULO QUAD4) };
 334
 335    Vertex* vdest[QUAD4] = { orient->v21, orient->v22,
 336                             dest->getVertex((n21+2) MODULO QUAD4),
 337                             dest->getVertex((n22+2) MODULO QUAD4) };
 338    if (db)
 339       {
 340       printf ("Quad nro %d : ", nquad);
 341       orig->printName (" est couple avec ");
 342       dest->printName ("\n");
 343       printf ("Orientation : (");
 344       for (int ii=0 ; ii<QUAD4 ; ii++) printf("%s ", vorig[ii]->getName());
 345       printf (")\n");
 346       printf ("           -> (");
 347       for (int ii=0 ; ii<QUAD4 ; ii++) printf("%s ", vdest[ii]->getName());
 348       printf (")\n");
 349       }
 350
 351                // ---------------- Les sommets + les aretes verticales
 352    for (int ns=0 ; ns< QUAD4 ; ns++)
 353        {
 354        Vertex* nd1 = vorig[ns];
 355        Vertex* nd2 = vdest[ns];
 356        int nref0   = nd1->getMark();
 357        tab_vertex [nroVertex (nquad,ns,0)] = nd1;
 358
 359        if (nref0==IS_NONE)
 360           {
 361           double px0 = nd1->getX();
 362           double py0 = nd1->getY();
 363           double pz0 = nd1->getZ();
 364
 365           double dxt = nd2->getX() - px0;
 366           double dyt = nd2->getY() - py0;
 367           double dzt = nd2->getZ() - pz0;
 368
 369           nd1->setMark (indVertex (nquad, ns, 0));
 370
 371           Vertex* nd = nd1;
 372           for (int nh=0 ; nh<gr_hauteur ; nh++)
 373               {
 374               double coeff = under_v6 ? ((double)(nh+1))/gr_hauteur : cum_values [nh];
 375               Vertex* ndp = nd;
 376               if (nh == gr_hauteur-1)
 377                  nd = nd2 ;
 378               else
 379                  nd = el_root->addVertex (px0 + coeff*dxt, py0 + coeff*dyt,
 380                                           pz0 + coeff*dzt);
 381               int nv  = indVertex (nquad, ns, nh);
 382               tab_vertex [nv] = nd;
 383               tab_edge   [nv] = el_root->addEdge (ndp, nd);
 384               if (db)
 385                  printf (" Edge vertical nro %d = %s = (%s, %s)\n", nv,
 386                          tab_edge[nv]->getName(),
 387                          ndp->getName(), nd->getName());
 388               }
 389           }
 390        else
 391           {
 392           for (int nh=0 ; nh<gr_hauteur ; nh++)
 393               {
 394               int nv  = indVertex (nquad, ns, nh);
 395               int nv0 = indVertex (nref0, nh);
 396               tab_vertex [nv] = tab_vertex [nv0];
 397               tab_edge   [nv] = tab_edge   [nv0];
 398               }
 399           }
 400        }
 401                // ---------------- Les quads verticaux + aretes horiz
 402    for (int ns=0 ; ns< QUAD4 ; ns++)
 403        {
 404        int next = (ns+1) MODULO QUAD4;
 405        Edge*   arete = orig->findEdge (vorig[ns], vorig[next]);
 406        int     nref0 = arete->getMark();
 407        int     nref  = indVertex (nquad, ns, 0);
 408                                   // Construction des faces & aretes H
 409        if (nref0==IS_NONE)
 410           {
 411           arete->setMark (nref);
 412           Edge* ea = arete;
 413           Edge *eb, *ec, *ed;
 414           int  nva, nvb, nha;
 415
 416           for (int nh=0 ; nh< gr_hauteur ; nh++)
 417               {
 418               nva = indVertex (nquad, ns,   nh);
 419               nvb = indVertex (nquad, next, nh);
 420               nha = nroEdgeH  (nva);
 421
 422               ed  = tab_edge [nva];
 423               ec  = ea;
 424               eb  = tab_edge [nvb];
 425               if (nh==gr_hauteur-1)
 426                  ea = dest->findEdge (vdest[ns], vdest[next]);
 427               else
 428                  ea = el_root->addEdge (tab_vertex [nva], tab_vertex [nvb]);
 429
 430               propagateAssociation (ec, ea, eb);
 431               tab_quad [nva] = newQuad (ea, eb, ec, ed);
 432               if (BadElement (tab_quad [nva]))
 433                   return HERR;
 434               tab_edge [nha] = ea;
 435               }
 436           }
 437                                   // L'arete a deja ete traitee
 438        else
 439           {
 440           for (int nh=0 ; nh<gr_hauteur ; nh++)
 441               {
 442               int nva = indVertex (nquad, ns, nh);
 443               int nha = nroEdgeH  (nva);
 444
 445               int nvb = indVertex (nref0, nh);
 446               int nhb = nroEdgeH  (nvb);
 447
 448               tab_quad [nva] = tab_quad [nvb];
 449               tab_edge [nha] = tab_edge [nhb];
 450               }
 451           }
 452        }
 453                // -------------------------------- Les Hexas
 454    Quad *fb = orig;
 455    Quad *fa, *fc, *fd, *fe, *ff;
 456    for (int nh=0 ; nh< gr_hauteur ; nh++)
 457        {
 458        fa = fb;
 459        if (nh == gr_hauteur-1)
 460           fb = dest;
 461        else
 462           fb = newQuad (tab_edge [nroEdgeH (nquad, E_A, nh)],
 463                         tab_edge [nroEdgeH (nquad, E_B, nh)],
 464                         tab_edge [nroEdgeH (nquad, E_C, nh)],
 465                         tab_edge [nroEdgeH (nquad, E_D, nh)]);
 466
 467        fc = tab_quad [indVertex (nquad, E_A, nh)];
 468        fd = tab_quad [indVertex (nquad, E_C, nh)];
 469        fe = tab_quad [indVertex (nquad, E_B, nh)];
 470        ff = tab_quad [indVertex (nquad, E_D, nh)];
 471
 472        tab_hexa [nroHexa(nquad,nh)] = newHexa (fa,fb,fc,fd,fe,ff);
 473        if (BadElement (tab_hexa [nroHexa(nquad,nh)]))
 474            return HERR;
 475        }
 476    return HOK;
 477 }
 478 // ====================================================== transform
 479 int Elements::transform (Matrix* matrice)
 480 {
 481                                            // -- 1) Raz Marques
 482    for (int nro=0 ; nro<nbr_hexas ; nro++)
 483        {
 484        Hexa* cell = tab_hexa[nro];
 485        if (cell!=NULL && cell->isHere())
 486            cell->razNodes ();
 487        }
 488                                            // -- 2) Move Nodes
 489    for (int nro=0 ; nro<nbr_hexas ; nro++)
 490        {
 491        Hexa* cell = tab_hexa[nro];
 492        if (cell!=NULL && cell->isHere())
 493            cell->moveNodes (matrice);
 494        }
 495
 496    if (nbr_hexas!=0 )
 497       return HOK;
 498                     // -- Cas pathologique : il n'y avait pas d'hexas
 499                     // -- On se rabat sur les sommets
 500
 501    for (int nro=0 ; nro<nbr_vertex ; nro++)
 502        {
 503        Vertex* node = tab_vertex[nro];
 504        if (node!=NULL && node->isHere())
 505            matrice->perform (node);
 506        }
 507
 508    return HOK;
 509 }
 510 // ====================================================== cutHexas
 511 int Elements::cutHexas  (const Edges& t_edges, int nbcuts)
 512 {
 513                                        // 1) marquage des hexas
 514    el_root->markAll (NO_USED);
 515                                        // 2) Memo noeuds
 516    std::vector <Quad*> q_amont;
 517    std::vector <Quad*> q_aval;
 518    std::map    <Vertex*, Vertex*> vis_a_vis;
 519
 520    int nbnodes = t_edges.size();
 521    std::vector <Vertex*> v_amont (nbnodes);
 522    std::vector <Vertex*> v_aval  (nbnodes);
 523
 524    int nbfaces = 0;
 525    for (int nro=0; nro<nbnodes ; nro++)
 526        {
 527        Edge* arete   = t_edges [nro];
 528        v_amont [nro] = arete->getAmont ();
 529        v_aval  [nro] = arete->getAval  ();
 530        if (db)
 531           {
 532           printf (" %3d : Edge = (", nro);
 533           v_amont[nro]->printName (", ");
 534           v_aval [nro]->printName (")\n");
 535           }
 536
 537        vis_a_vis [v_amont[nro]] = v_aval[nro];
 538        vis_a_vis [v_aval[nro]] = v_amont[nro];
 539        int nbcells   = arete->getNbrParents ();
 540
 541        for (int nq=0 ; nq<nbcells ; nq++)
 542            {
 543            Quad* quad = arete->getParent (nq);
 544            if (quad->getMark () != IS_USED)
 545               {
 546               quad->setMark (IS_USED);
 547               int nbcubes = quad->getNbrParents ();
 548               for (int nh=0 ; nh<nbcubes ; nh++)
 549                   {
 550                   Hexa* hexa = quad->getParent (nh);
 551                   if (hexa->getMark () != IS_USED)
 552                      {
 553                      hexa->setMark (IS_USED);
 554                      int namont = hexa->getBase (v_amont[nro], arete);
 555                      int naval  = glob->getOpposedQuad (namont);
 556                      q_amont.push_back (hexa->getQuad  (namont));
 557                      q_aval .push_back (hexa->getQuad  (naval ));
 558
 559                      if (db)
 560                         {
 561                         printf (" %3d : Quad = ", nbfaces);
 562                         hexa->printName (", ");
 563                         printf (" Faces = (");
 564                         hexa->getQuad (namont)->printName (", ");
 565                         hexa->getQuad (naval )->printName (")\n");
 566                         nbfaces ++;
 567                         }
 568                      }
 569                   }
 570               }
 571            }
 572        }
 573                    // ------------------- Dimensionnement
 574    int nbcells   = q_amont.size ();
 575    nbr_vertex    = nbnodes*(nbcuts+2);
 576    int nbpiliers = nbnodes*(nbcuts+1);            // aretes verticales
 577    int nbpoutres = nbcells*(nbcuts+2)*QUAD4;      // aretes horizontales
 578    nbr_edges     = nbpoutres;
 579    nbr_quads     = nbcells*(nbcuts+1)*QUAD4;      // faces Verticales
 580    nbr_hexas     = nbcells*(nbcuts+1);
 581
 582                    // ------------------- Les noeuds et les aretes verticales
 583    tab_quad.resize   (nbr_quads,  NULL);
 584    tab_edge.resize   (nbr_edges,  NULL);
 585    tab_hexa.resize   (nbr_hexas,  NULL);
 586    tab_vertex.resize (nbr_vertex, NULL);
 587    std::vector <Edge*>    tab_pilier (nbpiliers);
 588
 589    int nbinter = nbcuts + 1;
 590    for (int ned=0; ned<nbnodes ; ned++)
 591        {
 592        Edges  ass_edges;
 593        t_edges [ned]->remove ();
 594        Vertex* ndamont = v_amont [ned];
 595        Vertex* ndaval  = v_aval  [ned];
 596
 597        double lgx = ndaval->getX() - ndamont->getX();
 598        double lgy = ndaval->getY() - ndamont->getY();
 599        double lgz = ndaval->getZ() - ndamont->getZ();
 600
 601        Vertex* nd0 = tab_vertex [ned] = ndamont;
 602        for (int nc=0; nc<nbcuts ; nc++)
 603            {
 604            int nc1 = nc+1;
 605            double coeff = under_v6 ? ((double)(nc+1))/nbinter : cum_values[nc];
 606            Vertex* nd1 = el_root->addVertex (ndamont->getX() + coeff*lgx,
 607                                              ndamont->getY() + coeff*lgy,
 608                                              ndamont->getZ() + coeff*lgz);
 609            tab_vertex [nc1*nbnodes + ned] = nd1;
 610            tab_pilier [nc *nbnodes + ned] = newEdge (nd0, nd1);
 611            ass_edges.push_back (tab_pilier [nc *nbnodes + ned]);
 612            nd0 =  nd1;
 613            }
 614        tab_vertex [nbinter*nbnodes + ned] = ndaval;
 615        tab_pilier [nbcuts *nbnodes + ned] = newEdge (nd0, ndaval);
 616        ass_edges.push_back (tab_pilier[nbcuts *nbnodes + ned]);
 617        ndamont->setMark (ned);
 618        if (t_edges[ned]->isAssociated())
 619            cutAssociation (t_edges[ned], ass_edges);
 620        }
 621                    // ------------------- Les aretes horizontales
 622                    // ------------------- Les faces verticales
 623    HexDisplay (nbcells);
 624    int sizelig = nbcells*QUAD4;
 625    for (int nro=0; nro<nbcells ; nro++)
 626        {
 627        Quad* sol  = q_amont [nro];
 628        Quad* toit = q_aval  [nro];
 629        for (int ns=0; ns<QUAD4 ; ns++)
 630            {
 631            Edge* plinthe = sol->getEdge (ns);
 632            int   nmur  = nro*QUAD4 + ns;
 633            int   nmur0 = plinthe->getMark();
 634            if (nmur0 >= 0)
 635               {
 636               for (int nc=0 ; nc<nbinter ; nc++)
 637                   {
 638                   tab_edge [nc*sizelig + nmur] = tab_edge [nc*sizelig + nmur0];
 639                   tab_quad [nc*sizelig + nmur] = tab_quad [nc *sizelig + nmur0];
 640                   if (db)
 641                      {
 642                      printf (" %2d : %d quad_vertical [%02d] =", nro, ns,
 643                                                         nc*sizelig + nmur);
 644                      printf (" quad_vertical [%02d]\n",  nc*sizelig + nmur0);
 645                      }
 646                   }
 647               tab_edge [nbinter*sizelig+nmur] = tab_edge[nbinter*sizelig+nmur0];
 648               }
 649            else
 650               {
 651               plinthe->setMark (nmur);
 652               Vertex* vs1 = sol->getVertex (ns);
 653               Vertex* vs2 = sol->getVertex ((ns+1) MODULO QUAD4);
 654               int     nd1 = vs1->getMark ();
 655               int     nd2 = vs2->getMark ();
 656               Edge*   ed0 = tab_edge [nmur] = plinthe;
 657               Edge*   ed2 = NULL;
 658               Vertex* v1  = NULL;
 659               Vertex* v2  = NULL;
 660               for (int nc=0 ; nc<nbinter ; nc++)
 661                   {
 662                   int   nc1 = nc + 1;
 663                   if (nc<nbcuts)
 664                      {
 665                      v1 = tab_vertex [nc1*nbnodes + nd1];
 666                      v2 = tab_vertex [nc1*nbnodes + nd2];
 667                      ed2 = newEdge (v1, v2);
 668                      }
 669                   else
 670                      {
 671                      v1 = vis_a_vis [vs1];
 672                      v2 = vis_a_vis [vs2];
 673                      ed2 = toit->findEdge (v1, v2);
 674                      }
 675
 676                   tab_edge [nc1*sizelig + nmur] = ed2;
 677                   tab_quad [nc *sizelig + nmur] = newQuad (ed0,
 678                             tab_pilier [nc*nbnodes + nd1], ed2,
 679                             tab_pilier [nc*nbnodes + nd2]);
 680                   ed0 = ed2;
 681                   if (db)
 682                      {
 683                      printf (" %2d : %d quad_vertical [%02d] = ", nro, ns,
 684                                                          nc*sizelig + nmur);
 685                      PrintName (tab_quad [nc *sizelig + nmur]);
 686                      printf ("\n");
 687                      }
 688                   }
 689               }
 690            }
 691        }
 692                    // ------------------- Les faces horizontales
 693                    // ------------------- Les hexas
 694    // Rappel : sizelig = nbcells*QUAD4
 695    for (int nro=0; nro<nbcells ; nro++)
 696        {
 697        Quad* qa = q_amont [nro];
 698        for (int nc=0 ; nc<nbinter ; nc++)
 699            {
 700            int  quadv =  nc*nbcells*QUAD4 + nro*QUAD4;
 701            Quad* qb = q_aval[nro];
 702            if (nc<nbcuts)
 703               {
 704               int  edh   = (nc+1)*nbcells*QUAD4 + nro*QUAD4;
 705               qb = newQuad (tab_edge[edh],   tab_edge[edh+1],
 706                             tab_edge[edh+2], tab_edge[edh+3]);
 707               if (BadElement(qb))
 708                   return HERR;
 709               }
 710            Quad* qc = tab_quad [quadv];
 711            Quad* qd = tab_quad [quadv + 2];
 712            Quad* qe = tab_quad [quadv + 1];
 713            Quad* qf = tab_quad [quadv + 3];
 714            Hexa* hexa = newHexa (qa, qb, qc, qd, qe, qf);
 715            if (BadElement(hexa))
 716                return HERR;
 717            tab_hexa [nc*nbcells + nro] = hexa;
 718            qa = qb;
 719            }
 720        }
 721    return HOK;
 722 }
 723 // ====================================================== clear_associations
 724 void clear_associations (std::vector<Quad*>& table)
 725 {
 726    int nbelts = table.size();
 727    for (int nro=0 ; nro<nbelts ; nro++)
 728        {
 729        Quad* elt = table [nro];
 730        if (elt != NULL && elt->isValid())
 731            elt->clearAssociation ();
 732        }
 733 }
 734 // ====================================================== clear_associations
 735 void clear_associations (std::vector<Edge*>& table)
 736 {
 737    int nbelts = table.size();
 738    for (int nro=0 ; nro<nbelts ; nro++)
 739        {
 740        Edge* elt = table [nro];
 741        if (elt != NULL && elt->isValid())
 742            elt->clearAssociation ();
 743        }
 744 }
 745 // ====================================================== clear_associations
 746 void clear_associations (std::vector<Vertex*>& table)
 747 {
 748    int nbelts = table.size();
 749    for (int nro=0 ; nro<nbelts ; nro++)
 750        {
 751        Vertex* elt = table [nro];
 752        if (elt != NULL && elt->isValid())
 753            elt->clearAssociation ();
 754        }
 755 }
 756 // ====================================================== clearAssociation
 757 void Elements::clearAssociation  ()
 758 {
 759    clear_associations (tab_quad);
 760    clear_associations (tab_edge);
 761    clear_associations (tab_vertex);
 762
 763
 764    clear_associations (ker_hquad);
 765    clear_associations (ker_vquad);
 766    clear_associations (ker_hedge);
 767    clear_associations (ker_vedge);
 768 }
 769 // ============================================================ findVertex
 770 int Elements::findVertex (double vx, double vy, double vz)
 771 {
 772    double tol = el_root->getTolerance ();
 773    double xmin = vx - tol;
 774    double xmax = vx + tol;
 775    double ymin = vy - tol;
 776    double ymax = vy + tol;
 777    double zmin = vz - tol;
 778    double zmax = vz + tol;
 779
 780    int nbre = tab_vertex.size();
 781    for (int nro=0 ; nro<nbre ; nro++)
 782        {
 783        Vertex* node = tab_vertex [nro];
 784        if (node != NULL && node->isHere ()
 785                         && node->isin   (xmin, xmax, ymin, ymax, zmin, zmax))
 786               return nro;
 787        }
 788    return NOTHING;
 789 }
 790 // ============================================================ findQuad
 791 Quad* Elements::findQuad (Edge* e1, Edge* e2)
 792 {
 793    int nbre = tab_quad.size();
 794    for (int nro=0 ; nro<nbre ; nro++)
 795        {
 796        Quad* quad = tab_quad [nro];
 797        if (quad != NULL && quad->isHere ()
 798                         && quad->definedBy (e1, e2))
 799               return quad;
 800        }
 801    return NULL;
 802 }
 803
 804 END_NAMESPACE_HEXA