src/HEXABLOCK/HexCrossElements.cxx

   1
   2 // C++ : Gestion des cylindres croises
   3
   4 // Copyright (C) 2009-2013  CEA/DEN, EDF R&D
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include "HexCrossElements.hxx"
  24
  25 #include "HexDocument.hxx"
  26 #include "HexVector.hxx"
  27 #include "HexVertex.hxx"
  28 #include "HexHexa.hxx"
  29 #include "HexEdge.hxx"
  30
  31 #include "HexGlobale.hxx"
  32 #include "HexCylinder.hxx"
  33 #include "HexOldShape.hxx"
  34 #include "HexNewShape.hxx"
  35
  36 #include <stdlib.h>
  37
  38 static bool   db  = false;
  39
  40 static const double UnSur2pi = DEMI/M_PI;
  41
  42 BEGIN_NAMESPACE_HEXA
  43
  44 // ====================================================== Constructeur
  45 CrossElements::CrossElements (Document* doc, EnumGrid type)
  46              : Elements (doc)
  47 {
  48    cross_cyl1   = NULL;
  49    cross_cyl1   = NULL;
  50    cross_cyl2   = NULL;
  51    cross_center = NULL;
  52    grid_type    = type;
  53    angle_intermed = angle_inter [CylSmall] = angle_inter [CylBig] =  0;
  54    at_right  = at_left = true;
  55    is_filled = false;
  56    grid_geom = NULL;
  57 }
  58 // ====================================================== resize
  59 void CrossElements::resize ()
  60 {
  61    size_hx  = 2;
  62    size_hy  = S_MAXI;
  63    size_hz [CylSmall] = size_h1z;
  64    size_hz [CylBig]   = size_h2z;
  65
  66    size_vx  = size_hx + 1;
  67    size_vy  = size_hy;
  68    size_vz[CylSmall] = size_v1z;
  69    size_vz[CylBig]   = size_v2z;
  70
  71    nbr_vertex1 = size_vx*size_vy* size_v1z;
  72    nbr_quads1  = nbr_vertex1*DIM3;
  73    nbr_edges1  = nbr_quads1;
  74    nbr_hexas1  = size_hx * size_hy * size_h1z;
  75
  76    nbr_vertex  = nbr_vertex1 + size_vx * size_vy * size_v1z;
  77    nbr_quads   = nbr_vertex*DIM3;
  78    nbr_edges   = nbr_quads;
  79    nbr_hexas   = nbr_hexas1 + size_hx * size_hy * size_h2z;
  80
  81    tab_hexa  .resize (nbr_hexas*2);
  82    tab_quad  .resize (nbr_quads*2);
  83    tab_edge  .resize (nbr_edges*2);
  84    tab_vertex.resize (nbr_vertex*2);
  85
  86    for (int nc=0 ; nc< nbr_hexas  ; nc++) tab_hexa   [nc] = NULL;
  87    for (int nc=0 ; nc< nbr_quads  ; nc++) tab_quad   [nc] = NULL;
  88    for (int nc=0 ; nc< nbr_edges  ; nc++) tab_edge   [nc] = NULL;
  89    for (int nc=0 ; nc< nbr_vertex ; nc++) tab_vertex [nc] = NULL;
  90 }
  91 // ====================================================== indHexa
  92 int CrossElements::indHexa (int cyl, int nx, int ny, int nz)
  93 {
  94    if (cyl<0 || cyl>1)
  95        return NOTHING;
  96    if (   nx < 0 || nx >= size_hx || ny < 0 || ny >= size_hy
  97        || nz < 0 || nz >= size_hz[cyl]) return NOTHING;
  98
  99    int nro = cyl*nbr_hexas1 + nx + size_hx*ny + size_hx*size_hy*nz;
 100    return nro;
 101 }
 102 // ====================================================== indQuad
 103 int CrossElements::indQuad (int cyl, int dd, int nx, int ny, int nz)
 104 {
 105    if (cyl<0 || cyl>1 || dd <0 || dd >= DIM3)
 106        return NOTHING;
 107    if (   nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy
 108        || nz < 0 || nz >= size_vz[cyl]) return NOTHING;
 109
 110    int nro = cyl*nbr_quads1 + nx + size_vx*ny + size_vx*size_vy*nz
 111                                  + size_vx*size_vy*size_vz[cyl]*dd;
 112    return nro;
 113 }
 114 // ====================================================== indEdge
 115 int CrossElements::indEdge (int cyl, int dd, int nx, int ny, int nz)
 116 {
 117    return indQuad (cyl, dd, nx, ny, nz);
 118 }
 119 // ====================================================== indVertex
 120 int CrossElements::indVertex (int cyl, int nx, int ny, int nz)
 121 {
 122    if (cyl<0 || cyl>1)
 123       return NOTHING;
 124    if (   nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy
 125        || nz < 0 || nz >= size_vz[cyl])
 126       return NOTHING;
 127
 128    int nro = cyl*nbr_vertex1 + nx + size_vx*ny + size_vx*size_vy*nz;
 129    return nro;
 130 }
 131 // ------------------------------------------------------------------------
 132 // ====================================================== getHexaIJK
 133 Hexa* CrossElements::getHexaIJK (int cyl, int nx, int ny, int nz)
 134 {
 135    int nro = indHexa (cyl, nx, ny, nz);
 136
 137    if (nro >= 0 && tab_hexa[nro]!= NULL &&  tab_hexa[nro]->isHere ())
 138       return tab_hexa [nro];
 139    else
 140       return NULL;
 141 }
 142 // ====================================================== getQuadIJ
 143 Quad* CrossElements::getQuadIJ (int cyl, int nx, int ny, int nz)
 144 {
 145    int nro = indQuad (cyl, dir_z, nx, ny, nz);
 146    if (nro<0)
 147       return NULL;
 148
 149    return tab_quad [nro];
 150 }
 151 // ====================================================== getQuadJK
 152 Quad* CrossElements::getQuadJK (int cyl, int nx, int ny, int nz)
 153 {
 154    int nro = indQuad (cyl, dir_x, nx, ny, nz);
 155    if (nro<0)
 156       return NULL;
 157    return tab_quad [nro];
 158 }
 159 // ====================================================== getQuadIK
 160 Quad* CrossElements::getQuadIK (int cyl, int nx, int ny, int nz)
 161 {
 162    int nro = indQuad (cyl, dir_y, nx, ny, nz);
 163    if (nro<0)
 164       return NULL;
 165    return tab_quad [nro];
 166 }
 167 // ====================================================== getEdgeI
 168 Edge* CrossElements::getEdgeI (int cyl, int nx, int ny, int nz)
 169 {
 170    int nro = indEdge (cyl, dir_x, nx, ny, nz);
 171    if (nro<0)
 172       return NULL;
 173    return tab_edge [nro];
 174 }
 175 // ====================================================== getEdgeJ
 176 Edge* CrossElements::getEdgeJ (int cyl, int nx, int ny, int nz)
 177 {
 178    int nro = indEdge (cyl, dir_y, nx, ny, nz);
 179    if (nro<0)
 180       return NULL;
 181    return tab_edge [nro];
 182 }
 183 // ====================================================== getEdgeK
 184 Edge* CrossElements::getEdgeK (int cyl, int nx, int ny, int nz)
 185 {
 186    int nro = indEdge (cyl, dir_z, nx, ny, nz);
 187    if (nro<0)
 188       return NULL;
 189    return tab_edge [nro];
 190 }
 191 // ====================================================== getVertexIJK
 192 Vertex* CrossElements::getVertexIJK (int cyl, int nx, int ny, int nz)
 193 {
 194    int nro = indVertex (cyl, nx, ny, nz);
 195    if (nro<0)
 196       return NULL;
 197    return tab_vertex [nro];
 198 }
 199 // ------------------------------------------------------------------------
 200 // ====================================================== setHexa
 201 void CrossElements::setHexa (Hexa* elt, int cyl, int nx, int ny, int nz)
 202 {
 203    if (db)
 204       {
 205       printf ("tab_hexa [%d, %d,%d,%d] = ", cyl, nx, ny, nz);
 206       PrintName (elt);
 207       printf ("\n");
 208       }
 209
 210    int nro = indHexa (cyl, nx, ny, nz);
 211    if (nro<0)
 212       return;
 213    tab_hexa [nro] = elt;
 214 }
 215 // ====================================================== setQuad
 216 void CrossElements::setQuad (Quad* elt, int cyl, int dd, int nx, int ny, int nz)
 217 {
 218    if (db)
 219       {
 220       printf ("tab_quad [%d,%d, %d,%d,%d] = ", cyl, dd, nx, ny, nz);
 221       PrintName (elt);
 222       printf ("\n");
 223       }
 224
 225    int nro = indQuad (cyl, dd, nx, ny, nz);
 226    if (nro<0)
 227       return;
 228    tab_quad [nro] = elt;
 229 }
 230 // ====================================================== setEdge
 231 void CrossElements::setEdge (Edge* elt, int cyl, int dd, int nx, int ny, int nz)
 232 {
 233    if (db)
 234       {
 235       printf ("tab_edge [%d,%d, %d,%d,%d] = ", cyl, dd, nx, ny, nz);
 236       PrintName (elt);
 237       printf ("\n");
 238       }
 239
 240    int nro = indEdge (cyl, dd, nx, ny, nz);
 241    if (nro<0)
 242       return;
 243    tab_edge [nro] = elt;
 244 }
 245 // ====================================================== setVertex
 246 void CrossElements::setVertex (Vertex* elt, int cyl, int nx, int ny, int nz)
 247 {
 248    if (db)
 249       {
 250       printf ("tab_vertex [%d, %d,%d,%d] = ", cyl, nx, ny, nz);
 251       PrintName (elt);
 252       printf ("\n");
 253       }
 254
 255    int nro = indVertex (cyl, nx, ny, nz);
 256    if (nro<0)
 257       return;
 258    tab_vertex [nro] = elt;
 259 }
 260 // ====================================================== setVertex (2)
 261 inline bool isequals (double v1, double v2)
 262 {
 263    const double eps    = 0.01;
 264    bool         equals = v1 <= v2 + eps && v1 >= v2 - eps;
 265    return equals;
 266 }
 267 // ====================================================== setVertex (2)
 268 void CrossElements::setVertex (int cyl, int nx, int ny, int nz, double px,
 269                                                      double py, double pz)
 270 {
 271    if (isequals (px, 0) && isequals (py, -1.5)  && isequals (pz, 5))
 272       printf (" Vertex trouve : Cyl%d [%d,%d,%d] = (%g,%g,%g)\n",
 273                                           cyl, nx,ny,nz, px,py,pz);
 274
 275    int nro = indVertex (cyl, nx, ny, nz);
 276    if (nro<0)
 277       return;
 278    else if (tab_vertex[nro] != NULL)
 279       {
 280       Vertex* node = tab_vertex[nro];
 281       if (node->definedBy (px, py, pz))
 282          return;
 283
 284       printf (" ************ ATTENTION ****************\n");
 285       printf (" Creation d'un vertex : Cyl%d [%d,%d,%d] = (%g,%g,%g)\n",
 286                                           cyl, nx,ny,nz, px,py,pz);
 287       printf (" Indice %d deja occupe par le vertex (%g,%g,%g)\n", nro,
 288               node->getX(), node->getY(), node->getZ());
 289
 290       return;
 291       }
 292
 293    int trouve = findVertex (px, py, pz);
 294    if (trouve>=0)
 295       {
 296       printf (" Creation d'un vertex : Cyl%d [%d,%d,%d] = (%g,%g,%g)\n",
 297                                           cyl, nx,ny,nz, px,py,pz);
 298       printf (" Vertex %d present a l'indice %d\n", nro, trouve);
 299       tab_vertex [nro] = tab_vertex [trouve];
 300       return;
 301       }
 302    Vertex*    node = el_root->addVertex (px, py, pz);
 303    setVertex (node, cyl, nx, ny, nz);
 304 }
 305 // ====================================================== copyEdge
 306 void CrossElements::copyEdge (int d1, int i1, int j1, int k1, int d2, int i2,
 307                                       int j2, int k2)
 308 {
 309    Edge* edge = NULL;
 310    switch (d1)
 311           {
 312           case dir_x : edge = getEdgeI (CylSmall, i1, j1, k1);
 313                        break;
 314           case dir_y : edge = getEdgeJ (CylSmall, i1, j1, k1);
 315                        break;
 316           case dir_z : edge = getEdgeK (CylSmall, i1, j1, k1);
 317                        break;
 318           }
 319
 320    setEdge (edge, CylBig, d2, i2, j2, k2);
 321 }
 322 // ====================================================== copyQuad
 323 void CrossElements::copyQuad (int d1, int i1, int j1, int k1, int d2, int i2,
 324                                               int j2, int k2)
 325 {
 326    Quad* quad = NULL;
 327    switch (d1)
 328           {
 329           case dir_x : quad = getQuadJK (CylSmall, i1, j1, k1);
 330                        break;
 331           case dir_y : quad = getQuadIK (CylSmall, i1, j1, k1);
 332                        break;
 333           case dir_z : quad = getQuadIJ (CylSmall, i1, j1, k1);
 334                        break;
 335           }
 336
 337    setQuad (quad, CylBig, d2, i2, j2, k2);
 338 }
 339 // ====================================================== addEdge
 340 Edge* CrossElements::addEdge (Vertex* v1, Vertex* v2, int cyl, int dir,
 341                              int nx, int ny, int nz)
 342 {
 343    Edge* edge = NULL;
 344    if (v1==NULL || v2==NULL)
 345       return NULL;
 346
 347    else if (cyl==CylBig)
 348       {
 349       edge = findEdge1 (v1, v2);
 350       /* ************************************************
 351       if (edge != NULL)
 352          {
 353          printf (" Edge (%d, %d,%d,%d) trouve = ", dir, nx, ny, nz);
 354          edge->printName ("\n");
 355          }
 356          ************************************************ */
 357       }
 358
 359    if (edge == NULL)
 360       edge = newEdge (v1, v2);
 361
 362    setEdge (edge, cyl, dir, nx, ny, nz);
 363    return edge;
 364 }
 365 // ====================================================== addHexa
 366 Hexa* CrossElements::addHexa (Quad* q1, Quad* q2, Quad* q3, Quad* q4, Quad* q5,
 367                             Quad* q6, int cyl, int nx, int ny, int nz)
 368 {
 369 /* **************************
 370    if (cyl==CylSmall)
 371       {
 372       if (nx == 0)  return NULL;
 373       if (nz != 1)  return NULL;
 374       if (ny != 4)  return NULL;
 375       return NULL;
 376       }
 377
 378    if (cyl==CylBig)
 379       {
 380       if (nz >  2)  return NULL;
 381       if (nz <  1)  return NULL;
 382       // if (nx == 0)  return NULL;
 383       // if (ny != 4)  return NULL;
 384       }
 385    ************************** */
 386
 387    Hexa* hexa = NULL;
 388    if (cyl==CylBig)
 389       hexa = findHexa1 (q1, q2);
 390
 391    if (hexa == NULL)
 392        hexa = newHexa (q1, q2, q3, q4, q5, q6);
 393    else if (db)
 394       {
 395       printf (" Hexa (%d,%d,%d) trouve = ", nx, ny, nz);
 396       hexa->printName ("\n");
 397       }
 398
 399    setHexa (hexa, cyl, nx, ny, nz);
 400    return   hexa;
 401 }
 402 // ====================================================== addQuad
 403 Quad* CrossElements::addQuad (Edge* e1, Edge* e2, Edge* e3, Edge* e4, int cyl,
 404                              int dir, int nx, int ny, int nz)
 405 {
 406    Quad* quad = NULL;
 407    if (cyl==CylBig)
 408       quad = findQuad1 (e1, e3);
 409
 410    if (quad == NULL)
 411        quad = newQuad (e1, e2, e3, e4);
 412    else if (db)
 413       {
 414       printf (" Quad (%d, %d,%d,%d) trouve = ", dir, nx, ny, nz);
 415       quad->printName ("\n");
 416       }
 417
 418    setQuad (quad, cyl, dir, nx, ny, nz);
 419    return quad;
 420 }
 421 // ====================================================== findEdge1
 422 Edge* CrossElements::findEdge1 (Vertex* v1, Vertex* v2)
 423 {
 424    for (int nc=0; nc < nbr_edges1 ; nc++)
 425        if (tab_edge[nc] != NULL && tab_edge[nc]->definedBy (v1, v2))
 426           return tab_edge[nc];
 427
 428    return NULL;
 429 }
 430 // ====================================================== findHexa1
 431 Hexa* CrossElements::findHexa1 (Quad* e1, Quad* e2)
 432 {
 433    for (int nc=0; nc < nbr_hexas1 ; nc++)
 434        if (tab_hexa[nc] != NULL && tab_hexa[nc]->definedBy (e1, e2))
 435           return tab_hexa[nc];
 436
 437    return NULL;
 438 }
 439 // ====================================================== findQuad1
 440 Quad* CrossElements::findQuad1 (Edge* e1, Edge* e2)
 441 {
 442    for (int nc=0; nc < nbr_quads1 ; nc++)
 443        if (tab_quad[nc] != NULL && tab_quad[nc]->definedBy (e1, e2))
 444           return tab_quad[nc];
 445
 446    return NULL;
 447 }
 448 // ====================================================== fillGrid
 449 void CrossElements::fillGrid (int cyl, int deb, int fin)
 450 {
 451    int nz0 = deb >= 0 ? deb : 0;
 452    int nzn = fin >  0 ? fin : size_vz [cyl];
 453
 454    fillCenter (cyl, nz0, nzn);
 455
 456    for (int nz=nz0 ; nz<nzn ; nz++)
 457        {
 458        for (int nx=1 ; nx<size_hx ; nx++)
 459            {
 460                                  //   Edges horizontaux radiaux
 461                                  //   Edges horizontaux // cloisons ext
 462            for (int ny=0 ; ny<size_vy ; ny++)
 463                {
 464                int ny1  = (ny + 1) MODULO size_vy;
 465                Vertex* vc = getVertexIJK (cyl, nx,   ny, nz);
 466                Vertex* v1 = getVertexIJK (cyl, nx+1, ny, nz);
 467                Vertex* v2 = getVertexIJK (cyl, nx+1, ny1,nz);
 468
 469                addEdge (vc, v1, cyl, dir_x, nx,   ny, nz);
 470                addEdge (v1, v2, cyl, dir_y, nx+1, ny, nz);
 471                }
 472                                  //   Quads horizontaux
 473            for (int ny=0 ; ny<size_vy ; ny++)
 474                {
 475                int ny1  = (ny + 1) MODULO size_vy;
 476                Edge* e1 = getEdgeI (cyl, nx,   ny,  nz);
 477                Edge* e2 = getEdgeJ (cyl, nx+1, ny,  nz);
 478                Edge* e3 = getEdgeI (cyl, nx,   ny1, nz);
 479                Edge* e4 = getEdgeJ (cyl, nx,   ny,  nz);
 480
 481                addQuad (e1, e2, e3, e4, cyl, dir_z, nx, ny, nz);
 482                }
 483
 484            if (nz>0)
 485               {
 486                                  //   Edges verticaux + cloisons interieures
 487               for (int ny=0 ; ny<size_vy ; ny++)
 488                   {
 489                   Vertex* vp = getVertexIJK (cyl, nx+1, ny, nz-1);
 490                   Vertex* vv = getVertexIJK (cyl, nx+1, ny, nz);
 491
 492                   Edge* edge = addEdge (vp, vv, cyl, dir_z, nx+1, ny, nz-1);
 493
 494                   Edge* e0 = getEdgeI (cyl, nx, ny,  nz-1);
 495                   Edge* e2 = getEdgeI (cyl, nx, ny,  nz);
 496                   Edge* e3 = getEdgeK (cyl, nx, ny,  nz-1);
 497
 498                   addQuad (e0, edge, e2, e3, cyl, dir_y, nx, ny, nz-1);
 499                   }
 500                                  //   Cloisons exterieures ***
 501               for (int ny=0 ; ny<size_vy ; ny++)
 502                   {
 503                   int ny1  = (ny + 1) MODULO size_vy;
 504                   Edge* e0 = getEdgeJ (cyl, nx+1, ny,  nz-1);
 505                   Edge* e2 = getEdgeJ (cyl, nx+1, ny,  nz);
 506
 507                   Edge* e1 = getEdgeK (cyl, nx+1, ny,  nz-1);
 508                   Edge* e3 = getEdgeK (cyl, nx+1, ny1, nz-1);
 509                   addQuad (e0, e1, e2, e3, cyl, dir_x, nx+1, ny, nz-1);
 510                   }
 511                                  //   Hexas (8)
 512               if (is_filled || cyl==CylBig || (nz!=3 && nz != 4))
 513                   {
 514                   for (int ny=0 ; ny<size_hy ; ny++)
 515                       {
 516                       int ny1  = (ny + 1) MODULO size_vy;
 517                       // printf (" ---------- Hexa : nz=%d, ny=%d\n", ny, nz);
 518                       Quad* qa = getQuadIJ (cyl, nx, ny, nz-1);
 519                       Quad* qb = getQuadIJ (cyl, nx, ny, nz);
 520
 521                       Quad* qc = getQuadJK (cyl, nx+1, ny,  nz-1);
 522                       Quad* qd = getQuadJK (cyl, nx,   ny,  nz-1);
 523
 524                       Quad* qe = getQuadIK (cyl, nx,   ny1, nz-1);
 525                       Quad* qf = getQuadIK (cyl, nx,   ny,  nz-1);
 526
 527                       // Hexa* cell = newHexa (qa, qb, qc, qd, qe, qf);
 528                       // setHexa (cell, cyl, nx, ny, nz-1);
 529                       addHexa (qa, qb, qc, qd, qe, qf, cyl, nx, ny, nz-1);
 530                       }
 531                   }
 532               }
 533           }
 534        }
 535 }
 536 // ====================================================== fillCenter
 537 void CrossElements::fillCenter (int cyl, int nz0, int nzn)
 538 {
 539    int nx = 0;
 540
 541    for (int nz=nz0 ; nz<nzn ; nz++)
 542        {
 543        Vertex* center = getVertexIJK (cyl, 0, 0, nz);
 544                                  //   Edges horizontaux // cloisons ext
 545            for (int ny=0 ; ny<size_vy ; ny++)
 546                {
 547                Vertex* v1 = getVertexIJK (cyl, 1, ny, nz);
 548                Vertex* v2 = getVertexIJK (cyl, 1, (ny+1) MODULO size_vy, nz);
 549                addEdge (v1, v2, cyl, dir_y, nx+1, ny, nz);
 550                }
 551        if (is_filled)
 552           {
 553                                  //   Edges horizontaux radiaux
 554            for (int nc=0 ; nc<NbrIntCotes ; nc++)
 555                {
 556                Vertex* vv = getVertexIJK (cyl, 1, 2*nc, nz);
 557                addEdge (center, vv, cyl, dir_x, nx, nc, nz);
 558                }
 559                                  //   Quads horizontaux
 560            for (int nc=0; nc<NbrIntCotes ; nc++)
 561                {
 562                int nc1  = (nc + 1) MODULO NbrIntCotes;
 563                Edge* e1 = getEdgeI (cyl, nx, nc,    nz);
 564                Edge* e2 = getEdgeJ (cyl, nx+1, 2*nc,  nz);
 565                Edge* e3 = getEdgeJ (cyl, nx+1, 2*nc+1, nz);
 566                Edge* e4 = getEdgeI (cyl, nx, nc1,    nz);
 567
 568                addQuad (e1, e2, e3, e4, cyl, dir_z, nx, nc, nz);
 569                }
 570           }
 571
 572        if (nz>0)
 573           {
 574                                  //   Edges verticaux + cloisons interieures
 575           Vertex* vhaut  = getVertexIJK (cyl,nx,0, nz-1);
 576           Edge*   pilier = addEdge (center, vhaut, cyl, dir_z, nx, 0, nz-1);
 577
 578           for (int ny=0 ; ny<size_vy ; ny++)
 579               {
 580               Vertex* vp = getVertexIJK (cyl, 1, ny, nz-1);
 581               Vertex* vv = getVertexIJK (cyl, 1, ny, nz);
 582               Edge* edge = addEdge (vp, vv, cyl, dir_z, 1, ny, nz-1);
 583               if (is_filled && (ny MODULO 2) == 0)
 584                  {
 585                  Edge* e0 = getEdgeI (cyl, nx, ny/2,  nz-1);
 586                  Edge* e2 = getEdgeI (cyl, nx, ny/2,  nz);
 587                  addQuad (e0,edge, e2,pilier, cyl, dir_y, nx, ny/2, nz-1);
 588                  }
 589               }
 590                                  //   Cloisons exterieures
 591           for (int ny=0 ; ny<size_vy ; ny++)
 592               {
 593               int ny1  = (ny + 1) MODULO size_vy;
 594               Edge* e0 = getEdgeJ (cyl, 1, ny,  nz-1);
 595               Edge* e2 = getEdgeJ (cyl, 1, ny,  nz);
 596               Edge* e1 = getEdgeK (cyl, 1, ny,  nz-1);
 597               Edge* e3 = getEdgeK (cyl, 1, ny1, nz-1);
 598
 599               addQuad (e0, e1, e2, e3, cyl, dir_x, 1, ny, nz-1);
 600               }
 601                                  //   Hexas (4)
 602           if (is_filled)
 603              {
 604              for (int nc=0 ; nc < NbrIntCotes ; nc++)
 605                  {
 606                   // printf (" --------------- Hexa : nz=%d, nc=%d\n", nc, nz);
 607                  int nc1  = (nc + 1) MODULO NbrIntCotes;
 608                  Quad* qa = getQuadIJ (cyl, 0, nc, nz-1);
 609                  Quad* qb = getQuadIJ (cyl, 0, nc, nz);
 610
 611                  Quad* qc = getQuadJK (cyl, 1, 2*nc,   nz-1);
 612                  Quad* qe = getQuadJK (cyl, 1, 2*nc+1, nz-1);
 613
 614                  Quad* qd = getQuadIK (cyl, 0, nc1, nz-1);
 615                  Quad* qf = getQuadIK (cyl, 0, nc,  nz-1);
 616
 617                   // Hexa* cell = newHexa (qa, qb, qc, qd, qe, qf);
 618                   // setHexa (cell, cyl, 0, nc, nz-1);
 619                  addHexa (qa, qb, qc, qd, qe, qf, cyl, nx, nc, nz-1);
 620                  }
 621               }
 622           }
 623        }
 624 }
 625 // ====================================================== dump
 626 void CrossElements::dump ()
 627 {
 628    int sizey [2] = { 4, S_MAXI };
 629    for (int cyl=CylSmall ; cyl<=CylBig ; cyl++)
 630        {
 631        printf (" +++++++++ \n");
 632        printf (" +++++++++ Hexas du Cylindre nro %d\n", cyl);
 633        printf (" +++++++++ \n");
 634        for (int nz=0; nz<size_hz[cyl]; nz++)
 635            {
 636            for (int nx=0; nx<size_hx; nx++)
 637                {
 638                printf ("\n");
 639                for (int ny=0; ny<sizey[nx]; ny++)
 640                    {
 641                    Hexa* cell = getHexaIJK (cyl,nx,ny,nz);
 642                    int nro = indHexa (cyl, nx, ny, nz);
 643                    printf ("tab_hexa[%03d] (%d, %d,%d,%d) = ",
 644                            nro, cyl,nx,ny,nz);
 645                    if (cell!=NULL) cell->printName("\n");
 646                       else         printf ("NULL\n");
 647                    }
 648                }
 649            }
 650        }
 651 }
 652 // ====================================================== dumpVertex
 653 void CrossElements::dumpVertex ()
 654 {
 655    int sizey [3] = { 1, S_MAXI, S_MAXI };
 656
 657    for (int cyl=CylSmall ; cyl<=CylBig ; cyl++)
 658        {
 659        printf (" +++++++++ \n");
 660        printf (" +++++++++ Vertex du Cylindre nro %d\n", cyl);
 661        printf (" +++++++++ \n");
 662        for (int nz=0; nz<size_vz[cyl]; nz++)
 663            {
 664            for (int nx=0; nx<size_vx; nx++)
 665                {
 666                printf ("\n");
 667                for (int ny=0; ny<sizey[nx]; ny++)
 668                    {
 669                    Vertex* node = getVertexIJK (cyl,nx,ny,nz);
 670                    int nro = indVertex (cyl, nx, ny, nz);
 671                    printf ("tab_vertex[%03d] (%d, %d,%d,%d) = ",
 672                            nro, cyl,nx,ny,nz);
 673                    if (node!=NULL) node->printName("\n");
 674                       else         printf ("NULL\n");
 675                    }
 676                }
 677            }
 678        }
 679 }
 680 // ====================================================== dumpHexas
 681 void CrossElements::dumpHexas ()
 682 {
 683    int sizey [3] = { 1, 4, S_MAXI };
 684
 685    for (int cyl=CylSmall ; cyl<=CylBig ; cyl++)
 686        {
 687        printf (" +++++++++ \n");
 688        printf (" +++++++++ Hexaedres du Cylindre nro %d\n", cyl);
 689        printf (" +++++++++ \n");
 690        for (int nz=0; nz<size_hz[cyl]; nz++)
 691            {
 692            for (int nx=0; nx<size_hx; nx++)
 693                {
 694                printf ("\n");
 695                for (int ny=0; ny<sizey[nx]; ny++)
 696                    {
 697                    Hexa* elt = getHexaIJK (cyl,nx,ny,nz);
 698                    if (elt!=NULL)
 699                       {
 700                       int nro = indHexa (cyl, nx, ny, nz);
 701                       printf ("tab_hexa[%03d] (%d, %d,%d,%d) = ",
 702                               nro, cyl,nx,ny,nz);
 703                       elt->printName("\n");
 704                       }
 705                    }
 706                }
 707            }
 708        }
 709 }
 710 // ====================================================== calcul_centre
 711 double calcul_centre (Vertex* orig, Vertex* inter)
 712 {
 713    double dx = inter->getX () - orig->getX ();
 714    double dy = inter->getY () - orig->getY ();
 715    double dz = inter->getZ () - orig->getZ ();
 716    double dd = sqrt (dx*dx + dy*dy + dz*dz);
 717    return dd;
 718 }
 719 // ===================================================== assoSlice
 720 void CrossElements::assoSlice (int cyl, double* base, double* normal, int nx,
 721                                                                       int nzs)
 722 {
 723    Real3  center, pnt1, pnt2;
 724    // string brep;
 725
 726    Vertex* v_n  = getVertexIJK (cyl, nx, S_N , nzs);
 727    Vertex* v_s  = getVertexIJK (cyl, nx, S_S , nzs);
 728
 729    v_s->getPoint (pnt1);
 730    v_n->getPoint (pnt2);
 731
 732    double rayon  = calc_distance (pnt1, pnt2)/2;
 733    for (int nro=0 ; nro<DIM3 ; nro++)
 734        center[nro] = (pnt1[nro] + pnt2[nro])/2;
 735
 736    int subid = grid_geom->addCircle (center, rayon, normal, base);
 737
 738    for (int ny=0 ; ny<S_MAXI ; ny++)
 739        {
 740        assoArc (cyl, nx, ny, nzs, subid);
 741        }
 742 }
 743 // ===================================================== assoArc
 744 void CrossElements::assoArc (int cyl, int nx, int ny, int nz, int subid)
 745 {
 746     double angle1 = getAngle (cyl, ny,   nz);
 747     double angle2 = getAngle (cyl, ny+1, nz);
 748     Edge*  edge   = getEdgeJ (cyl, nx, ny, nz);
 749     if (edge==NULL)
 750        return;
 751
 752     grid_geom->addAssociation (edge, subid, angle1*UnSur2pi, angle2*UnSur2pi);
 753
 754     // Assurer avec les vertex
 755     int     ny2   = ny>=S_MAXI-1 ? 0 : ny+1;
 756     Vertex* node1 = getVertexIJK (cyl, nx, ny,  nz);
 757     Vertex* node2 = getVertexIJK (cyl, nx, ny2, nz);
 758
 759     grid_geom->addAssociation (node1, subid, angle1*UnSur2pi);
 760     grid_geom->addAssociation (node2, subid, angle2*UnSur2pi);
 761 }
 762 // ====================================================== getAngleInter
 763 // ==== Angle intersection intermediaire
 764 double CrossElements::getAngleInter (int cyl, int nz)
 765 {
 766    if (cyl==CylBig && (nz==1 || nz==3))
 767       return angle_intermed;
 768    else
 769       return angle_inter[cyl];
 770 }
 771 // ====================================================== getAngle
 772 double CrossElements::getAngle (int cyl, int nj, int nz)
 773 {
 774    switch (nj)
 775       {
 776       case S_E    : return 0;
 777       case S_NE   : return angle_inter[cyl];
 778       case S_N    : return M_PI/2;
 779       case S_NW   : return M_PI - angle_inter[cyl];
 780       case S_W    : return M_PI;
 781       case S_SW   : return M_PI + angle_inter[cyl];
 782       case S_S    : return 3*M_PI/2;
 783       case S_SE   : return 2*M_PI - angle_inter[cyl];
 784       case S_MAXI : return 2*M_PI;
 785       default     : break;
 786       }
 787    return 0;
 788 }
 789 // ====================================================== addVertex
 790 void CrossElements::addVertex (int cyl, int ni, int nj, int nk, double px,
 791                                double rayon)
 792 {
 793    if (rayon<0.0)
 794        rayon = cross_rayon [cyl][ni];
 795    double theta = getAngle (cyl, nj);
 796
 797    if (cyl==CylSmall)
 798        setVertex (cyl, ni, nj, nk,  px, rayon*cos(theta), rayon*sin(theta));
 799    else
 800        setVertex (cyl, ni, nj, nk,  rayon*cos(theta), rayon*sin(theta), px);
 801 }
 802 // ====================================================== addSlice
 803 void CrossElements::addSlice (int cyl, int ni, int nk, double px, double rayon)
 804 {
 805    for (int nj=0 ; nj < S_MAXI ; nj++)
 806        addVertex (cyl, ni, nj, nk, px, rayon);
 807 }
 808 END_NAMESPACE_HEXA