src/HEXABLOCK/HexElements.cxx

   1
   2 // C++ : Grilles
   3
   4 // Copyright (C) 2009-2013  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 #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 static bool db=false;
  35
  36 BEGIN_NAMESPACE_HEXA
  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 }
  58 // ====================================================== Constructeur
  59 Elements::Elements (Document* doc, int nx, int ny, int nz)
  60         : EltBase (doc, EL_GRID)
  61 {
  62    glob  = Globale::getInstance ();
  63
  64    grid_type  = GR_NONE;
  65    size_qx = size_ex = size_vx = size_hx = 0;
  66    size_qy = size_ey = size_vy = size_hy = 0;
  67    size_qz = size_ez = size_vz = size_hz = 0;
  68    size_qvplus = size_qhplus = size_ehplus = size_evplus = size_hplus = 0;
  69
  70    nbr_hexas  = nbr_quads  = nbr_edges  = nbr_vertex = 0;
  71    cyl_closed = true;
  72    cyl_fill   = false;
  73    cyl_dispo  = CYL_NOFILL;
  74
  75    resize (GR_CYLINDRIC, nx, ny, nz);
  76    cyl_closed = true;
  77 }
  78 // ====================================================== Constructeur (clonage)
  79 Elements::Elements (Elements* orig) : EltBase (orig->el_root)
  80 {
  81    glob = Globale::getInstance ();
  82
  83    grid_type  = orig->grid_type;
  84    cyl_closed = orig->cyl_closed;
  85    cyl_fill   = orig->cyl_fill;
  86    cyl_dispo  = orig->cyl_dispo;
  87
  88    size_qx = size_ex = size_vx = size_hx = 0;
  89    size_qy = size_ey = size_vy = size_hy = 0;
  90    size_qz = size_ez = size_vz = size_hz = 0;
  91    size_qvplus = size_qhplus = size_ehplus = size_evplus = size_hplus = 0;
  92    nbr_hexas  = nbr_quads  = nbr_edges  = nbr_vertex = 0;
  93
  94    resize (orig->grid_type, orig->size_hx, orig->size_hy, orig->size_hz);
  95    cyl_closed = orig->cyl_closed;
  96 }
  97 // ====================================================== resize
  98 void Elements::resize (EnumGrid type, int nx, int ny, int nz, int nplus)
  99 {
 100    grid_type   = type;
 101    grid_nocart = true;
 102
 103    switch (grid_type)
 104       {
 105       case GR_CARTESIAN :
 106       case GR_CYLINDRIC :
 107       default :
 108            grid_nocart = false;
 109            size_hx = std::max (nx, 1);
 110            size_hy = std::max (ny, 1);
 111            size_hz = std::max (nz, 1);
 112
 113            size_qx = size_ex = size_vx = size_hx + 1;
 114            size_qy = size_ey = size_vy = size_hy + 1;
 115            size_qz = size_ez = size_vz = size_hz + 1;
 116
 117            nbr_hexas  = size_hx * size_hy * size_hz;
 118            nbr_quads  = size_qx * size_qy * size_qz * DIM3;
 119            nbr_edges  = size_ex * size_ey * size_ez * DIM3;
 120            nbr_vertex = size_vx * size_vy * size_vz;
 121            break;
 122
 123       case GR_SPHERIC :
 124            size_qx = size_ex = size_vx = size_hx = 0;
 125            size_qy = size_ey = size_vy = size_hy = 0;
 126            size_qz = size_ez = size_vz = size_hz = 0;
 127            nbr_quads  = nbr_edges  = nbr_vertex = 0;
 128            gr_rayon = std::max (nx, 1);
 129
 130            nbr_hexas = 1 +  gr_rayon*HQ_MAXI;
 131            tab_hexa.clear ();
 132            ker_vertex = nbr_vertex;
 133            return;
 134
 135       case GR_JOINT :
 136            nbr_orig = std::max (nx, 1);
 137            gr_hauteur  = std::max (ny, 1);
 138            size_hx  = nbr_orig;
 139            size_hy  = 1;
 140            size_hz  = gr_hauteur;
 141
 142            nbr_hexas  = nbr_orig * gr_hauteur;
 143            nbr_vertex = nbr_hexas * QUAD4;
 144            nbr_vertex = nbr_orig * (gr_hauteur+1)*QUAD4;
 145            nbr_quads  = nbr_vertex;
 146            nbr_edges  = 2*nbr_vertex;
 147            break;
 148
 149       case GR_REPLACE :
 150            nbr_orig    = std::max (nx, 1);    // nb quads du pattern
 151            gr_hauteur  = ny + 1;              // Hauteur des hexas
 152            pat_nbvertex  = std::max (nz, 1);    // nb vertex du pattern
 153            pat_nbedges   = std::max (nplus, 1); // nb edges du pattern
 154            size_hx  = nbr_orig;
 155            size_hy  = 1;
 156            size_hz  = gr_hauteur;
 157
 158            nbr_hexas  = nbr_orig  * gr_hauteur;
 159            nbr_vertex = pat_nbvertex * (gr_hauteur+1);
 160            nbr_edges  = pat_nbedges * (gr_hauteur+1) + pat_nbvertex*gr_hauteur;
 161            nbr_quads  = nbr_orig    * (gr_hauteur+1) + pat_nbedges *gr_hauteur;
 162            break;
 163
 164       case GR_BICYL :
 165            cyl_closed = true;
 166            size_hx = 1;
 167            size_hy = 8;
 168            size_hz = 4;
 169
 170            size_qx = size_ex = size_vx = size_hx + 1;
 171            size_qy = size_ey = size_vy = size_hy + 1;
 172            size_qz = size_ez = size_vz = size_hz + 1;
 173
 174            nbr_hexas  = size_hx * size_hy * size_hz;
 175            nbr_quads  = size_qx * size_qy * size_qz * DIM3;
 176            nbr_edges  = size_ex * size_ey * size_ez * DIM3;
 177            nbr_vertex = size_vx * size_vy * size_vz;
 178            break;
 179       }
 180
 181    tab_hexa  .resize (nbr_hexas);
 182    tab_quad  .resize (nbr_quads);
 183    tab_edge  .resize (nbr_edges);
 184    tab_vertex.resize (nbr_vertex);
 185
 186    ker_vertex = nbr_vertex;
 187
 188    for (int nc=0 ; nc< nbr_hexas  ; nc++) tab_hexa   [nc] = NULL;
 189    for (int nc=0 ; nc< nbr_quads  ; nc++) tab_quad   [nc] = NULL;
 190    for (int nc=0 ; nc< nbr_edges  ; nc++) tab_edge   [nc] = NULL;
 191    for (int nc=0 ; nc< nbr_vertex ; nc++) tab_vertex [nc] = NULL;
 192 }
 193
 194 // ====================================================== makeCartesianGrid
 195 int Elements::makeCartesianGrid (Vertex* orig, Vector* v1, Vector* v2,
 196                    Vector* v3, int px, int py, int pz, int mx, int my, int mz)
 197 {
 198    if (BadElement (orig) || BadElement(v1) || BadElement(v2) || BadElement(v3)
 199                          || px<=0 || py<=0 || pz <= 0
 200        || v1->getNorm () <= Epsil || v2->getNorm () <= Epsil
 201        || v3->getNorm () <= Epsil)
 202       {
 203       setError ();
 204       return HERR;
 205       }
 206
 207    resize (GR_CARTESIAN, px+mx, py+my, pz+mz);
 208
 209    makeCartesianNodes (orig, v1, v2, v3, px, py, pz, mx, my, mz);
 210
 211    fillGrid ();
 212    return HOK;
 213 }
 214 // ====================================================== makeCylindricalGrid
 215 int Elements::makeCylindricalGrid (Vertex* c, Vector* b, Vector* h,
 216          double dr, double da, double dl, int nr, int na, int nl, bool fill)
 217 {
 218    if (BadElement (c) || BadElement(b) || BadElement(h)
 219        || nr<=0 || na<=0 || nl <= 0 || dr < Epsil || da < Epsil || dl < Epsil
 220        || b->getNorm () <= Epsil || h->getNorm () <= Epsil)
 221       {
 222       setError ();
 223       return HERR;
 224       }
 225    resize (GR_CYLINDRIC, nr, na, nl);
 226    cyl_closed = da >= 360.0;
 227    makeCylindricalNodes (c, b, h, dr, da, dl, nr, na, nl, fill);
 228    fillGrid ();
 229    assoCylinder (c, h, da);
 230    return HOK;
 231 }
 232 // ====================================================== makeSphericalGrid
 233 int Elements::makeSphericalGrid (Vertex* c, Vector* dv, int nb, double  k)
 234 {
 235    if (BadElement (c) || BadElement(dv)
 236        || nb<=0 || k < Epsil
 237        || dv->getDx()<=Epsil || dv->getDy()<=Epsil || dv->getDz()<=Epsil)
 238       {
 239       setError ();
 240       return HERR;
 241       }
 242
 243    resize (GR_SPHERIC, nb);
 244
 245    Vertex* i_node [HV_MAXI];    // Les noeuds de l'hexa englobant
 246    Edge*   i_edge [HE_MAXI];    // Les noeuds de l'hexa englobant
 247    Quad*   i_quad [HQ_MAXI];    // Les noeuds de l'hexa englobant
 248
 249    for (int nro=0 ; nro<HV_MAXI; nro++)
 250        {
 251        double dx = glob->CoordVertex (nro, dir_x) * dv->getDx();
 252        double dy = glob->CoordVertex (nro, dir_y) * dv->getDy();
 253        double dz = glob->CoordVertex (nro, dir_z) * dv->getDz();
 254
 255        i_node [nro] = el_root->addVertex (c->getX ()+dx, c->getY ()+dy,
 256                                           c->getZ ()+dz);
 257        }
 258
 259    for (int nro=0 ; nro<HE_MAXI; nro++)
 260        {
 261        int v1 = glob->EdgeVertex (nro, V_AMONT);
 262        int v2 = glob->EdgeVertex (nro, V_AVAL);
 263        i_edge[nro] = newEdge (i_node[v1], i_node[v2]);
 264
 265        if (db)
 266           {
 267           char nm0[8], nm1 [8], nm2 [8];
 268           printf (" %2d : %s = %s = [%s, %s] = [%d,%d] = [%s,%s]\n", nro,
 269                  glob->namofHexaEdge(nro), i_edge[nro]->getName(nm0),
 270                  glob->namofHexaVertex(v1), glob->namofHexaVertex(v2), v1, v2,
 271                  i_node[v1]->getName(nm1), i_node[v2]->getName(nm2));
 272           }
 273        }
 274
 275    for (int nro=0 ; nro<HQ_MAXI; nro++)
 276        i_quad[nro] = newQuad (i_edge[glob->QuadEdge (nro, E_A)],
 277                               i_edge[glob->QuadEdge (nro, E_B)],
 278                               i_edge[glob->QuadEdge (nro, E_C)],
 279                               i_edge[glob->QuadEdge (nro, E_D)]);
 280
 281    tab_hexa.push_back (newHexa (i_quad[Q_A], i_quad[Q_B], i_quad[Q_C],
 282                                 i_quad[Q_D], i_quad[Q_E], i_quad[Q_F]));
 283    double lambda  = 1;
 284    double dcell   = 1;
 285    for (int niv=0; niv<gr_rayon ; niv++)
 286        {
 287        double lambda0 = lambda;
 288        dcell  *= k;
 289        lambda += dcell;
 290        addStrate (i_quad, i_edge, i_node, c,  lambda/lambda0);
 291        }
 292
 293    return HOK;
 294 }
 295 // ====================================================== addStrate
 296 int Elements::addStrate (Quad* i_quad[], Edge* i_edge[], Vertex* i_node[],
 297                         Vertex* center,  double lambda)
 298 {
 299    Vertex* e_node [HV_MAXI];    // Les noeuds de l'hexa englobant
 300    Edge*   e_edge [HE_MAXI];    // Les noeuds de l'hexa englobant
 301    Quad*   e_quad [HQ_MAXI];    // Les noeuds de l'hexa englobant
 302
 303    Edge*   d_edge [HV_MAXI];    // Les aretes diagonales (1 par sommet)
 304    Quad*   d_quad [HE_MAXI];    // Les faces  diagonales (1 par arete)
 305
 306                                           // Les sommets
 307                                           //  + les aretes diagonales
 308    for (int nv=0 ; nv<HV_MAXI ; nv++)
 309        {
 310        double px0 = center->getX ();
 311        double py0 = center->getY ();
 312        double pz0 = center->getZ ();
 313        e_node[nv] = el_root->addVertex (px0+lambda*(i_node[nv]->getX()-px0),
 314                                         py0+lambda*(i_node[nv]->getY()-py0),
 315                                         pz0+lambda*(i_node[nv]->getZ()-pz0));
 316
 317        d_edge[nv] = newEdge (i_node[nv], e_node[nv]);
 318        }
 319                                           // Les aretes exterieures
 320                                           //  + les faces diagonales
 321    for (int nro=0 ; nro<HE_MAXI ; nro++)
 322        {
 323        int nv0  = glob->EdgeVertex (nro, V_AMONT);
 324        int nv1  = glob->EdgeVertex (nro, V_AVAL );
 325        e_edge[nro] = newEdge (e_node [nv0], e_node [nv1]);
 326        d_quad[nro] = newQuad (i_edge [nro], d_edge [nv0],
 327                               e_edge [nro], d_edge [nv1]);
 328        }
 329                                           // Les faces exterieures
 330                                           //  + les hexas
 331    Hexa* strate = NULL;
 332    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 333        {
 334        int ne0 = glob->QuadEdge (nro, E_A);
 335        int ne1 = glob->QuadEdge (nro, E_B);
 336        int ne2 = glob->QuadEdge (nro, E_C);
 337        int ne3 = glob->QuadEdge (nro, E_D);
 338
 339        e_quad[nro] = newQuad (e_edge[ne0], e_edge[ne1],
 340                               e_edge[ne2], e_edge[ne3]);
 341        strate  = newHexa (i_quad[nro], e_quad[nro], d_quad[ne0],
 342                           d_quad[ne2], d_quad[ne1], d_quad[ne3]);
 343        tab_hexa.push_back (strate);
 344        }
 345
 346    for (int nv=0 ; nv<HV_MAXI ; nv++) i_node  [nv] = e_node [nv];
 347    for (int ns=0 ; ns<HE_MAXI ; ns++) i_edge  [ns] = e_edge [ns];
 348    for (int nq=0 ; nq<HQ_MAXI ; nq++) i_quad  [nq] = e_quad [nq];
 349
 350    return HOK;
 351 }
 352 // ====================================================== saveVtk
 353 int Elements::saveVtk (cpchar nomfic)
 354 {
 355                                            // -- 1) Raz numerotation precedente
 356    for (int nro=0 ; nro<nbr_hexas ; nro++)
 357        {
 358        Hexa* cell = tab_hexa[nro];
 359        if (cell!=NULL && cell->isHere())
 360            cell->razNodes ();
 361        }
 362
 363    int nbnodes = 0;
 364    int nbcells = 0;
 365                                            // -- 2) Comptage
 366    for (int nro=0 ; nro<nbr_hexas ; nro++)
 367        {
 368        Hexa* cell = tab_hexa[nro];
 369        if (cell!=NULL && cell->isHere())
 370           {
 371           nbcells ++;
 372           nbnodes += cell->countNodes ();
 373           }
 374        }
 375
 376    pfile    vtk = fopen (nomfic, "w");
 377    fprintf (vtk, "# vtk DataFile Version 3.1\n");
 378    fprintf (vtk, "%s \n", nomfic);
 379    fprintf (vtk, "ASCII\n");
 380    fprintf (vtk, "DATASET UNSTRUCTURED_GRID\n");
 381    fprintf (vtk, "POINTS %d float\n", nbnodes);
 382
 383                                            // -- 2) Les noeuds
 384    int nronode = 0;
 385    for (int nro=0 ; nro<nbr_hexas ; nro++)
 386        {
 387        Hexa* cell = tab_hexa[nro];
 388        if (cell!=NULL && cell->isHere())
 389           cell->printNodes (vtk, nronode);
 390        }
 391                                            // -- 2) Les hexas
 392
 393    fprintf (vtk, "CELLS %d %d\n", nbcells, nbcells*(HV_MAXI+1));
 394
 395    for (int nro=0 ; nro<nbr_hexas ; nro++)
 396        {
 397        Hexa* cell = tab_hexa[nro];
 398        if (cell!=NULL && cell->isHere())
 399           cell->printHexa (vtk);
 400        }
 401
 402    fprintf (vtk, "CELL_TYPES %d\n", nbcells);
 403    for (int nro=0 ; nro<nbcells ; nro++)
 404        fprintf (vtk, "%d\n", HE_MAXI);
 405
 406    fclose (vtk);
 407    return HOK;
 408 }
 409 // ============ = = = = = = = = = Vertices ..........
 410 // ====================================================== newEdge
 411 Edge* Elements::newEdge (Vertex* v1, Vertex* v2)
 412 {
 413    if (v1==NULL || v2==NULL)
 414       return NULL;
 415
 416    Edge* elt = new Edge (v1, v2);
 417    return elt;
 418 }
 419 // ====================================================== newQuad
 420 Quad* Elements::newQuad (Edge* e1, Edge* e2, Edge* e3, Edge* e4)
 421 {
 422    if (e1==NULL || e2==NULL || e3==NULL|| e4==NULL)
 423       return NULL;
 424
 425    Quad* elt = new Quad (e1, e2, e3, e4);
 426    return elt;
 427 }
 428 // ====================================================== newHexa
 429 Hexa* Elements::newHexa (Quad* qa, Quad* qb, Quad* qc, Quad* qd,
 430                                              Quad* qe, Quad* qf)
 431 {
 432    if (qa==NULL || qb==NULL || qc==NULL|| qd==NULL || qe==NULL|| qf==NULL)
 433       return NULL;
 434
 435    Hexa*  elt = new Hexa (qa, qb, qc, qd, qe, qf);
 436    return elt;
 437 }
 438 // ====================================================== joinQuads
 439 int Elements::joinQuads (Quads& orig, int nb, Vertex* v1, Vertex* v2,
 440                   Vertex* v3, Vertex* v4, Quad* cible)
 441 {
 442    resize (GR_JOINT, orig.size(), nb);
 443
 444    el_root->markAll (IS_NONE);
 445    db = on_debug();
 446    // db = el_root->debug ();
 447
 448    gr_hauteur  = nb;
 449    nbr_orig = orig.size();
 450
 451    for (int nro=0 ; nro<nbr_orig ; nro++)
 452        {
 453        Quad* face = orig[nro];
 454        if (face==NULL)
 455           {
 456           printf ("\n");
 457           printf (" *** joinQuads : donnees incorrectes\n");
 458           printf (" *** le %deme quadrangle de depart est NULL\n", nro);
 459           setError ();
 460           return HERR;
 461           }
 462        else if (face->getNbrParents()>1)
 463           {
 464           printf ("\n");
 465           printf (" *** joinQuads : donnees incorrectes\n");
 466           printf (" *** le %deme quadrangle de depart n'est pas une "
 467                   "face externe\n", nro);
 468           face->dump ();
 469           setError ();
 470           return HERR;
 471           }
 472        orig [nro]->setMark (nro);
 473        tab_orig.push_back (orig[nro]);
 474        }
 475
 476    Edge* e_orig = tab_orig[0] -> findEdge (v1, v3);
 477    Edge* e_dest = cible       -> findEdge (v2, v4);
 478
 479    if (e_orig==NULL)
 480       {
 481       printf ("\n");
 482       printf (" *** joinQuads : donnees incorrectes\n");
 483       printf (" *** Les vertex v1 et v3 passes en argument ne definissent\n");
 484       printf (" *** pas une arete du 1er quadrangle de depart\n");
 485       printf ("\n");
 486       HexDump (v1);
 487       HexDump (v3);
 488       HexDump (orig[0]);
 489       }
 490
 491    if (e_dest==NULL)
 492       {
 493       printf ("\n");
 494       printf (" *** joinQuads : donnees incorrectes\n");
 495       printf (" *** Les vertex v2 et v4 passes en argument ne definissent\n");
 496       printf (" *** pas une arete du quadrangle cible\n");
 497       printf ("\n");
 498       HexDump (v2);
 499       HexDump (v4);
 500       HexDump (cible);
 501       }
 502
 503    if (e_orig==NULL || e_dest==NULL)
 504       {
 505       setError ();
 506       return HERR;
 507       }
 508
 509    StrOrient orient (v1, v3, v2, v4);
 510    int ier =  this   ->coupler (0, cible, &orient);
 511    if (ier!=HOK)
 512       {
 513       setError ();
 514       return HERR;
 515       }
 516    ier = orig[0]->coupler (cible, &orient, this);
 517    return ier;
 518 }
 519 // ======================================================== coupler
 520 int Elements::coupler (int nquad, Quad* dest, StrOrient* orient)
 521 {
 522    Quad* orig = tab_orig [nquad];
 523
 524    setQuad (orig, dir_z, 0, 0, 0);
 525    setQuad (dest, dir_z, nquad, 0, 0);
 526
 527    int n11 = orig->indexVertex (orient->v11);
 528    int n12 = orig->indexVertex (orient->v12);
 529    int n21 = dest->indexVertex (orient->v21);
 530    int n22 = dest->indexVertex (orient->v22);
 531
 532                // ---------------- Les 4 sommets initiaux
 533    Vertex* vorig[QUAD4] = { orient->v11, orient->v12,
 534                             orig->getVertex((n11+2) MODULO QUAD4),
 535                             orig->getVertex((n12+2) MODULO QUAD4) };
 536
 537    Vertex* vdest[QUAD4] = { orient->v21, orient->v22,
 538                             dest->getVertex((n21+2) MODULO QUAD4),
 539                             dest->getVertex((n22+2) MODULO QUAD4) };
 540    if (db)
 541       {
 542       printf ("Quad nro %d : ", nquad);
 543       orig->printName (" est couple avec ");
 544       dest->printName ("\n");
 545       printf ("Orientation : (");
 546       for (int ii=0 ; ii<QUAD4 ; ii++) printf("%s ", vorig[ii]->getName());
 547       printf (")\n");
 548       printf ("           -> (");
 549       for (int ii=0 ; ii<QUAD4 ; ii++) printf("%s ", vdest[ii]->getName());
 550       printf (")\n");
 551       }
 552
 553                // ---------------- Les sommets + les aretes verticales
 554    for (int ns=0 ; ns< QUAD4 ; ns++)
 555        {
 556        Vertex* nd1 = vorig[ns];
 557        Vertex* nd2 = vdest[ns];
 558        int nref0   = nd1->getMark();
 559        tab_vertex [nroVertex (nquad,ns,0)] = nd1;
 560
 561        if (nref0==IS_NONE)
 562           {
 563           double px0 = nd1->getX();
 564           double py0 = nd1->getY();
 565           double pz0 = nd1->getZ();
 566
 567           double dx = (nd2->getX() - nd1->getX()) / gr_hauteur;
 568           double dy = (nd2->getY() - nd1->getY()) / gr_hauteur;
 569           double dz = (nd2->getZ() - nd1->getZ()) / gr_hauteur;
 570
 571           nd1->setMark (indVertex (nquad, ns, 0));
 572
 573           Vertex* nd = nd1;
 574           for (int nh=0 ; nh<gr_hauteur ; nh++)
 575               {
 576               int nh1 = nh + 1;
 577               Vertex* ndp = nd;
 578               if (nh == gr_hauteur-1)
 579                  nd = nd2 ;
 580               else
 581                  nd = el_root->addVertex (px0 + nh1*dx, py0 + nh1*dy,
 582                                           pz0 + nh1*dz);
 583               int nv  = indVertex (nquad, ns, nh);
 584               tab_vertex [nv] = nd;
 585               tab_edge   [nv] = el_root->addEdge (ndp, nd);
 586               if (db)
 587                  printf (" Edge vertical nro %d = %s = (%s, %s)\n", nv,
 588                          tab_edge[nv]->getName(),
 589                          ndp->getName(), nd->getName());
 590               }
 591           }
 592        else
 593           {
 594           for (int nh=0 ; nh<gr_hauteur ; nh++)
 595               {
 596               int nv  = indVertex (nquad, ns, nh);
 597               int nv0 = indVertex (nref0, nh);
 598               tab_vertex [nv] = tab_vertex [nv0];
 599               tab_edge   [nv] = tab_edge   [nv0];
 600               }
 601           }
 602        }
 603                // ---------------- Les quads verticaux + aretes horiz
 604    for (int ns=0 ; ns< QUAD4 ; ns++)
 605        {
 606        int next = (ns+1) MODULO QUAD4;
 607        Edge*   arete = orig->findEdge (vorig[ns], vorig[next]);
 608        int     nref0 = arete->getMark();
 609        int     nref  = indVertex (nquad, ns, 0);
 610                                   // Construction des faces & aretes H
 611        if (nref0==IS_NONE)
 612           {
 613           arete->setMark (nref);
 614           Edge* ea = arete;
 615           Edge *eb, *ec, *ed;
 616           int  nva, nvb, nha;
 617
 618           for (int nh=0 ; nh< gr_hauteur ; nh++)
 619               {
 620               nva = indVertex (nquad, ns,   nh);
 621               nvb = indVertex (nquad, next, nh);
 622               nha = nroEdgeH  (nva);
 623
 624               ed  = tab_edge [nva];
 625               ec  = ea;
 626               eb  = tab_edge [nvb];
 627               if (nh==gr_hauteur-1)
 628                  ea = dest->findEdge (vdest[ns], vdest[next]);
 629               else
 630                  ea = el_root->addEdge (tab_vertex [nva], tab_vertex [nvb]);
 631
 632               propagateAssociation (ec, ea, eb);
 633               tab_quad [nva] = newQuad (ea, eb, ec, ed);
 634               if (BadElement (tab_quad [nva]))
 635                   return HERR;
 636               tab_edge [nha] = ea;
 637               }
 638           }
 639                                   // L'arete a deja ete traitee
 640        else
 641           {
 642           for (int nh=0 ; nh<gr_hauteur ; nh++)
 643               {
 644               int nva = indVertex (nquad, ns, nh);
 645               int nha = nroEdgeH  (nva);
 646
 647               int nvb = indVertex (nref0, nh);
 648               int nhb = nroEdgeH  (nvb);
 649
 650               tab_quad [nva] = tab_quad [nvb];
 651               tab_edge [nha] = tab_edge [nhb];
 652               }
 653           }
 654        }
 655                // -------------------------------- Les Hexas
 656    Quad *fb = orig;
 657    Quad *fa, *fc, *fd, *fe, *ff;
 658    for (int nh=0 ; nh< gr_hauteur ; nh++)
 659        {
 660        fa = fb;
 661        if (nh == gr_hauteur-1)
 662           fb = dest;
 663        else
 664           fb = newQuad (tab_edge [nroEdgeH (nquad, E_A, nh)],
 665                         tab_edge [nroEdgeH (nquad, E_B, nh)],
 666                         tab_edge [nroEdgeH (nquad, E_C, nh)],
 667                         tab_edge [nroEdgeH (nquad, E_D, nh)]);
 668
 669        fc = tab_quad [indVertex (nquad, E_A, nh)];
 670        fd = tab_quad [indVertex (nquad, E_C, nh)];
 671        fe = tab_quad [indVertex (nquad, E_B, nh)];
 672        ff = tab_quad [indVertex (nquad, E_D, nh)];
 673
 674        tab_hexa [nroHexa(nquad,nh)] = newHexa (fa,fb,fc,fd,fe,ff);
 675        if (BadElement (tab_hexa [nroHexa(nquad,nh)]))
 676            return HERR;
 677        }
 678    return HOK;
 679 }
 680 // ====================================================== makeCartesianNodes
 681 int Elements::makeCartesianNodes (Vertex* orig, Vector* v1, Vector* v2,
 682                    Vector* v3, int px, int py, int pz, int mx, int my, int mz)
 683 {
 684    double dx  = v1->getDx() + v2->getDx() + v3->getDx();
 685    double dy  = v1->getDy() + v2->getDy() + v3->getDy();
 686    double dz  = v1->getDz() + v2->getDz() + v3->getDz();
 687
 688    double px0 = orig->getX () - mx * dx;
 689    double py0 = orig->getY () - my * dy;
 690    double pz0 = orig->getZ () - mz * dz;
 691
 692    int nbre= 0;
 693
 694    for (int nz=0 ; nz<size_vz ; nz++)
 695        for (int ny=0 ; ny<size_vy ; ny++)
 696            for (int nx=0 ; nx<size_vx ; nx++)
 697                {
 698                Vertex* node = orig;
 699                if (nx!=mx || ny!=my || nz!=mz)
 700                   node = el_root->addVertex (px0 + nx*dx,  py0 + ny*dy,
 701                                                           pz0 + nz*dz);
 702                setVertex (node, nx, ny, nz);
 703                nbre++;
 704                }
 705    return HOK;
 706 }
 707 // ====================================================== makeCylindricalNodes
 708 int Elements::makeCylindricalNodes (Vertex* orig, Vector* base, Vector* haut,
 709             double dr, double da, double dl, int nr, int na, int nl, bool fill)
 710 {
 711    int ier = makeBasicCylinder (dr, da, dl, nr, na, nl, fill);
 712    if (ier!=HOK)
 713        return ier;
 714
 715    transfoVertices  (orig,  base, haut);
 716    return HOK;
 717 }
 718 // ====================================================== transfoVertices
 719 void Elements::transfoVertices (Vertex* orig, Vector* base, Vector* haut)
 720 {
 721    Vector* iprim = new Vector (base);
 722    Vector* jprim = new Vector (base);
 723    Vector* kprim = new Vector (haut);
 724
 725    int ier = kprim->renormer ();
 726    if (ier!=HOK)
 727        return;
 728
 729    jprim->vectoriel (kprim, base);
 730    ier = jprim->renormer ();
 731    if (ier!=HOK)
 732        return;
 733
 734    iprim->vectoriel (jprim, kprim);
 735    transfoVertices  (orig,  iprim, jprim, kprim);
 736 }
 737 // ====================================================== transfoVertices
 738 void Elements::transfoVertices (Vertex* orig, Vector* iprim, Vector* jprim,
 739                                 Vector* kprim)
 740 {
 741    double matrice[DIM3][DIM3]={{iprim->getDx(),jprim->getDx(),kprim->getDx()},
 742                                {iprim->getDy(),jprim->getDy(),kprim->getDy()},
 743                                {iprim->getDz(),jprim->getDz(),kprim->getDz()}};
 744
 745    double matkx = orig->getX();
 746    double matky = orig->getY();
 747    double matkz = orig->getZ();
 748
 749    int nbre = tab_vertex.size ();
 750    for (int nro=0 ; nro<nbre ; nro++)
 751        {
 752        if (tab_vertex[nro] != NULL)
 753            tab_vertex[nro]->setMark (NO_USED);
 754        }
 755
 756    for (int nro=0 ; nro<nbre ; nro++)
 757        {
 758        Vertex* node =  tab_vertex[nro];
 759        if (node != NULL  && node->getMark() == NO_USED)
 760           {
 761           double point [DIM3] = {node->getX(), node->getY(), node->getZ()};
 762           double result[DIM3] = {matkx, matky, matkz};
 763
 764           for (int ni=0 ; ni<DIM3; ni++)
 765               for (int nj=0 ; nj<DIM3; nj++)
 766                   result [ni] += matrice[ni][nj] * point[nj];
 767
 768           node->setCoord (result[dir_x], result[dir_y], result[dir_z]);
 769           node->setMark (IS_USED);
 770           }
 771        }
 772 }
 773 // ====================================================== transform
 774 int Elements::transform (Matrix* matrice)
 775 {
 776                                            // -- 1) Raz Marques
 777    for (int nro=0 ; nro<nbr_hexas ; nro++)
 778        {
 779        Hexa* cell = tab_hexa[nro];
 780        if (cell!=NULL && cell->isHere())
 781            cell->razNodes ();
 782        }
 783                                            // -- 2) Move Nodes
 784    for (int nro=0 ; nro<nbr_hexas ; nro++)
 785        {
 786        Hexa* cell = tab_hexa[nro];
 787        if (cell!=NULL && cell->isHere())
 788            cell->moveNodes (matrice);
 789        }
 790
 791    if (nbr_hexas!=0 )
 792       return HOK;
 793                     // -- Cas pathologique : il n'y avait pas d'hexas
 794                     // -- On se rabat sur les sommets
 795
 796    for (int nro=0 ; nro<nbr_vertex ; nro++)
 797        {
 798        Vertex* node = tab_vertex[nro];
 799        if (node!=NULL && node->isHere())
 800            matrice->perform (node);
 801        }
 802
 803    return HOK;
 804 }
 805 // ====================================================== cutHexas
 806 int Elements::cutHexas  (const Edges& t_edges, int nbcuts)
 807 {
 808                                        // 1) marquage des hexas
 809    el_root->markAll (NO_USED);
 810                                        // 2) Memo noeuds
 811    vector <Quad*> q_amont;
 812    vector <Quad*> q_aval;
 813    map    <Vertex*, Vertex*> vis_a_vis;
 814
 815    int nbnodes = t_edges.size();
 816    vector <Vertex*> v_amont (nbnodes);
 817    vector <Vertex*> v_aval  (nbnodes);
 818
 819    int nbfaces = 0;
 820    for (int nro=0; nro<nbnodes ; nro++)
 821        {
 822        Edge* arete   = t_edges [nro];
 823        v_amont [nro] = arete->getAmont ();
 824        v_aval  [nro] = arete->getAval  ();
 825        if (db)
 826           {
 827           printf (" %3d : Edge = (", nro);
 828           v_amont[nro]->printName (", ");
 829           v_aval [nro]->printName (")\n");
 830           }
 831
 832        vis_a_vis [v_amont[nro]] = v_aval[nro];
 833        vis_a_vis [v_aval[nro]] = v_amont[nro];
 834        int nbcells   = arete->getNbrParents ();
 835
 836        for (int nq=0 ; nq<nbcells ; nq++)
 837            {
 838            Quad* quad = arete->getParent (nq);
 839            if (quad->getMark () != IS_USED)
 840               {
 841               quad->setMark (IS_USED);
 842               int nbcubes = quad->getNbrParents ();
 843               for (int nh=0 ; nh<nbcubes ; nh++)
 844                   {
 845                   Hexa* hexa = quad->getParent (nh);
 846                   if (hexa->getMark () != IS_USED)
 847                      {
 848                      hexa->setMark (IS_USED);
 849                      int namont = hexa->getBase (v_amont[nro], arete);
 850                      int naval  = glob->getOpposedQuad (namont);
 851                      q_amont.push_back (hexa->getQuad  (namont));
 852                      q_aval .push_back (hexa->getQuad  (naval ));
 853
 854                      if (db)
 855                         {
 856                         printf (" %3d : Quad = ", nbfaces);
 857                         hexa->printName (", ");
 858                         printf (" Faces = (");
 859                         hexa->getQuad (namont)->printName (", ");
 860                         hexa->getQuad (naval )->printName (")\n");
 861                         nbfaces ++;
 862                         }
 863                      }
 864                   }
 865               }
 866            }
 867        }
 868                    // ------------------- Dimensionnement
 869    int nbcells   = q_amont.size ();
 870    nbr_vertex    = nbnodes*(nbcuts+2);
 871    int nbpiliers = nbnodes*(nbcuts+1);            // aretes verticales
 872    int nbpoutres = nbcells*(nbcuts+2)*QUAD4;      // aretes horizontales
 873    nbr_edges     = nbpoutres;
 874    nbr_quads     = nbcells*(nbcuts+1)*QUAD4;      // faces Verticales
 875    nbr_hexas     = nbcells*(nbcuts+1);
 876
 877                    // ------------------- Les noeuds et les aretes verticales
 878    tab_quad.resize   (nbr_quads);
 879    tab_edge.resize   (nbr_edges);
 880    tab_hexa.resize   (nbr_hexas);
 881    tab_vertex.resize (nbr_vertex);
 882    vector <Edge*>    tab_pilier (nbpiliers);
 883
 884    int nbinter = nbcuts + 1;
 885    for (int ned=0; ned<nbnodes ; ned++)
 886        {
 887        Shapes ass_shapes = t_edges[ned] -> getAssociations ();
 888        Edges  ass_edges;
 889        t_edges [ned]->remove ();
 890        Vertex* ndamont = v_amont [ned];
 891        Vertex* ndaval  = v_aval  [ned];
 892
 893        double dx = (ndaval->getX() - ndamont->getX()) / nbinter;
 894        double dy = (ndaval->getY() - ndamont->getY()) / nbinter;
 895        double dz = (ndaval->getZ() - ndamont->getZ()) / nbinter;
 896
 897        Vertex* nd0 = tab_vertex [ned] = ndamont;
 898        for (int nc=0; nc<nbcuts ; nc++)
 899            {
 900            int nc1 = nc+1;
 901            Vertex* nd1 = el_root->addVertex (ndamont->getX() + nc1*dx,
 902                                              ndamont->getY() + nc1*dy,
 903                                              ndamont->getZ() + nc1*dz);
 904            tab_vertex [nc1*nbnodes + ned] = nd1;
 905            tab_pilier [nc *nbnodes + ned] = newEdge (nd0, nd1);
 906            ass_edges.push_back (tab_pilier [nc *nbnodes + ned]);
 907            nd0 =  nd1;
 908            }
 909        tab_vertex [nbinter*nbnodes + ned] = ndaval;
 910        tab_pilier [nbcuts *nbnodes + ned] = newEdge (nd0, ndaval);
 911        ass_edges.push_back (tab_pilier[nbcuts *nbnodes + ned]);
 912        ndamont->setMark (ned);
 913        cutAssociation (ass_shapes, ass_edges);
 914        }
 915                    // ------------------- Les aretes horizontales
 916                    // ------------------- Les faces verticales
 917    HexDisplay (nbcells);
 918    int sizelig = nbcells*QUAD4;
 919    for (int nro=0; nro<nbcells ; nro++)
 920        {
 921        Quad* sol  = q_amont [nro];
 922        Quad* toit = q_aval  [nro];
 923        for (int ns=0; ns<QUAD4 ; ns++)
 924            {
 925            Edge* plinthe = sol->getEdge (ns);
 926            int   nmur  = nro*QUAD4 + ns;
 927            int   nmur0 = plinthe->getMark();
 928            if (nmur0 >= 0)
 929               {
 930               for (int nc=0 ; nc<nbinter ; nc++)
 931                   {
 932                   tab_edge [nc*sizelig + nmur] = tab_edge [nc*sizelig + nmur0];
 933                   tab_quad [nc*sizelig + nmur] = tab_quad [nc *sizelig + nmur0];
 934                   if (db)
 935                      {
 936                      printf (" %2d : %d quad_vertical [%02d] =", nro, ns,
 937                                                         nc*sizelig + nmur);
 938                      printf (" quad_vertical [%02d]\n",  nc*sizelig + nmur0);
 939                      }
 940                   }
 941               tab_edge [nbinter*sizelig+nmur] = tab_edge[nbinter*sizelig+nmur0];
 942               }
 943            else
 944               {
 945               plinthe->setMark (nmur);
 946               Vertex* vs1 = sol->getVertex (ns);
 947               Vertex* vs2 = sol->getVertex ((ns+1) MODULO QUAD4);
 948               int     nd1 = vs1->getMark ();
 949               int     nd2 = vs2->getMark ();
 950               Edge*   ed0 = tab_edge [nmur] = plinthe;
 951               Edge*   ed2 = NULL;
 952               Vertex* v1  = NULL;
 953               Vertex* v2  = NULL;
 954               for (int nc=0 ; nc<nbinter ; nc++)
 955                   {
 956                   int   nc1 = nc + 1;
 957                   if (nc<nbcuts)
 958                      {
 959                      v1 = tab_vertex [nc1*nbnodes + nd1];
 960                      v2 = tab_vertex [nc1*nbnodes + nd2];
 961                      ed2 = newEdge (v1, v2);
 962                      }
 963                   else
 964                      {
 965                      v1 = vis_a_vis [vs1];
 966                      v2 = vis_a_vis [vs2];
 967                      ed2 = toit->findEdge (v1, v2);
 968                      }
 969
 970                   tab_edge [nc1*sizelig + nmur] = ed2;
 971                   tab_quad [nc *sizelig + nmur] = newQuad (ed0,
 972                             tab_pilier [nc*nbnodes + nd1], ed2,
 973                             tab_pilier [nc*nbnodes + nd2]);
 974                   ed0 = ed2;
 975                   if (db)
 976                      {
 977                      printf (" %2d : %d quad_vertical [%02d] = ", nro, ns,
 978                                                          nc*sizelig + nmur);
 979                      PrintName (tab_quad [nc *sizelig + nmur]);
 980                      printf ("\n");
 981                      }
 982                   }
 983               }
 984            }
 985        }
 986                    // ------------------- Les faces horizontales
 987                    // ------------------- Les hexas
 988    // Rappel : sizelig = nbcells*QUAD4
 989    for (int nro=0; nro<nbcells ; nro++)
 990        {
 991        Quad* qa = q_amont [nro];
 992        for (int nc=0 ; nc<nbinter ; nc++)
 993            {
 994            int  quadv =  nc*nbcells*QUAD4 + nro*QUAD4;
 995            Quad* qb = q_aval[nro];
 996            if (nc<nbcuts)
 997               {
 998               int  edh   = (nc+1)*nbcells*QUAD4 + nro*QUAD4;
 999               qb = newQuad (tab_edge[edh],   tab_edge[edh+1],
1000                             tab_edge[edh+2], tab_edge[edh+3]);
1001               if (BadElement(qb))
1002                   return HERR;
1003               }
1004            Quad* qc = tab_quad [quadv];
1005            Quad* qd = tab_quad [quadv + 2];
1006            Quad* qe = tab_quad [quadv + 1];
1007            Quad* qf = tab_quad [quadv + 3];
1008            Hexa* hexa = newHexa (qa, qb, qc, qd, qe, qf);
1009            if (BadElement(hexa))
1010                return HERR;
1011            tab_hexa [nc*nbcells + nro] = hexa;
1012            qa = qb;
1013            }
1014        }
1015    return HOK;
1016 }
1017 // ====================================================== clear_associations
1018 void clear_associations (std::vector<Quad*>& table)
1019 {
1020    int nbelts = table.size();
1021    for (int nro=0 ; nro<nbelts ; nro++)
1022        {
1023        Quad* elt = table [nro];
1024        if (elt != NULL && elt->isValid())
1025            elt->clearAssociation ();
1026        }
1027 }
1028 // ====================================================== clear_associations
1029 void clear_associations (std::vector<Edge*>& table)
1030 {
1031    int nbelts = table.size();
1032    for (int nro=0 ; nro<nbelts ; nro++)
1033        {
1034        Edge* elt = table [nro];
1035        if (elt != NULL && elt->isValid())
1036            elt->clearAssociation ();
1037        }
1038 }
1039 // ====================================================== clear_associations
1040 void clear_associations (std::vector<Vertex*>& table)
1041 {
1042    int nbelts = table.size();
1043    for (int nro=0 ; nro<nbelts ; nro++)
1044        {
1045        Vertex* elt = table [nro];
1046        if (elt != NULL && elt->isValid())
1047            elt->clearAssociation ();
1048        }
1049 }
1050 // ====================================================== clearAssociation
1051 void Elements::clearAssociation  ()
1052 {
1053    clear_associations (tab_quad);
1054    clear_associations (tab_edge);
1055    clear_associations (tab_vertex);
1056
1057
1058    clear_associations (ker_hquad);
1059    clear_associations (ker_vquad);
1060    clear_associations (ker_hedge);
1061    clear_associations (ker_vedge);
1062 }
1063 // ============================================================ findVertex
1064 int Elements::findVertex (double vx, double vy, double vz)
1065 {
1066    double tol = el_root->getTolerance ();
1067    double xmin = vx - tol;
1068    double xmax = vx + tol;
1069    double ymin = vy - tol;
1070    double ymax = vy + tol;
1071    double zmin = vz - tol;
1072    double zmax = vz + tol;
1073
1074    int nbre = tab_vertex.size();
1075    for (int nro=0 ; nro<nbre ; nro++)
1076        {
1077        Vertex* node = tab_vertex [nro];
1078        if (node != NULL && node->isHere ()
1079                         && node->isin   (xmin, xmax, ymin, ymax, zmin, zmax))
1080               return nro;
1081        }
1082    return NOTHING;
1083 }
1084
1085 END_NAMESPACE_HEXA