src/HEXABLOCK/HexElements.cxx

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