src/HEXABLOCK/HexElements.cxx

   1
   2 // C++ : Table des noeuds
   3
   4 #include "HexElements.hxx"
   5
   6 #include "HexDocument.hxx"
   7 #include "HexVector.hxx"
   8 #include "HexVertex.hxx"
   9 #include "HexHexa.hxx"
  10 #include "HexEdge.hxx"
  11
  12 #include "HexGlobale.hxx"
  13
  14 #include <cmath>
  15
  16 static bool db=true;
  17
  18 BEGIN_NAMESPACE_HEXA
  19
  20 // ====================================================== Constructeur
  21 Elements::Elements (Document* doc) : EltBase (doc)
  22 {
  23    glob  = Globale::getInstance ();
  24
  25    grid_type  = GR_NONE;
  26    size_qx = size_ex = size_vx = size_hx = 0;
  27    size_qy = size_ey = size_vy = size_hy = 0;
  28    size_qz = size_ez = size_vz = size_hz = 0;
  29    size_qplus = size_eplus = size_vplus = size_hplus = 0;
  30
  31    nbr_hexas  = nbr_quads  = nbr_edges  = nbr_vertex = 0;
  32    ker_vertex = 0;
  33    cyl_closed = false;
  34    cyl_fill   = false;
  35    cyl_dispo  = CYL_NOFILL;
  36 }
  37 // ====================================================== Constructeur
  38 Elements::Elements (Document* doc, int nx, int ny, int nz) : EltBase (doc)
  39 {
  40    glob  = Globale::getInstance ();
  41
  42    grid_type  = GR_NONE;
  43    size_qx = size_ex = size_vx = size_hx = 0;
  44    size_qy = size_ey = size_vy = size_hy = 0;
  45    size_qz = size_ez = size_vz = size_hz = 0;
  46    size_qplus = size_eplus = size_vplus = size_hplus = 0;
  47
  48    nbr_hexas  = nbr_quads  = nbr_edges  = nbr_vertex = 0;
  49    cyl_closed = true;
  50    cyl_fill   = false;
  51    cyl_dispo  = CYL_NOFILL;
  52
  53    resize (GR_CYLINDRIC, nx, ny, nz);
  54    cyl_closed = true;
  55 }
  56 // ====================================================== Constructeur (clonage)
  57 Elements::Elements (Elements* orig) : EltBase (orig->el_root)
  58 {
  59    glob = Globale::getInstance ();
  60
  61    grid_type  = orig->grid_type;
  62    cyl_closed = orig->cyl_closed;
  63    cyl_fill   = orig->cyl_fill;
  64    cyl_dispo  = orig->cyl_dispo;
  65
  66    size_qx = size_ex = size_vx = size_hx = 0;
  67    size_qy = size_ey = size_vy = size_hy = 0;
  68    size_qz = size_ez = size_vz = size_hz = 0;
  69    size_qplus = size_eplus = size_vplus = size_hplus = 0;
  70    nbr_hexas  = nbr_quads  = nbr_edges  = nbr_vertex = 0;
  71
  72    resize (orig->grid_type, orig->size_hx, orig->size_hy, orig->size_hz);
  73    cyl_closed = orig->cyl_closed;
  74 }
  75 // ====================================================== resize
  76 void Elements::resize (EnumGrid type, int nx, int ny, int nz)
  77 {
  78    grid_type = type;
  79
  80    switch (grid_type)
  81       {
  82       case GR_CARTESIAN :
  83       case GR_CYLINDRIC :
  84       default :
  85            size_hx = std::max (nx, 1);
  86            size_hy = std::max (ny, 1);
  87            size_hz = std::max (nz, 1);
  88
  89            size_qx = size_ex = size_vx = size_hx + 1;
  90            size_qy = size_ey = size_vy = size_hy + 1;
  91            size_qz = size_ez = size_vz = size_hz + 1;
  92
  93            nbr_hexas  = size_hx * size_hy * size_hz;
  94            nbr_quads  = size_qx * size_qy * size_qz * DIM3;
  95            nbr_edges  = size_ex * size_ey * size_ez * DIM3;
  96            nbr_vertex = size_vx * size_vy * size_vz;
  97            break;
  98
  99       case GR_SPHERIC :
 100            size_qx = size_ex = size_vx = size_hx = 0;
 101            size_qy = size_ey = size_vy = size_hy = 0;
 102            size_qz = size_ez = size_vz = size_hz = 0;
 103            nbr_quads  = nbr_edges  = nbr_vertex = 0;
 104            gr_rayon = std::max (nx, 1);
 105
 106            nbr_hexas = 1 +  gr_rayon*HQ_MAXI;
 107            tab_hexa.clear ();
 108            ker_vertex = nbr_vertex;
 109            return;
 110
 111       case GR_JOINT :
 112            nbr_orig = std::max (nx, 1);
 113            gr_hauteur  = std::max (ny, 1);
 114            size_hx  = nbr_orig;
 115            size_hy  = 1;
 116            size_hz  = gr_hauteur;
 117
 118            nbr_hexas  = nbr_orig * gr_hauteur;
 119            nbr_vertex = nbr_hexas * QUAD4;
 120            nbr_vertex = nbr_orig * (gr_hauteur+1)*QUAD4;
 121            nbr_quads  = nbr_vertex;
 122            nbr_edges  = 2*nbr_vertex;
 123            break;
 124
 125       case GR_BICYL :
 126            cyl_closed = true;
 127            size_hx = 1;
 128            size_hy = 8;
 129            size_hz = 4;
 130
 131            size_qx = size_ex = size_vx = size_hx + 1;
 132            size_qy = size_ey = size_vy = size_hy + 1;
 133            size_qz = size_ez = size_vz = size_hz + 1;
 134
 135            nbr_hexas  = size_hx * size_hy * size_hz;
 136            nbr_quads  = size_qx * size_qy * size_qz * DIM3;
 137            nbr_edges  = size_ex * size_ey * size_ez * DIM3;
 138            nbr_vertex = size_vx * size_vy * size_vz;
 139            break;
 140       }
 141
 142    tab_hexa  .resize (nbr_hexas);
 143    tab_quad  .resize (nbr_quads);
 144    tab_edge  .resize (nbr_edges);
 145    tab_vertex.resize (nbr_vertex);
 146
 147    ker_vertex = nbr_vertex;
 148
 149    for (int nc=0 ; nc< nbr_hexas  ; nc++) tab_hexa   [nc] = NULL;
 150    for (int nc=0 ; nc< nbr_quads  ; nc++) tab_quad   [nc] = NULL;
 151    for (int nc=0 ; nc< nbr_edges  ; nc++) tab_edge   [nc] = NULL;
 152    for (int nc=0 ; nc< nbr_vertex ; nc++) tab_vertex [nc] = NULL;
 153 }
 154
 155 // ====================================================== makeCartesianGrid
 156 int Elements::makeCartesianGrid (Vertex* orig, Vector* v1, Vector* v2,
 157                    Vector* v3, int px, int py, int pz, int mx, int my, int mz)
 158 {
 159    resize (GR_CARTESIAN, px+mx, py+my, pz+mz);
 160
 161    makeCartesianNodes (orig, v1, v2, v3, px, py, pz, mx, my, mz);
 162
 163    fillGrid ();
 164    return HOK;
 165 }
 166 // ====================================================== makeCylindricalGrid
 167 int Elements::makeCylindricalGrid (Vertex* c, Vector* b, Vector* h,
 168          double dr, double da, double dl, int nr, int na, int nl, bool fill)
 169 {
 170    resize (GR_CYLINDRIC, nr, na, nl);
 171    cyl_closed = da >= 360.0;
 172    makeCylindricalNodes (c, b, h, dr, da, dl, nr, na, nl, fill);
 173    fillGrid ();
 174    return HOK;
 175 }
 176 // ====================================================== makeSphericalGrid
 177 int Elements::makeSphericalGrid (Vertex* c, Vector* dv, int nb, double  k)
 178 {
 179    resize (GR_SPHERIC, nb);
 180
 181    if (nb<0)
 182       return HERR;
 183    else if (dv->getDx()<=ZEROR || dv->getDy()<=ZEROR || dv->getDz()<=ZEROR)
 184       return HERR;
 185
 186    Vertex* i_node [HV_MAXI];    // Les noeuds de l'hexa englobant
 187    Edge*   i_edge [HE_MAXI];    // Les noeuds de l'hexa englobant
 188    Quad*   i_quad [HQ_MAXI];    // Les noeuds de l'hexa englobant
 189
 190    for (int nro=0 ; nro<HV_MAXI; nro++)
 191        {
 192        double dx = glob->CoordVertex (nro, dir_x) * dv->getDx();
 193        double dy = glob->CoordVertex (nro, dir_y) * dv->getDy();
 194        double dz = glob->CoordVertex (nro, dir_z) * dv->getDz();
 195
 196        i_node [nro] = el_root->addVertex (c->getX ()+dx, c->getY ()+dy,
 197                                           c->getZ ()+dz);
 198        }
 199
 200    for (int nro=0 ; nro<HE_MAXI; nro++)
 201        i_edge[nro] = newEdge (i_node[glob->EdgeVertex (nro, V_AMONT)],
 202                               i_node[glob->EdgeVertex (nro, V_AVAL) ]);
 203
 204    for (int nro=0 ; nro<HQ_MAXI; nro++)
 205        i_quad[nro] = newQuad (i_edge[glob->QuadEdge (nro, E_A)],
 206                               i_edge[glob->QuadEdge (nro, E_B)],
 207                               i_edge[glob->QuadEdge (nro, E_C)],
 208                               i_edge[glob->QuadEdge (nro, E_D)]);
 209
 210    tab_hexa.push_back (newHexa (i_quad[Q_A], i_quad[Q_B], i_quad[Q_C],
 211                                 i_quad[Q_D], i_quad[Q_E], i_quad[Q_F]));
 212    double lambda  = 1;
 213    double dcell   = 1;
 214    for (int niv=0; niv<gr_rayon ; niv++)
 215        {
 216        double lambda0 = lambda;
 217        dcell  *= k;
 218        lambda += dcell;
 219        addStrate (i_quad, i_edge, i_node, c,  lambda/lambda0);
 220        }
 221
 222    return HOK;
 223 }
 224 // ====================================================== addStrate
 225 int Elements::addStrate (Quad* i_quad[], Edge* i_edge[], Vertex* i_node[],
 226                         Vertex* center,  double lambda)
 227 {
 228    Vertex* e_node [HV_MAXI];    // Les noeuds de l'hexa englobant
 229    Edge*   e_edge [HE_MAXI];    // Les noeuds de l'hexa englobant
 230    Quad*   e_quad [HQ_MAXI];    // Les noeuds de l'hexa englobant
 231
 232    Edge*   d_edge [HV_MAXI];    // Les aretes diagonales (1 par sommet)
 233    Quad*   d_quad [HE_MAXI];    // Les faces  diagonales (1 par arete)
 234
 235                                           // Les sommets
 236                                           //  + les aretes diagonales
 237    for (int nv=0 ; nv<HV_MAXI ; nv++)
 238        {
 239        double px0 = center->getX ();
 240        double py0 = center->getY ();
 241        double pz0 = center->getZ ();
 242        e_node[nv] = el_root->addVertex (px0+lambda*(i_node[nv]->getX()-px0),
 243                                         py0+lambda*(i_node[nv]->getY()-py0),
 244                                         pz0+lambda*(i_node[nv]->getZ()-pz0));
 245
 246        d_edge[nv] = newEdge (i_node[nv], e_node[nv]);
 247        }
 248                                           // Les aretes exterieures
 249                                           //  + les faces diagonales
 250    for (int nro=0 ; nro<HE_MAXI ; nro++)
 251        {
 252        int nv0  = glob->EdgeVertex (nro, V_AMONT);
 253        int nv1  = glob->EdgeVertex (nro, V_AVAL );
 254        e_edge[nro] = newEdge (e_node [nv0], e_node [nv1]);
 255        d_quad[nro] = newQuad (i_edge [nro], d_edge [nv0],
 256                               e_edge [nro], d_edge [nv1]);
 257        }
 258                                           // Les faces exterieures
 259                                           //  + les hexas
 260    Hexa* strate = NULL;
 261    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 262        {
 263        int ne0 = glob->QuadEdge (nro, E_A);
 264        int ne1 = glob->QuadEdge (nro, E_B);
 265        int ne2 = glob->QuadEdge (nro, E_C);
 266        int ne3 = glob->QuadEdge (nro, E_D);
 267
 268        e_quad[nro] = newQuad (e_edge[ne0], e_edge[ne1],
 269                               e_edge[ne2], e_edge[ne3]);
 270        strate  = newHexa (i_quad[nro], e_quad[nro], d_quad[ne0],
 271                           d_quad[ne2], d_quad[ne1], d_quad[ne3]);
 272        tab_hexa.push_back (strate);
 273        }
 274
 275    for (int nv=0 ; nv<HV_MAXI ; nv++) i_node  [nv] = e_node [nv];
 276    for (int ns=0 ; ns<HE_MAXI ; ns++) i_edge  [ns] = e_edge [ns];
 277    for (int nq=0 ; nq<HE_MAXI ; nq++) i_quad  [nq] = e_quad [nq];
 278
 279    return HOK;
 280 }
 281 // ====================================================== saveVtk
 282 int Elements::saveVtk (cpchar nomfic)
 283 {
 284                                            // -- 1) Raz numerotation precedente
 285    for (int nro=0 ; nro<nbr_hexas ; nro++)
 286        {
 287        Hexa* cell = tab_hexa[nro];
 288        if (cell!=NULL && cell->isHere())
 289            cell->razNodes ();
 290        }
 291
 292    int nbnodes = 0;
 293    int nbcells = 0;
 294                                            // -- 2) Comptage
 295    for (int nro=0 ; nro<nbr_hexas ; nro++)
 296        {
 297        Hexa* cell = tab_hexa[nro];
 298        if (cell!=NULL && cell->isHere())
 299           {
 300           nbcells ++;
 301           nbnodes += cell->countNodes ();
 302           }
 303        }
 304
 305    pfile    vtk = fopen (nomfic, "w");
 306    fprintf (vtk, "# vtk DataFile Version 3.1\n");
 307    fprintf (vtk, "%s \n", nomfic);
 308    fprintf (vtk, "ASCII\n");
 309    fprintf (vtk, "DATASET UNSTRUCTURED_GRID\n");
 310    fprintf (vtk, "POINTS %d float\n", nbnodes);
 311
 312                                            // -- 2) Les noeuds
 313    int nronode = 0;
 314    for (int nro=0 ; nro<nbr_hexas ; nro++)
 315        {
 316        Hexa* cell = tab_hexa[nro];
 317        if (cell!=NULL && cell->isHere())
 318           cell->printNodes (vtk, nronode);
 319        }
 320                                            // -- 2) Les hexas
 321
 322    fprintf (vtk, "CELLS %d %d\n", nbcells, nbcells*(HV_MAXI+1));
 323
 324    for (int nro=0 ; nro<nbr_hexas ; nro++)
 325        {
 326        Hexa* cell = tab_hexa[nro];
 327        if (cell!=NULL && cell->isHere())
 328           cell->printHexa (vtk);
 329        }
 330
 331    fprintf (vtk, "CELL_TYPES %d\n", nbcells);
 332    for (int nro=0 ; nro<nbcells ; nro++)
 333        fprintf (vtk, "%d\n", HE_MAXI);
 334
 335    fclose (vtk);
 336    return HOK;
 337 }
 338 // ============ = = = = = = = = = Vertices ..........
 339 // ====================================================== newEdge
 340 Edge* Elements::newEdge (Vertex* v1, Vertex* v2)
 341 {
 342    if (v1==NULL || v2==NULL)
 343       return NULL;
 344
 345    Edge* elt = new Edge (v1, v2);
 346    return elt;
 347 }
 348 // ====================================================== newQuad
 349 Quad* Elements::newQuad (Edge* e1, Edge* e2, Edge* e3, Edge* e4)
 350 {
 351    if (e1==NULL || e2==NULL || e3==NULL|| e4==NULL)
 352       return NULL;
 353
 354    Quad* elt = new Quad (e1, e2, e3, e4);
 355    return elt;
 356 }
 357 // ====================================================== newHexa
 358 Hexa* Elements::newHexa (Quad* qa, Quad* qb, Quad* qc, Quad* qd,
 359                                              Quad* qe, Quad* qf)
 360 {
 361    if (qa==NULL || qb==NULL || qc==NULL|| qd==NULL || qe==NULL|| qf==NULL)
 362       return NULL;
 363
 364    Hexa*  elt = new Hexa (qa, qb, qc, qd, qe, qf);
 365    return elt;
 366 }
 367 // ====================================================== joinQuads
 368 int Elements::joinQuads (Quads& orig, int nb, Vertex* v1, Vertex* v2,
 369                   Vertex* v3, Vertex* v4, Quad* cible)
 370 {
 371    resize (GR_JOINT, orig.size(), nb);
 372
 373    el_root->markAll (IS_NONE);
 374
 375    gr_hauteur  = nb;
 376    nbr_orig = orig.size();
 377
 378    for (int nro=0 ; nro<nbr_orig ; nro++)
 379        {
 380        Quad* face = orig[nro];
 381        if (face==NULL)
 382           {
 383           printf ("\n");
 384           printf (" *** joinQuads : donnees incorrectes\n");
 385           printf (" *** le %deme quadrangle de depart est NULL\n", nro);
 386           return HERR;
 387           }
 388        else if (face->getNbrParents()>1)
 389           {
 390           printf ("\n");
 391           printf (" *** joinQuads : donnees incorrectes\n");
 392           printf (" *** le %deme quadrangle de depart n'est pas une "
 393                   "face externe\n", nro);
 394           face->dump ();
 395           return HERR;
 396           }
 397        orig [nro]->setMark (nro);
 398        tab_orig.push_back (orig[nro]);
 399        }
 400
 401    Edge* e_orig = tab_orig[0] -> findEdge (v1, v3);
 402    Edge* e_dest = cible        -> findEdge (v2, v4);
 403    HexDump (e_orig);
 404    HexDump (e_dest);
 405
 406    if (e_orig==NULL)
 407       {
 408       printf ("\n");
 409       printf (" *** joinQuads : donnees incorrectes\n");
 410       printf (" *** Les vertex v1 et v3 passes en argument ne definissent\n");
 411       printf (" *** pas une arete du 1er quadrangle de depart\n");
 412       printf ("\n");
 413       HexDump (v1);
 414       HexDump (v3);
 415       HexDump (orig[0]);
 416       }
 417
 418    if (e_dest==NULL)
 419       {
 420       printf ("\n");
 421       printf (" *** joinQuads : donnees incorrectes\n");
 422       printf (" *** Les vertex v2 et v4 passes en argument ne definissent\n");
 423       printf (" *** pas une arete du quadrangle cible\n");
 424       printf ("\n");
 425       HexDump (v2);
 426       HexDump (v4);
 427       HexDump (cible);
 428       }
 429
 430    if (e_orig==NULL || e_dest==NULL)
 431       return HERR;
 432
 433    StrOrient orient (v1, v3, v2, v4);
 434    this   ->coupler (0, cible, &orient);
 435    orig[0]->coupler (cible, &orient, this);
 436    return HOK;
 437 }
 438 // ======================================================== coupler
 439 void Elements::coupler (int nquad, Quad* dest, StrOrient* orient)
 440 {
 441    Quad* orig = tab_orig [nquad];
 442
 443    setQuad (orig, dir_z, 0, 0, 0);
 444    setQuad (dest, dir_z, nquad, 0, 0);
 445
 446    int n11 = orig->indexVertex (orient->v11);
 447    int n12 = orig->indexVertex (orient->v12);
 448    int n21 = dest->indexVertex (orient->v21);
 449    int n22 = dest->indexVertex (orient->v22);
 450
 451    printf ("Quad nro %d : ", nquad);
 452    orig->printName (" est couple avec ");
 453    dest->printName (", ");
 454    printf ("Orientation=(%d,%d) (%d,%d)\n", n11, n12, n21, n22);
 455                // ---------------- Les 4 sommets initiaux
 456    Vertex* vorig[QUAD4] = { orient->v11, orient->v12,
 457                             orig->getVertex((n11+2) MODULO QUAD4),
 458                             orig->getVertex((n12+2) MODULO QUAD4) };
 459
 460    Vertex* vdest[QUAD4] = { orient->v21, orient->v22,
 461                             dest->getVertex((n21+2) MODULO QUAD4),
 462                             dest->getVertex((n22+2) MODULO QUAD4) };
 463
 464                // ---------------- Les sommets + les aretes verticales
 465 #define Put(e) if (e) { printf ( " ... " #e " = ") ; e->printName("\n"); } \
 466    else printf ( " ... " #e " = 0x0\n")
 467 #define Putn(e,n) { printf ( " ... " #e "[%d] = ",n) ; if(e[n]) e[n]->printName("\n"); \
 468    else printf ( "0x0\n")
 469
 470    for (int ns=0 ; ns< QUAD4 ; ns++)
 471        {
 472        Vertex* nd1 = vorig[ns];
 473        Vertex* nd2 = vdest[ns];
 474        int nref0   = nd1->getMark();
 475        tab_vertex [nroVertex (nquad,ns,0)] = nd1;
 476
 477        if (nref0==IS_NONE)
 478           {
 479           double px0 = nd1->getX();
 480           double py0 = nd1->getY();
 481           double pz0 = nd1->getZ();
 482
 483           double dx = (nd2->getX() - nd1->getX()) / gr_hauteur;
 484           double dy = (nd2->getY() - nd1->getY()) / gr_hauteur;
 485           double dz = (nd2->getZ() - nd1->getZ()) / gr_hauteur;
 486
 487           nd1->setMark (indVertex (nquad, ns, 0));
 488
 489           Vertex* nd = nd1;
 490           for (int nh=0 ; nh<gr_hauteur ; nh++)
 491               {
 492               int nh1 = nh + 1;
 493               Vertex* ndp = nd;
 494               if (nh == gr_hauteur-1)
 495                  nd = nd2 ;
 496               else
 497                  nd = el_root->addVertex (px0 + nh1*dx, py0 + nh1*dy,
 498                                           pz0 + nh1*dz);
 499               int nv  = indVertex (nquad, ns, nh);
 500               tab_vertex [nv] = nd;
 501               tab_edge   [nv] = el_root->addEdge (ndp, nd);
 502               }
 503           }
 504        else
 505           {
 506           for (int nh=0 ; nh<gr_hauteur ; nh++)
 507               {
 508               int nv  = indVertex (nquad, ns, nh);
 509               int nv0 = indVertex (nref0, nh);
 510               tab_vertex [nv] = tab_vertex [nv0];
 511               tab_edge   [nv] = tab_edge   [nv0];
 512               }
 513           }
 514        }
 515                // ---------------- Les quads verticaux + aretes horiz
 516    for (int ns=0 ; ns< QUAD4 ; ns++)
 517        {
 518        int next = (ns+1) MODULO QUAD4;
 519        Edge*   arete = orig->findEdge (vorig[ns], vorig[next]);
 520        int     nref0 = arete->getMark();
 521        int     nref  = indVertex (nquad, ns, 0);
 522                                   // Construction des faces & aretes H
 523        if (nref0==IS_NONE)
 524           {
 525           arete->setMark (nref);
 526           Edge* ea = arete;
 527           Edge *eb, *ec, *ed;
 528           int  nva, nvb, nha;
 529
 530           for (int nh=0 ; nh< gr_hauteur ; nh++)
 531               {
 532               nva = indVertex (nquad, ns,   nh);
 533               nvb = indVertex (nquad, next, nh);
 534               nha = nroEdgeH  (nva);
 535
 536               ed  = tab_edge [nva];
 537               ec  = ea;
 538               eb  = tab_edge [nvb];
 539               if (nh==gr_hauteur-1)
 540                  ea = dest->findEdge (vdest[ns], vdest[next]);
 541               else
 542                  ea = el_root->addEdge (tab_vertex [nva], tab_vertex [nvb]);
 543
 544               tab_quad [nva] = newQuad (ea, eb, ec, ed);
 545               tab_edge [nha] = ea;
 546               }
 547           }
 548                                   // L'arete a deja ete traitee
 549        else
 550           {
 551           for (int nh=0 ; nh<gr_hauteur ; nh++)
 552               {
 553               int nva = indVertex (nquad, ns, nh);
 554               int nha = nroEdgeH  (nva);
 555
 556               int nvb = indVertex (nref0, nh);
 557               int nhb = nroEdgeH  (nvb);
 558
 559               tab_quad [nva] = tab_quad [nvb];
 560               tab_edge [nha] = tab_edge [nhb];
 561               }
 562           }
 563        }
 564                // -------------------------------- Les Hexas
 565    Quad *fb = orig;
 566    Quad *fa, *fc, *fd, *fe, *ff;
 567    for (int nh=0 ; nh< gr_hauteur ; nh++)
 568        {
 569        fa = fb;
 570        if (nh == gr_hauteur-1)
 571           fb = dest;
 572        else
 573           fb = newQuad (tab_edge [nroEdgeH (nquad, E_A, nh)],
 574                         tab_edge [nroEdgeH (nquad, E_B, nh)],
 575                         tab_edge [nroEdgeH (nquad, E_C, nh)],
 576                         tab_edge [nroEdgeH (nquad, E_D, nh)]);
 577
 578        fc = tab_quad [indVertex (nquad, E_A, nh)];
 579        fd = tab_quad [indVertex (nquad, E_C, nh)];
 580        fe = tab_quad [indVertex (nquad, E_B, nh)];
 581        ff = tab_quad [indVertex (nquad, E_D, nh)];
 582
 583        tab_hexa [nroHexa(nquad,nh)] = newHexa (fa,fb,fc,fd,fe,ff);
 584        }
 585 }
 586 // ====================================================== makeCartesianNodes
 587 int Elements::makeCartesianNodes (Vertex* orig, Vector* v1, Vector* v2,
 588                    Vector* v3, int px, int py, int pz, int mx, int my, int mz)
 589 {
 590    double dx  = v1->getDx() + v2->getDx() + v3->getDx();
 591    double dy  = v1->getDy() + v2->getDy() + v3->getDy();
 592    double dz  = v1->getDz() + v2->getDz() + v3->getDz();
 593
 594    double px0 = orig->getX () - mx * dx;
 595    double py0 = orig->getY () - my * dy;
 596    double pz0 = orig->getZ () - mz * dz;
 597
 598    int nbre= 0;
 599
 600    for (int nz=0 ; nz<size_vz ; nz++)
 601        for (int ny=0 ; ny<size_vy ; ny++)
 602            for (int nx=0 ; nx<size_vx ; nx++)
 603                {
 604                Vertex* node = orig;
 605                if (nx!=mx || ny!=my || nz!=mz)
 606                   node = el_root->addVertex (px0 + nx*dx,  py0 + ny*dy,
 607                                                           pz0 + nz*dz);
 608                setVertex (node, nx, ny, nz);
 609                nbre++;
 610                }
 611    return HOK;
 612 }
 613 // ====================================================== makeCylindricalNodes
 614 int Elements::makeCylindricalNodes (Vertex* orig, Vector* base, Vector* haut,
 615             double dr, double da, double dl, int nr, int na, int nl, bool fill)
 616 {
 617    int ier = makeBasicCylinder (dr, da, dl, nr, na, nl, fill);
 618    if (ier!=HOK)
 619        return ier;
 620
 621    Vector* iprim = new Vector (base);
 622    Vector* jprim = new Vector (base);
 623    Vector* kprim = new Vector (haut);
 624
 625    ier = kprim->renormer ();
 626    if (ier!=HOK)
 627        return ier;
 628
 629    jprim->vectoriel (kprim, base);
 630    ier = jprim->renormer ();
 631    if (ier!=HOK)
 632        return ier;
 633
 634    iprim->vectoriel (jprim, kprim);
 635    transfoVertices  (orig,  iprim, jprim, kprim);
 636    return HOK;
 637 }
 638 // ====================================================== transfoVertices
 639 void Elements::transfoVertices (Vertex* orig, Vector* iprim, Vector* jprim,
 640                                 Vector* kprim)
 641 {
 642    double matrice[DIM3][DIM3]={{iprim->getDx(),jprim->getDx(),kprim->getDx()},
 643                                {iprim->getDy(),jprim->getDy(),kprim->getDy()},
 644                                {iprim->getDz(),jprim->getDz(),kprim->getDz()}};
 645
 646    double matkx = orig->getX();
 647    double matky = orig->getY();
 648    double matkz = orig->getZ();
 649
 650    int nbre = tab_vertex.size ();
 651    for (int nro=0 ; nro<nbre ; nro++)
 652        {
 653        if (tab_vertex[nro] != NULL)
 654            tab_vertex[nro]->setMark (NO_USED);
 655        }
 656
 657    for (int nro=0 ; nro<nbre ; nro++)
 658        {
 659        Vertex* node =  tab_vertex[nro];
 660        if (node != NULL  && node->getMark() == NO_USED)
 661           {
 662           double point [DIM3] = {node->getX(), node->getY(), node->getZ()};
 663           double result[DIM3] = {matkx, matky, matkz};
 664
 665           for (int ni=0 ; ni<DIM3; ni++)
 666               for (int nj=0 ; nj<DIM3; nj++)
 667                   result [ni] += matrice[ni][nj] * point[nj];
 668
 669           node->setCoord (result[dir_x], result[dir_y], result[dir_z]);
 670           node->setMark (IS_USED);
 671           }
 672        }
 673 }
 674 // ====================================================== transform
 675 int Elements::transform (Matrix* matrice)
 676 {
 677                                            // -- 1) Raz Marques
 678    for (int nro=0 ; nro<nbr_hexas ; nro++)
 679        {
 680        Hexa* cell = tab_hexa[nro];
 681        if (cell!=NULL && cell->isHere())
 682            cell->razNodes ();
 683        }
 684                                            // -- 2) Move Nodes
 685    for (int nro=0 ; nro<nbr_hexas ; nro++)
 686        {
 687        Hexa* cell = tab_hexa[nro];
 688        if (cell!=NULL && cell->isHere())
 689            cell->moveNodes (matrice);
 690        }
 691
 692    return HOK;
 693 }
 694 // ====================================================== cutHexas
 695 int Elements::cutHexas  (const Edges& t_edges, int nbcuts)
 696 {
 697                                        // 1) marquage des hexas
 698    el_root->markAll (NO_USED);
 699                                        // 2) Memo noeuds
 700    vector <Quad*> q_amont;
 701    vector <Quad*> q_aval;
 702
 703    int nbnodes = t_edges.size();
 704    vector <Vertex*> v_amont (nbnodes);
 705    vector <Vertex*> v_aval  (nbnodes);
 706
 707    for (int nro=0; nro<nbnodes ; nro++)
 708        {
 709        Edge* arete   = t_edges [nro];
 710        v_amont [nro] = arete->getAmont ();
 711        v_aval  [nro] = arete->getAval  ();
 712        int nbcells = arete->getNbrParents ();
 713
 714        for (int nq=0 ; nq<nbcells ; nq++)
 715            {
 716            Quad* quad = arete->getParent (nq);
 717            if (quad->getMark () != IS_USED)
 718               {
 719               quad->setMark (IS_USED);
 720               int nbcubes = quad->getNbrParents ();
 721               for (int nh=0 ; nh<nbcubes ; nh++)
 722                   {
 723                   Hexa* hexa = quad->getParent (nh);
 724                   if (hexa->getMark () != IS_USED)
 725                      {
 726                      hexa->setMark (IS_USED);
 727                      int namont = hexa->getBase (v_amont[nro], arete);
 728                      int naval  = glob->getOpposedQuad (namont);
 729                      q_amont.push_back (hexa->getQuad  (namont));
 730                      q_aval .push_back (hexa->getQuad  (naval ));
 731
 732                      if (db)
 733                         {
 734                         printf (" Quad = ");
 735                         hexa->printName (", ");
 736                         printf (" Faces = (");
 737                         hexa->getQuad (namont)->printName (", ");
 738                         hexa->getQuad (naval )->printName (")\n");
 739                         }
 740                      }
 741                   }
 742               }
 743            }
 744        }
 745                    // ------------------- Dimensionnement
 746    int nbcells   = q_amont.size ();
 747    nbr_vertex    = nbnodes*(nbcuts+2);
 748    int nbpiliers = nbnodes*(nbcuts+1);            // aretes verticales
 749    int nbpoutres = nbcells*(nbcuts+2)*QUAD4;      // aretes horizontales
 750    nbr_edges     = nbpoutres;
 751    nbr_quads     = nbcells*(nbcuts+1)*QUAD4;      // faces Verticales
 752    nbr_hexas     = nbcells*(nbcuts+1);
 753
 754                    // ------------------- Les noeuds et les aretes verticales
 755    tab_quad.resize   (nbr_quads);
 756    tab_edge.resize   (nbr_edges);
 757    tab_hexa.resize   (nbr_hexas);
 758    tab_vertex.resize (nbr_vertex);
 759    vector <Edge*>    tab_pilier (nbpiliers);
 760
 761    int nbinter = nbcuts + 1;
 762    for (int ned=0; ned<nbnodes ; ned++)
 763        {
 764        t_edges [ned]->remove ();
 765        Vertex* ndamont = v_amont [ned];
 766        Vertex* ndaval  = v_aval  [ned];
 767
 768        double dx = (ndaval->getX() - ndamont->getX()) / nbinter;
 769        double dy = (ndaval->getY() - ndamont->getY()) / nbinter;
 770        double dz = (ndaval->getZ() - ndamont->getZ()) / nbinter;
 771
 772        Vertex* nd0 = tab_vertex [ned] = ndamont;
 773        for (int nc=0; nc<nbcuts ; nc++)
 774            {
 775            int nc1 = nc+1;
 776            Vertex* nd1 = el_root->addVertex (ndamont->getX() + nc1*dx,
 777                                              ndamont->getY() + nc1*dy,
 778                                              ndamont->getZ() + nc1*dz);
 779            tab_vertex [nc1*nbnodes + ned] = nd1;
 780            tab_pilier [nc *nbnodes + ned] = newEdge (nd0, nd1);
 781            nd0 =  nd1;
 782            }
 783        tab_vertex [nbinter*nbnodes + ned] = ndaval;
 784        tab_pilier [nbcuts *nbnodes + ned] = newEdge (nd0, ndaval);
 785        ndamont->setMark (ned);
 786        }
 787                    // ------------------- Les aretes horizontales
 788                    // ------------------- Les faces verticales
 789    int sizelig = nbcells*QUAD4;
 790    for (int nro=0; nro<nbcells ; nro++)
 791        {
 792        Quad* sol  = q_amont [nro];
 793        Quad* toit = q_aval  [nro];
 794        // Display (nro);
 795        for (int ns=0; ns<QUAD4 ; ns++)
 796            {
 797            // Display (ns);
 798            Edge* plinthe = sol->getEdge (ns);
 799            int   nmur  = nro*QUAD4 + ns;
 800            int   nmur0 = plinthe->getMark();
 801            if (nmur0 >= 0)
 802               {
 803               for (int nc=0 ; nc<nbinter ; nc++)
 804                   {
 805                   tab_edge [nc*sizelig + nmur] = tab_edge [nc*sizelig + nmur0];
 806                   tab_quad [nc*sizelig + nmur] = tab_quad [nc *sizelig + nmur0];
 807                   printf (" quad_vertical [%02d] = ", nc*sizelig + nmur);
 808                   printf (" quad_vertical [%02d]\n",  nc*sizelig + nmur0);
 809                   printf ("\n");
 810                   }
 811               tab_edge [nbinter*sizelig+nmur] = tab_edge[nbinter*sizelig+nmur0];
 812               }
 813            else
 814               {
 815               plinthe->setMark (nmur);
 816               Vertex* v1  = sol->getVertex (ns);
 817               Vertex* v2  = sol->getVertex ((ns+1) MODULO QUAD4);
 818               int     nd1 = v1->getMark ();
 819               int     nd2 = v2->getMark ();
 820               Edge*   ed0 = tab_edge [nmur] = plinthe;
 821               Edge*   edtoit = toit->getEdge (ns);
 822               for (int nc=0 ; nc<nbinter ; nc++)
 823                   {
 824                   int   nc1 = nc + 1;
 825                   Edge* ed2 = edtoit;
 826                   if (nc<nbcuts)
 827                      {
 828                      v1  = tab_vertex [nc1*nbnodes + nd1];
 829                      v2  = tab_vertex [nc1*nbnodes + nd2];
 830                      ed2 = newEdge (v1, v2);
 831                      }
 832                   tab_edge [nc1*sizelig + nmur] = ed2;
 833                   tab_quad [nc *sizelig + nmur] = newQuad (ed0,
 834                             tab_pilier [nc*nbnodes + nd1], ed2,
 835                             tab_pilier [nc*nbnodes + nd2]);
 836                   ed0 = ed2;
 837                   printf (" quad_vertical [%02d] = ", nc*sizelig + nmur);
 838                   PrintName (tab_quad [nc *sizelig + nmur]);
 839                   printf ("\n");
 840                   }
 841               }
 842            }
 843        }
 844                    // ------------------- Les faces horizontales
 845                    // ------------------- Les hexas
 846    // Rappel : sizelig = nbcells*QUAD4
 847    for (int nro=0; nro<nbcells ; nro++)
 848        {
 849        Quad* qa = q_amont [nro];
 850        for (int nc=0 ; nc<nbinter ; nc++)
 851            {
 852            int  quadv =  nc*nbcells*QUAD4 + nro*QUAD4;
 853            Quad* qb = q_aval[nro];
 854            if (nc<nbcuts)
 855               {
 856               int  edh   = (nc+1)*nbcells*QUAD4 + nro*QUAD4;
 857               qb = newQuad (tab_edge[edh],   tab_edge[edh+1],
 858                             tab_edge[edh+2], tab_edge[edh+3]);
 859               }
 860            Quad* qc = tab_quad [quadv];
 861            Quad* qd = tab_quad [quadv + 2];
 862            Quad* qe = tab_quad [quadv + 1];
 863            Quad* qf = tab_quad [quadv + 3];
 864            tab_hexa [nc*nbcells + nro] = newHexa (qa, qb, qc, qd, qe, qf);
 865            qa = qb;
 866            }
 867        }
 868    return HOK;
 869 }
 870 END_NAMESPACE_HEXA