src/HEXABLOCK/HexBiCylinder.cxx

   1
   2 // C++ : Gestion des cylindres croises
   3
   4 // Copyright (C) 2009-2023  CEA, EDF
   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, or (at your option) any later version.
  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 "HexBiCylinder.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 "HexNewShape.hxx"
  33 #include "HexBiCylinderShape.hxx"
  34
  35 ///  static const double UnSur2pi = DEMI/M_PI;
  36
  37 BEGIN_NAMESPACE_HEXA
  38
  39 static bool db = on_debug ();
  40 static const double cos45 = cos (M_PI/4);
  41 #define NameOf(x) ((x)?x->getName():"Null")
  42
  43 // ====================================================== Constructeur
  44 BiCylinder::BiCylinder (Document* doc)
  45           : Elements (doc)
  46 {
  47    cyl_fill     = false;
  48    grid_type    = GR_BIPIPE;
  49    at_right  = at_left = true;
  50
  51    tab_hexa.push_back   (NULL);
  52    tab_quad.push_back   (NULL);
  53    tab_edge.push_back   (NULL);
  54    tab_vertex.push_back (NULL);
  55    nbr_vertex = nbr_edges = nbr_quads = nbr_hexas = 1;
  56
  57    grid_geom = NULL;
  58 }
  59 // ====================================================== getHexaIJK
  60 Hexa* BiCylinder::getHexaIJK (int cyl, int nx, int ny, int nz)
  61 {
  62    int key = getKey (cyl, nx, ny, nz);
  63    int nro = map_hexa [key];
  64    return tab_hexa [nro];
  65 }
  66 // ====================================================== getQuadIJ
  67 Quad* BiCylinder::getQuadIJ (int cyl, int nx, int ny, int nz)
  68 {
  69    int key = getKey (cyl, dir_z, nx, ny, nz);
  70    int nro = map_quad [key];
  71    return tab_quad [nro];
  72 }
  73 // ====================================================== getQuadJK
  74 Quad* BiCylinder::getQuadJK (int cyl, int nx, int ny, int nz)
  75 {
  76    int key = getKey (cyl, dir_x, nx, ny, nz);
  77    int nro = map_quad [key];
  78    return tab_quad [nro];
  79 }
  80 // ====================================================== getQuadIK
  81 Quad* BiCylinder::getQuadIK (int cyl, int nx, int ny, int nz)
  82 {
  83    int key = getKey (cyl, dir_y, nx, ny, nz);
  84    int nro = map_quad [key];
  85    return tab_quad [nro];
  86 }
  87 // ====================================================== getEdgeI
  88 Edge* BiCylinder::getEdgeI (int cyl, int nx, int ny, int nz)
  89 {
  90    int key = getKey (cyl, dir_x, nx, ny, nz);
  91    int nro = map_edge [key];
  92    return tab_edge [nro];
  93 }
  94 // ====================================================== getEdgeJ
  95 Edge* BiCylinder::getEdgeJ (int cyl, int nx, int ny, int nz)
  96 {
  97    int key = getKey (cyl, dir_y, nx, ny, nz);
  98    int nro = map_edge [key];
  99    return tab_edge [nro];
 100 }
 101 // ====================================================== getEdgeK
 102 Edge* BiCylinder::getEdgeK (int cyl, int nx, int ny, int nz)
 103 {
 104    int key = getKey (cyl, dir_z, nx, ny, nz);
 105    int nro = map_edge [key];
 106    return tab_edge [nro];
 107 }
 108 // ====================================================== getVertexIJK
 109 Vertex* BiCylinder::getVertexIJK (int cyl, int nx, int ny, int nz)
 110 {
 111    int key = getKey (cyl, nx, ny, nz);
 112    int nro = map_vertex [key];
 113    if (db)
 114       {
 115       printf ("getVertexIJK (%d, %d,%d,%d) = V%d = ", cyl, nx, ny, nz, nro);
 116       if (tab_vertex[nro] == NULL)
 117          printf ("NULL\n");
 118       else
 119          printf ("%s = (%g, %g, %g)\n", NameOf (tab_vertex[nro]),
 120                  tab_vertex[nro]->getX(), tab_vertex[nro]->getY(),
 121                  tab_vertex[nro]->getZ());
 122       }
 123    return tab_vertex [nro];
 124 }
 125 // ------------------------------------------------------------------------
 126 // ====================================================== addVertex
 127 Vertex* BiCylinder::addVertex (double px, double py, double pz, int cyl,
 128                                int nx, int ny, int nz)
 129 {
 130    Vertex* node = el_root->addVertex (px, py, pz);
 131    int     key  = getKey (cyl, nx, ny, nz);
 132    tab_vertex.push_back (node);
 133    map_vertex [key] = nbr_vertex;
 134
 135    if (db)
 136       {
 137       printf (" tab_vertex [%d, %d,%d,%d] = V%d = ", cyl, nx, ny, nz,
 138                                                      nbr_vertex);
 139       if (node == NULL)
 140          printf ("NULL\n");
 141       else
 142          printf ("%s = (%g, %g, %g)\n", NameOf (node), node->getX(),
 143                                         node->getY(),    node->getZ());
 144       }
 145    nbr_vertex ++;
 146    return  node;
 147 }
 148 // ====================================================== addEdge
 149 Edge* BiCylinder::addEdge (Vertex* v1, Vertex* v2, int cyl, int dir, int nx,
 150                            int ny, int nz)
 151 {
 152    int key = getKey  (cyl, dir, nx, ny, nz);
 153    int nro = map_edge [key];
 154    if (nro>0)
 155       {
 156       if (db)
 157          {
 158          Edge* edge = tab_edge [nro];
 159
 160          printf ("pres_edge [%d,%d, %d,%d,%d] = E%d = ", cyl, dir, nx, ny, nz,
 161                                                         nbr_edges);
 162          printf ("%s  ((%s, %s)) = (%s,%s)\n", NameOf(edge),
 163                  NameOf(edge->getVertex(0)), NameOf(edge->getVertex(1)),
 164                  NameOf(v1), NameOf(v2));
 165          if (edge->definedBy (v1,v2))
 166             return tab_edge [nro];
 167          else
 168             printf (" **** Incoherence !!\n");
 169          }
 170       }
 171
 172    if (v1==NULL || v2==NULL)
 173       return NULL;
 174
 175    Edge* edge = findEdge (v1, v2);
 176    if (edge!=NULL)
 177        return edge;
 178
 179    edge = newEdge (v1, v2);
 180    map_edge [key] = nbr_edges;
 181    tab_edge.push_back (edge);
 182    nbr_edges ++;
 183
 184    if (db)
 185       {
 186       printf (" tab_edge [%d,%d, %d,%d,%d] = E%d = ", cyl, dir, nx, ny, nz,
 187                                                       nbr_edges-1);
 188       if (edge == NULL)
 189          printf ("NULL\n");
 190       else
 191          printf ("%s = (%s, %s)\n", NameOf(edge),
 192                                     NameOf(edge->getVertex(0)),
 193                                     NameOf(edge->getVertex(1)));
 194       }
 195    return edge;
 196 }
 197 // ====================================================== addQuad
 198 Quad* BiCylinder::addQuad (Edge* e1, Edge* e2, Edge* e3, Edge* e4, int cyl,
 199                            int dir, int nx, int ny, int nz)
 200 {
 201    int key = getKey (cyl, dir, nx, ny, nz);
 202    int nro = map_quad [key];
 203    if (nro>0)
 204       {
 205       if (db)
 206          {
 207          Quad* quad = tab_quad [nro];
 208
 209          printf ("pres_quad [%d,%d, %d,%d,%d] = Q%d = ", cyl, dir, nx, ny, nz,
 210                                                         nbr_quads);
 211          printf ("%s  (%s, %s, %s, %s)\n", NameOf(quad),
 212                  NameOf(e1), NameOf(e2), NameOf(e3), NameOf(e4));
 213          if (quad->definedBy (e1,e3))
 214             return tab_quad [nro];
 215          else
 216             {
 217             printf (" **** Incoherence !!\n");
 218             printf ("%s =  (%s, %s, %s, %s)\n", NameOf(quad),
 219                  NameOf(quad->getEdge(0)), NameOf(quad->getEdge(1)),
 220                  NameOf(quad->getEdge(1)), NameOf(quad->getEdge(2)));
 221             printf ("%s =  (%s, %s, %s, %s)\n", NameOf(quad),
 222                  NameOf(quad->getVertex(0)), NameOf(quad->getVertex(1)),
 223                  NameOf(quad->getVertex(1)), NameOf(quad->getVertex(2)));
 224             }
 225          }
 226       }
 227
 228    Quad* quad = Elements::findQuad (e1, e2);
 229    if (quad!=NULL)
 230       return quad;
 231
 232    quad = newQuad (e1, e2, e3, e4);
 233    map_quad [key] = nbr_quads;
 234    tab_quad.push_back (quad);
 235    nbr_quads ++;
 236
 237    if (db)
 238       {
 239       printf (" tab_quad [%d,%d, %d,%d,%d] = Q%d = ", cyl, dir, nx, ny, nz,
 240                                                       nbr_quads-1);
 241       if (quad == NULL)
 242          printf ("NULL\n");
 243       else
 244          printf ("%s = (%s, %s, %s, %s)\n",  NameOf(quad),
 245                  NameOf(quad->getEdge(0)), NameOf(quad->getEdge(1)),
 246                  NameOf(quad->getEdge(2)), NameOf(quad->getEdge(3)));
 247       }
 248    return quad;
 249 }
 250 // ====================================================== addHexa
 251 Hexa* BiCylinder::addHexa (Quad* q1, Quad* q2, Quad* q3, Quad* q4, Quad* q5,
 252                            Quad* q6, int cyl, int nx, int ny, int nz)
 253 {
 254    int key = getKey (cyl, nx, ny, nz);
 255    int nro = map_hexa [key];
 256    if (nro>0)
 257       {
 258       if (tab_hexa [nro]->definedBy (q1,q2))
 259          {
 260          printf (" tab_hexa [%d, %d,%d,%d] = H%d est deja la\n ",
 261                           cyl, nx, ny, nz, nbr_hexas);
 262          return tab_hexa [nro];
 263          }
 264       printf (" Incoherence : tab_hexa [%d, %d,%d,%d] = H%d = %s\n ",
 265                     cyl, nx, ny, nz, nbr_hexas, tab_hexa [nro]->getName ());
 266       }
 267
 268    Hexa* hexa = newHexa (q1, q2, q3, q4, q5, q6);
 269    map_hexa [key] = nbr_hexas;
 270    tab_hexa.push_back (hexa);
 271
 272    if (db)
 273       {
 274       printf (" tab_hexa [%d, %d,%d,%d] = H%d = ", cyl, nx, ny, nz, nbr_hexas);
 275       if (hexa == NULL)
 276          printf ("NULL\n");
 277       else
 278          printf ("%s = (%s, %s, %s, %s, %s, %s)\n",  NameOf(hexa),
 279                  NameOf(hexa->getQuad(0)), NameOf(hexa->getQuad(1)),
 280                  NameOf(hexa->getQuad(2)), NameOf(hexa->getQuad(3)),
 281                  NameOf(hexa->getQuad(4)), NameOf(hexa->getQuad(5)));
 282       }
 283    nbr_hexas ++;
 284    return   hexa;
 285 }
 286 // ------------------------------------------------------------------------
 287 // ====================================================== findVertex
 288 Vertex* BiCylinder::findVertex (double px, double py, double pz,
 289                                int nx, int ny, int nz)
 290 {
 291    for (it_map=map_vertex.begin() ; it_map!=map_vertex.end() ; ++it_map)
 292        {
 293        int     nro = it_map->second;
 294        Vertex* elt = tab_vertex[nro];
 295        if (elt != NULL && elt->definedBy (px, py, pz))
 296           {
 297           int key = getKey (CylBig, nx, ny, nz);
 298           map_vertex [key] = nro;
 299           if (db)
 300              printf ("findVertex [Big,%d,%d,%d] = V%d = '%d' = %s\n",
 301                        nx, ny, nz, nro, it_map->first,
 302                        NameOf(tab_vertex[nro]));
 303           return elt;
 304           }
 305        }
 306    printf ("**** Recherche du vertex (%g, %g, %g)\n", px, py, pz);
 307    fatal_error ("HexBiCylinder : Vertex non trouve");
 308    return NULL;
 309 }
 310 // ====================================================== findEdge
 311 Edge* BiCylinder::findEdge (Vertex* v1, Vertex* v2)
 312 {
 313    int nbedges = tab_edge.size();
 314    for (int nro = 0 ; nro<nbedges ; nro++)
 315        {
 316        Edge* elt = tab_edge[nro];
 317        if (elt != NULL && elt->definedBy (v1, v2))
 318           return elt;
 319        }
 320    return NULL;
 321 }
 322 // ====================================================== findEdge
 323 Edge* BiCylinder::findEdge (Vertex* v1, Vertex* v2, int dir, int nx,
 324                            int ny, int nz)
 325 {
 326    for (it_map=map_edge.begin() ; it_map!=map_edge.end() ; ++it_map)
 327        {
 328        int   nro = it_map->second;
 329        Edge* elt = tab_edge[nro];
 330        if (elt != NULL && elt->definedBy (v1, v2))
 331           {
 332           int key = getKey (CylBig, dir, nx, ny, nz);
 333           map_edge [key] = nro;
 334           if (db)
 335              printf ("findEdge [%d, %d,%d,%d] = E%d = '%d' = %s\n",
 336                        dir, nx, ny, nz, nro, it_map->first,
 337                        NameOf(elt));
 338           return elt;
 339           }
 340        }
 341    fatal_error ("HexBiCylinder : Edge non trouve");
 342    return NULL;
 343 }
 344 // ====================================================== findQuad
 345 Quad* BiCylinder::findQuad (Vertex* v1, Vertex* v2,
 346                            int dir, int nx, int ny, int nz)
 347 {
 348    for (it_map=map_quad.begin() ; it_map!=map_quad.end() ; ++it_map)
 349        {
 350        int   nro = it_map->second;
 351        Quad* elt = tab_quad[nro];
 352        if (elt != NULL && elt->definedBy (v1, v2))
 353           {
 354           int key = getKey (CylBig, dir, nx, ny, nz);
 355           map_quad [key] = nro;
 356           if (db)
 357              printf ("findQuad [%d, %d,%d,%d] = E%d = '%d' = %s\n",
 358                        dir, nx, ny, nz, nro, it_map->first,
 359                        NameOf(elt));
 360           return elt;
 361           }
 362        }
 363    fatal_error ("HexBiCylinder : Quad non trouve");
 364    return NULL;
 365 }
 366 // ====================================================== findQuad
 367 Quad* BiCylinder::findQuad (Edge* e1, Edge* e2, int dir, int nx, int ny, int nz)
 368 {
 369    for (it_map=map_quad.begin() ; it_map!=map_quad.end() ; ++it_map)
 370        {
 371        int   nro = it_map->second;
 372        Quad* elt = tab_quad[nro];
 373        if (elt != NULL && elt->definedBy (e1, e2))
 374           {
 375           int key = getKey (CylBig, dir, nx, ny, nz);
 376           map_quad [key] = nro;
 377           if (db)
 378              printf ("findQuad [%d, %d,%d,%d] = E%d = '%d' = %s\n",
 379                        dir, nx, ny, nz, nro, it_map->first, NameOf(elt));
 380           return elt;
 381           }
 382        }
 383    fatal_error ("HexBiCylinder : Quad non trouve");
 384    return NULL;
 385 }
 386 /*******************************************
 387 // ====================================================== findQuad
 388 Quad* BiCylinder::findQuad (Edge* e1, Edge* e2)
 389 {
 390    for (it_map=map_quad.begin() ; it_map!=map_quad.end() ; ++it_map)
 391        {
 392        int   nro = it_map->second;
 393        Quad* elt = tab_quad[nro];
 394        if (elt != NULL && elt->definedBy (e1, e2))
 395           {
 396           int key = getKey (CylBig, dir, nx, ny, nz);
 397           map_quad [key] = nro;
 398           if (db)
 399              printf ("findQuad [%d, %d,%d,%d] = E%d = '%d' = %s\n",
 400                        dir, nx, ny, nz, nro, it_map->first, NameOf(elt));
 401           return elt;
 402           }
 403        }
 404    fatal_error ("HexBiCylinder : Quad non trouve");
 405    return NULL;
 406 }
 407  ********************************************/
 408 // ====================================================== findHexa
 409 Hexa* BiCylinder::findHexa (Quad* q1, Quad* q2, int nx, int ny, int nz)
 410 {
 411    for (it_map=map_hexa.begin() ; it_map!=map_hexa.end() ; ++it_map)
 412        {
 413        int   nro = it_map->second;
 414        Hexa* elt = tab_hexa[nro];
 415        if (elt != NULL && elt->definedBy (q1, q2))
 416           {
 417           int key = getKey (CylBig, nx, ny, nz);
 418           map_hexa [key] = nro;
 419           if (db)
 420              printf ("findHexa [%d,%d,%d] = H%d = '%d'\n",
 421                        nx, ny, nz, nro, it_map->first);
 422           return elt;
 423           }
 424        }
 425    fatal_error ("HexBiCylinder : Hexa non trouve");
 426    return NULL;
 427 }
 428 // ====================================================== createLittlePipe
 429 void BiCylinder::createLittlePipe ()
 430 {
 431    Real lg  =  cross_hauteur[CylSmall];
 432    Real h1  =  calc_distance (cross_center, cross_orismall);
 433    Real h2  =  cross_rayext[CylBig]*cos45;
 434    Real h3  = (cross_rayext[CylBig] - cross_rayint[CylBig])*cos45;
 435    Real h5  =  lg-h1;
 436
 437 // double t_haut [NbrVslices] = { -h1, -h2, -h2+h3, h2-h3, h2, lg-h1 };
 438
 439    if (at_left)
 440       {
 441       addSlice (CylSmall, 0, 0, -h1, cross_rayint [CylSmall]);
 442       addSlice (CylSmall, 1, 0, -h1, cross_rayext [CylSmall]);
 443       addSlice (CylSmall, 1, 1, -h2, cross_rayext [CylSmall]);
 444       addSlice (CylSmall, 0, 2, h3-h2, cross_rayint [CylSmall]);
 445       }
 446    else if (cyl_fill)
 447       {
 448       addSlice (CylSmall, 0, 1, -h2, cross_rayint [CylSmall]);
 449
 450       addSlice (CylSmall, 1, 1, -h2, cross_rayext [CylSmall]);
 451       addSlice (CylSmall, 0, 2, h3-h2, cross_rayint [CylSmall]);
 452       }
 453
 454
 455    addSlice (CylSmall, 2, 1, -h2, cross_rayext [CylBig]);
 456    addSlice (CylSmall, 1, 2, h3-h2, cross_rayint [CylBig]);
 457    addSlice (CylSmall, 1, 3, h2-h3, cross_rayint [CylBig]);
 458    addSlice (CylSmall, 2, 4,  h2, cross_rayext [CylBig]);
 459
 460    if (at_right)
 461       {
 462       addSlice (CylSmall, 0, 3, h2-h3, cross_rayint [CylSmall]);
 463       addSlice (CylSmall, 1, 4,  h2, cross_rayext [CylSmall]);
 464       addSlice (CylSmall, 0, 5,  h5, cross_rayint [CylSmall]);
 465       addSlice (CylSmall, 1, 5,  h5, cross_rayext [CylSmall]);
 466       }
 467    else if (cyl_fill)
 468       {
 469       addSlice (CylSmall, 0, 3, h2-h3, cross_rayint [CylSmall]);
 470       addSlice (CylSmall, 1, 4,  h2, cross_rayext [CylSmall]);
 471
 472       addSlice (CylSmall, 0, 4,  h2, cross_rayint [CylSmall]);
 473       }
 474                  //--------------------------------- Remplissage
 475                     //    ka kb kc kd
 476    if (at_left)
 477       {
 478       fillSlice (CylSmall, 0,0, 0,2, 1,1, 1,0);
 479       fillSlice (CylSmall, 0,2, 1,2, 2,1, 1,1);
 480       if (cyl_fill)
 481          {
 482          addCube (CylSmall, 0,0,2);
 483          }
 484       }
 485     else if (cyl_fill)
 486       {
 487       fillSlice (CylSmall, 0,1, 0,2, 1,2, 2,1);
 488       addCube (CylSmall, 1,0);
 489       }
 490     else
 491       addCube ( CylSmall, 1, 2);
 492
 493    fillSlice     (CylSmall, 1,2, 1,3, 2,4, 2,1, true);
 494    if (cyl_fill)
 495       {
 496       addCube   (CylSmall, 2,0, 3);
 497       fillSlice (CylSmall, 0,2, 0,3, 1,3, 1,2);
 498       }
 499
 500    if (at_right)
 501       {
 502       fillSlice (CylSmall, 1,3, 0,3, 1,4, 2,4);
 503       fillSlice (CylSmall, 0,3, 0,5, 1,5, 1,4);
 504       if (cyl_fill)
 505           addCube (CylSmall, 3, 0, 5);
 506       }
 507     else if (cyl_fill)
 508       {
 509       fillSlice (CylSmall, 0,3, 0,4, 2,4, 1,3);
 510       addCube   (CylSmall, 3, 0);
 511       }
 512     else
 513       addCube ( CylSmall, 3, 1);
 514 }
 515 // ====================================================== createBigPipe
 516 void BiCylinder::createBigPipe ()
 517 {
 518    const int cyl = CylBig;
 519    Real lg   = cross_hauteur[cyl];
 520    Real rext = cross_rayext [cyl];
 521    Real rint = cross_rayint [cyl];
 522
 523    Real h1   = calc_distance (cross_center, cross_oribig);
 524    Real h2   = rext * cos45;
 525    Real h3   = lg - h1;
 526    Real dh   = (rext - rint)*cos45;
 527
 528    addSlice (CylBig, 0, 0, -h1,    rint);
 529    addSlice (CylBig, 1, 0, -h1,    rext);
 530    addSlice (CylBig, 0, 1, -h2+dh, rint, true);
 531    addSlice (CylBig, 1, 1, -h2,    rext, true);
 532    addSlice (CylBig, 0, 2,  h2-dh, rint, true);
 533    addSlice (CylBig, 1, 2,  h2,    rext, true);
 534    addSlice (CylBig, 0, 3,  h3,    rint);
 535    addSlice (CylBig, 1, 3,  h3,    rext);
 536
 537                   //   A    B    C    D
 538    if (cyl_fill)
 539       addCube (CylBig, 0);
 540    fillSlice (CylBig, 0,0, 0,1, 1,1, 1,0);
 541    if (cyl_fill)
 542       addCube (CylBig, 2);
 543    fillSlice (CylBig, 0,2, 0,3, 1,3, 1,2);
 544 }
 545 // ====================================================== adjustLittleSlice
 546 void BiCylinder::adjustLittleSlice (int ni, int nk)
 547 {
 548    Vertex* node = getVertexIJK (CylSmall, ni, 0, nk);
 549    if (node==NULL)
 550       return;
 551
 552    double grayon2 = cross_rayext[CylBig] * cross_rayext[CylBig];
 553    double prayon  = cross_rayext[CylSmall];
 554    if (ni==0)
 555       {
 556       if (NOT cyl_fill)
 557           grayon2 = cross_rayint[CylBig] * cross_rayint[CylBig];
 558       prayon  = cross_rayint[CylSmall];
 559       }
 560
 561    for (int nj=0;  nj<NbrCotes ; nj++)
 562        {
 563        node  = getVertexIJK (CylSmall, ni, nj, nk);
 564        double angle = getAngle (nj);
 565        double py = prayon * cos (angle);
 566        double pz = prayon * sin (angle);
 567        double px = sqrt (grayon2-py*py);
 568
 569        double qx = node->getX();
 570        if (qx>=0) node->setX ( px);
 571        else       node->setX (-px);
 572        node->setY (py);
 573        node->setZ (pz);
 574        }
 575 }
 576 // ====================================================== addSlice
 577 void BiCylinder::addSlice (int cyl, int ni, int nk, double hauteur,
 578                            double rayon, bool find)
 579 {
 580    for (int nj=0 ; nj<NbrCotes ; nj++)
 581        {
 582        double theta = getAngle (nj);
 583        double px = rayon*cos(theta);
 584        double py = rayon*sin(theta);
 585                   // Find = forcement le gros cylindre
 586        if (find)
 587           findVertex (px, py, hauteur, ni, nj, nk);
 588        else if (cyl==CylBig)
 589           addVertex  (px, py, hauteur, CylBig, ni, nj, nk);
 590        else    // CylSmall
 591           addVertex (hauteur, px, py, CylSmall, ni, nj, nk);
 592        }
 593 }
 594 // ====================================================== addCarre
 595 void BiCylinder::addCarre (int cyl, int nk, double hauteur, double rayon,
 596                            bool find)
 597 {
 598    for (int nj=0 ; nj<NbrCotes ; nj++)
 599        {
 600        double theta = getAngle (nj);
 601        double px = rayon*cos(theta);
 602        double py = rayon*sin(theta);
 603        if (find)
 604           findVertex (px, py, hauteur, 0, nj, nk);
 605        else if (cyl==CylBig)
 606           addVertex  (px, py, hauteur, CylBig, 0, nj, nk);
 607        else    // CylSmall
 608           addVertex (hauteur, px, py, cyl, 0, nj, nk);
 609        }
 610 }
 611 // ====================================================== adjustLittleSquare
 612 void BiCylinder::adjustLittleSquare (int nk)
 613 {
 614    double grayon2 = cross_rayext[CylBig] * cross_rayext[CylBig];
 615
 616    for (int nj=0 ; nj<NbrCotes ; nj++)
 617        {
 618        Vertex* node = getVertexIJK (CylSmall, 0, nj, nk);
 619        double  py   = node->getY();
 620        double  px   = sqrt (grayon2 - py*py);
 621
 622        if (node->getX() > 0) node->setX ( px);
 623        else                  node->setX (-px);
 624        }
 625 }
 626 // ====================================================== addCube
 627 void BiCylinder::addCube (int cyl, int nk0, int nvi0, int k1)
 628 {
 629 /* *****************************************************************
 630        H=bed  +----bd-----+ bdf=G......0
 631              /|          /|
 632            be |   B    bf |
 633            /  |        /  |
 634     E=bce +----bc-----+...|...bcf=F
 635           |  de     D |   df
 636           | E |       | F |             J
 637          ce   | C     cf  |             ^
 638   D=ade...|...+----ad-|---+ adf=C.....3 |   I
 639           |  /        |  /              |  /
 640           | ae    A   | af              | /
 641           |/          |/                |/
 642     A=ace +----ac-----+ acf=B.....2     +----->  K
 643
 644    *****************************************************************  */
 645    enum { nj0, nj1, nj2, nj3, ni0=0 };
 646    const int nk1  = k1>0 ? k1 : nk0+1;
 647    const int nvi1 = nvi0==1  ? 2
 648                   : nvi0==2  ? 1 : 0;
 649
 650                            // La face F
 651    Vertex* vacf = getVertexIJK (cyl, nvi1, nj2, nk1);
 652    Vertex* vadf = getVertexIJK (cyl, nvi1, nj3, nk1);
 653    Vertex* vbdf = getVertexIJK (cyl, nvi1, nj0, nk1);
 654    Vertex* vbcf = getVertexIJK (cyl, nvi1, nj1, nk1);
 655
 656    Edge* eaf = addEdge (vacf, vadf, cyl, dir_y, ni0, nj2, nk1);
 657    Edge* edf = addEdge (vadf, vbdf, cyl, dir_y, ni0, nj3, nk1);
 658    Edge* ebf = addEdge (vbcf, vbdf, cyl, dir_y, ni0, nj0, nk1);
 659    Edge* ecf = addEdge (vacf, vbcf, cyl, dir_y, ni0, nj1, nk1);
 660
 661    Quad* qf  = addQuad (eaf, edf, ebf, ecf, cyl, dir_z, ni0 , nj0, nk1);
 662
 663    Vertex* vace = getVertexIJK (cyl, nvi0, nj2, nk0);
 664    Vertex* vade = getVertexIJK (cyl, nvi0, nj3, nk0);
 665    Vertex* vbde = getVertexIJK (cyl, nvi0, nj0, nk0);
 666    Vertex* vbce = getVertexIJK (cyl, nvi0, nj1, nk0);
 667
 668    Edge* eac = addEdge (vace, vacf, cyl, dir_z, ni0, nj2, nk0);
 669    Edge* ead = addEdge (vade, vadf, cyl, dir_z, ni0, nj3, nk0);
 670    Edge* ebd = addEdge (vbde, vbdf, cyl, dir_z, ni0, nj0, nk0);
 671    Edge* ebc = addEdge (vbce, vbcf, cyl, dir_z, ni0, nj1, nk0);
 672
 673    Edge* eae = addEdge (vace, vade, cyl, dir_y, ni0, nj2, nk0);
 674    Edge* ede = addEdge (vade, vbde, cyl, dir_y, ni0, nj3, nk0);
 675    Edge* ebe = addEdge (vbce, vbde, cyl, dir_y, ni0, nj0, nk0);
 676    Edge* ece = addEdge (vace, vbce, cyl, dir_y, ni0, nj1, nk0);
 677
 678 // Quad* qe  = addQuad (eae, ede, ebe, ece, cyl, dir_z, ni0 , nj0, nk0);
 679    Quad* qe  = addQuad (eae, ece, ebe, ede, cyl, dir_z, ni0 , nj0, nk0);
 680
 681 // Quad* qa  = addQuad (eac, eaf, ead, eae, cyl, dir_y, ni0 , nj2, nk0);
 682    Quad* qa  = addQuad (eac, eae, ead, eaf, cyl, dir_y, ni0 , nj2, nk0);
 683 // Quad* qd  = addQuad (ead, edf, ebd, ede, cyl, dir_x, ni0 , nj3, nk0);
 684    Quad* qd  = addQuad (ead, ede, ebd, edf, cyl, dir_x, ni0 , nj3, nk0);
 685    Quad* qb  = addQuad (ebc, ebf, ebd, ebe, cyl, dir_y, ni0 , nj0, nk0);
 686    Quad* qc  = addQuad (eac, ecf, ebc, ece, cyl, dir_x, ni0 , nj1, nk0);
 687
 688    addHexa (qa, qb, qc, qd, qe, qf, cyl, ni0, nj0, nk0);
 689 }
 690 // ====================================================== fillSlice
 691 /* *****************************************************************
 692        H=bde  +----bd-----+ bdf=G
 693              /|          /|
 694            be |   B    bf |
 695            /  |        /  |
 696     E=bce +----bc-----+...|...bcf=F
 697           |  de     D |   df
 698           | E |       | F |             J
 699          ce   | C     cf  |             ^
 700   D=ade...|...+----ad-|---+ adf=C       |   I
 701           |  /        |  /              |  /
 702           | ae    A   | af              | /
 703           |/          |/                |/
 704     A=ace +----ac-----+ acf=B           +----->  K
 705
 706    ***************************************************************** */
 707 void BiCylinder::fillSlice (int cyl, int nia, int nka, int nib, int nkb,
 708                                      int nic, int nkc, int nid, int nkd,
 709                                      bool meddle)
 710 {
 711    for (int nj0=0; nj0<NbrCotes ; nj0++)
 712        {
 713        if (meddle) nj0++;
 714        int nj1 = nj0+1;
 715        if (nj1>=NbrCotes) nj1=0;
 716        Vertex* vace = getVertexIJK (cyl, nia, nj0, nka);
 717        Vertex* vacf = getVertexIJK (cyl, nib, nj0, nkb);
 718        Vertex* vadf = getVertexIJK (cyl, nic, nj0, nkc);
 719        Vertex* vade = getVertexIJK (cyl, nid, nj0, nkd);
 720
 721        Vertex* vbce = getVertexIJK (cyl, nia, nj1, nka);
 722        Vertex* vbcf = getVertexIJK (cyl, nib, nj1, nkb);
 723        Vertex* vbdf = getVertexIJK (cyl, nic, nj1, nkc);
 724        Vertex* vbde = getVertexIJK (cyl, nid, nj1, nkd);
 725
 726        Edge* eac = addEdge (vace, vacf, cyl, dir_z, nia, nj0, nka);
 727        Edge* ead = addEdge (vadf, vade, cyl, dir_z, nid, nj0, nkd);
 728        Edge* eae = addEdge (vade, vace, cyl, dir_x, nia, nj0, nka);
 729        Edge* eaf = addEdge (vacf, vadf, cyl, dir_x, nib, nj0, nkb);
 730
 731        Edge* ebc = addEdge (vbce, vbcf, cyl, dir_z, nia, nj1, nka);
 732        Edge* ebd = addEdge (vbdf, vbde, cyl, dir_z, nid, nj1, nkd);
 733        Edge* ebe = addEdge (vbde, vbce, cyl, dir_x, nia, nj1, nka);
 734        Edge* ebf = addEdge (vbcf, vbdf, cyl, dir_x, nib, nj1, nkb);
 735
 736        Edge* ece = addEdge (vace, vbce, cyl, dir_y, nia, nj0, nka);
 737        Edge* ecf = addEdge (vacf, vbcf, cyl, dir_y, nib, nj0, nkb);
 738        Edge* edf = addEdge (vadf, vbdf, cyl, dir_y, nic, nj0, nkc);
 739        Edge* ede = addEdge (vade, vbde, cyl, dir_y, nid, nj0, nkd);
 740
 741   //   Quad* qa  = addQuad (eac, eaf, ead, eae, cyl, dir_y, nia+1, nj0, nka);
 742        Quad* qa  = addQuad (eac, eae, ead, eaf, cyl, dir_y, nia+1, nj0, nka);
 743        Quad* qb  = addQuad (ebc, ebf, ebd, ebe, cyl, dir_y, nia+1, nj1, nka);
 744
 745        Quad* qc  = addQuad (eac, ecf, ebc, ece, cyl, dir_x, nia+1, nj0, nka);
 746  //    Quad* qd  = addQuad (ead, edf, ebd, ede, cyl, dir_x, nid+1, nj0, nkd);
 747        Quad* qd  = addQuad (ead, ede, ebd, edf, cyl, dir_x, nid+1, nj0, nkd);
 748
 749   //   Quad* qe  = addQuad (eae, ede, ebe, ece, cyl, dir_z, nia+1, nj0, nka);
 750        Quad* qe  = addQuad (eae, ece, ebe, ede, cyl, dir_z, nia+1, nj0, nka);
 751        Quad* qf  = addQuad (eaf, edf, ebf, ecf, cyl, dir_z, nib+1, nj0, nkb);
 752
 753        addHexa (qa, qb, qc, qd, qe, qf, cyl, nia+1, nj0, nka);
 754        }
 755 }
 756 // ====================================================== assoCylinders
 757 void BiCylinder::assoCylinders ()
 758 {
 759    char     name [12];
 760    sprintf (name, "grid_%02d", el_id);
 761    grid_geom = el_root->addShape (name, SH_INTER);
 762    grid_geom -> openShape();
 763
 764    int s_kmax  = 5;
 765    int imin    = cyl_fill ? 1 : 0;
 766    int g_ksize = 4;
 767
 768    for (int ni=imin ; ni<2 ; ni++)
 769        {
 770        assoSlice (CylSmall, ni, 0, cross_dirsmall);
 771        assoSlice (CylSmall, ni, s_kmax, cross_dirsmall);
 772
 773        for (int nk=0 ; nk<g_ksize ; nk++)
 774             assoSlice (CylBig, ni, nk, cross_dirbig);
 775       }
 776
 777
 778    if (at_left)
 779       {
 780       assoIntersection (1, 1);
 781       if (NOT cyl_fill)
 782           assoIntersection (0, 2);
 783       }
 784
 785    if (at_right)
 786       {
 787       assoIntersection (1, 4);
 788       if (NOT cyl_fill)
 789           assoIntersection (0, 3);
 790       }
 791
 792    grid_geom -> closeShape();
 793 }
 794 // ====================================================== assoSlice
 795 void BiCylinder::assoSlice (int cyl, int nx, int nzs, double* normal)
 796 {
 797    Real3  center, pnt0, pnt1, pnt2, vbase;
 798    int ny0 = 0;
 799    int ny1 = 1;
 800    int ny2 = NbrCotes/2;
 801
 802    Vertex* v0 = getVertexIJK (cyl, nx, ny0 , nzs);
 803    Vertex* v1 = getVertexIJK (cyl, nx, ny1 , nzs);
 804    Vertex* v2 = getVertexIJK (cyl, nx, ny2 , nzs);
 805
 806    if (v0==NULL || v1==NULL || v2==NULL)
 807       return;
 808
 809    v0->getPoint (pnt0);
 810    v1->getPoint (pnt1);
 811    v2->getPoint (pnt2);
 812
 813    double rayon = 0;
 814    for (int nro=0 ; nro<DIM3 ; nro++)
 815        {
 816        center[nro] = (pnt0[nro] + pnt2[nro])/2;
 817        vbase [nro] = (pnt0[nro] - pnt1[nro]);
 818        rayon += carre (center [nro] - pnt0 [nro]);
 819        }
 820
 821    int subid = grid_geom->addCircle (center, sqrt(rayon), normal, vbase);
 822    for (int ny=0 ; ny<NbrCotes ; ny++)
 823        {
 824        int ny1 = (ny+1) MODULO NbrCotes;
 825        Vertex* va = getVertexIJK (cyl, nx, ny , nzs);
 826        Vertex* vb = getVertexIJK (cyl, nx, ny1, nzs);
 827
 828        assoArc (cyl, ny, va, vb, subid);
 829        }
 830 }
 831 // ===================================================== assoArc
 832 void BiCylinder::assoArc (int cyl, int ny, Vertex* v1, Vertex* v2, int subid)
 833 {
 834     const double Decal = 1.0 / NbrCotes;
 835     const double Start = Decal / 2;
 836
 837     int    ny2    = ny+1;
 838     double posit  = Start + ny*Decal;
 839     Edge*  edge   = findEdge (v1, v2);
 840     if (edge==NULL)
 841        return;
 842                                         // Vertex
 843     grid_geom->addAssociation (v1, subid, posit);
 844
 845     if (ny2 < NbrCotes)
 846        {
 847        grid_geom->addAssociation (edge, subid, posit, posit+Decal);
 848        }
 849     else
 850        {
 851        grid_geom->addAssociation (edge, subid, posit, 1.0);
 852        grid_geom->addAssociation (edge, subid, 0,   Start);
 853        }
 854 }
 855 // ===================================================== assoArc
 856 void BiCylinder::assoArc (int cyl, int nx, int ny, int nz, int subid)
 857 {
 858     const double Decal = 1.0 / NbrCotes;
 859     const double Start = Decal / 2;
 860
 861     int    ny2    = ny+1;
 862     double posit  = Start + ny*Decal;
 863     Edge*  edge   = getEdgeJ (cyl, nx, ny, nz);
 864     if (edge==NULL)
 865        return;
 866
 867                                         // Vertex
 868     Vertex* node = getVertexIJK (cyl, nx, ny,  nz);
 869     grid_geom->addAssociation (node, subid, posit);
 870
 871     if (ny2 < NbrCotes)
 872        {
 873        grid_geom->addAssociation (edge, subid, posit, posit+Decal);
 874        // node = getVertexIJK (cyl, nx, ny2, nz);
 875        // grid_geom->addAssociation (node, subid, angle2*UnSur2pi);
 876        }
 877     else
 878        {
 879        grid_geom->addAssociation (edge, subid, posit, 1.0);
 880        grid_geom->addAssociation (edge, subid, 0,   Start);
 881        }
 882 }
 883 // ===================================================== assoIntersection
 884 int BiCylinder::assoIntersection (int nxs, int nzs)
 885 {
 886    int ier = HOK;
 887    Real3  pse, psw, sorig, sbase;
 888    Real3  pbe, pbw, borig, bbase;
 889    std::string brep;
 890    int ny0 = 0;
 891    int nyd = NbrCotes/2;
 892    int MiddleSlice1 = 3;
 893
 894    int nz = nzs < MiddleSlice1 ? 0 : 5;
 895
 896    Vertex* vse = getVertexIJK (CylSmall, nxs, ny0 , nz);
 897    Vertex* vsw = getVertexIJK (CylSmall, nxs, nyd , nz);
 898    Vertex* vbe = getVertexIJK (CylBig,   nxs, ny0 , 0);
 899    Vertex* vbw = getVertexIJK (CylBig,   nxs, nyd , 0);
 900
 901    if (vse==NULL || vsw==NULL || vbe==NULL || vbw==NULL)
 902       return HERR;
 903
 904    vse->getPoint (pse);
 905    vsw->getPoint (psw);
 906    vbe->getPoint (pbe);
 907    vbw->getPoint (pbw);
 908
 909    double srayon = calc_distance (pse, psw)/2;
 910    double brayon = calc_distance (pbe, pbw)/2;
 911
 912    calc_milieu  (psw, pse, sorig);
 913    calc_milieu  (pbw, pbe, borig);
 914    calc_vecteur (psw, pse, sbase);
 915    calc_vecteur (pbw, pbe, bbase);
 916
 917    double shaut = calc_distance (cross_center, sorig);
 918    double bhaut = calc_distance (cross_center, borig)*2;
 919    double* orig = nzs < MiddleSlice1 ? sorig : cross_center; // Pb orientation
 920
 921    BiCylinderShape bicyl_shape (el_root);
 922    ier = bicyl_shape.defineCyls (borig, cross_dirbig,   bbase, brayon, bhaut,
 923                                   orig, cross_dirsmall, sbase, srayon, shaut);
 924    if (ier != HOK)
 925       return ier;
 926
 927    for (int ny=0 ; ny<NbrCotes ; ny++)
 928        {
 929        Vertex* node = getVertexIJK (CylSmall, nxs, ny, nzs);
 930        if (node!=NULL)
 931            node->clearAssociation ();
 932        // Edge*   edge = getEdgeJ     (CylSmall, nxs, ny, nzs);
 933        // if (edge!=NULL) edge->clearAssociation ();
 934        }
 935
 936    for (int ny=0 ; ny<NbrCotes ; ny++)
 937        {
 938        int ny1 = (ny+1) MODULO QUAD4;
 939        Vertex* node0 = getVertexIJK (CylSmall, nxs, ny, nzs);
 940        Vertex* node1 = getVertexIJK (CylSmall, nxs, ny1, nzs);
 941  //    Edge*   edge  = getEdgeJ (CylSmall, nxs, ny, nzs);
 942        Edge*   edge  = findEdge (node0, node1);
 943        bicyl_shape.associate (edge);
 944        }
 945
 946    return HOK;
 947 }
 948 // ====================================================== makeCylinders
 949 int BiCylinder::makeCylinders (Vertex* ori1, double* vz1, double r1, double h1,
 950                                Vertex* ori2, double* vz2, double r2, double h2)
 951 {
 952    grid_type = GR_BICYL;
 953    int ier = makePipes (ori1,vz1, r1/2, r1,h1, ori2,vz2, r2/2, r2,h2, true);
 954    return ier;
 955 }
 956 // ====================================================== makePipes
 957 int BiCylinder::makePipes (Vertex* ori1, double* vz1, double rint1,
 958                            double rext1, double h1, Vertex* ori2, double* vz2,
 959                            double rint2, double rext2, double h2, bool fill)
 960 {
 961    cyl_fill = fill;
 962    cross_hauteur [CylBig]   = h2;
 963    cross_rayext  [CylBig  ] = rext2;
 964    cross_rayint  [CylBig  ] = rint2;
 965
 966    cross_hauteur [CylSmall] = h1;
 967    cross_rayext  [CylSmall] = rext1;
 968    cross_rayint  [CylSmall] = rint1;
 969
 970    ori1->getPoint (cross_orismall);
 971    ori2->getPoint (cross_oribig  );
 972    copy_vecteur (vz2, cross_dirbig);
 973    copy_vecteur (vz1, cross_dirsmall);
 974
 975 //   Intersection des 2 cylindres :
 976    int ier = checkInter (cross_oribig,   vz2, rext2, h2,
 977                          cross_orismall, vz1, rext1, h1,
 978                          cross_center, at_left, at_right);
 979    if (ier != HOK)
 980       return ier;
 981
 982    createLittlePipe ();
 983    createBigPipe    ();
 984
 985    if (at_left)
 986       {
 987       adjustLittleSlice (1, 1);
 988       adjustLittleSlice (0, 2);
 989       }
 990    else if (cyl_fill)
 991       {
 992       adjustLittleSlice (0, 1);
 993       }
 994
 995    if (at_right)
 996       {
 997       adjustLittleSlice (0, 3);
 998       adjustLittleSlice (1, 4);
 999       }
1000    else if (cyl_fill)
1001       {
1002       adjustLittleSlice (0, 4);
1003       }
1004
1005
1006    transfoVertices (cross_center, cross_dirsmall, cross_dirbig);
1007    assoCylinders   ();
1008
1009    //  assoResiduelle ();
1010    // el_root->reorderQuads ();
1011    return HOK;
1012 }
1013 END_NAMESPACE_HEXA