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