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