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