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/ 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 "HexCylinder.hxx"
  33
  34 static bool   db  = false;
  35
  36 BEGIN_NAMESPACE_HEXA
  37
  38 static const double cos45 = cos (M_PI/4);
  39
  40
  41 void geom_define_line (string& brep);
  42 void geom_asso_point  (double angle, Vertex* node);
  43
  44 void geom_create_circle (double* milieu, double rayon, double* normale,
  45                          double* base, string& brep);
  46 int  geom_create_cylcyl (double* borig, double* bnorm, double* bbase,
  47                          double  bray,  double  bhaut,
  48                          double* sorig, double* snorm, double* sbase,
  49                          double  sray,  double  shaut);
  50 int  geom_asso_cylcyl   (Edge* edge);
  51
  52 // ====================================================== Constructeur
  53 BiCylinder::BiCylinder (Document* doc)
  54           : Elements (doc)
  55 {
  56    cyl_fill     = false;
  57    cross_cyl1   = NULL;
  58    cross_cyl1   = NULL;
  59    cross_cyl2   = NULL;
  60    grid_type    = GR_BICYL;
  61    at_right  = at_left = true;
  62
  63    tab_hexa.push_back   (NULL);
  64    tab_quad.push_back   (NULL);
  65    tab_edge.push_back   (NULL);
  66    tab_vertex.push_back (NULL);
  67    nbr_vertex = nbr_edges = nbr_quads = nbr_hexas = 1;
  68 }
  69 // ====================================================== getHexaIJK
  70 Hexa* BiCylinder::getHexaIJK (int cyl, int nx, int ny, int nz)
  71 {
  72    int key = getKey (cyl, nx, ny, nz);
  73    int nro = map_hexa [key];
  74    return tab_hexa [nro];
  75 }
  76 // ====================================================== getQuadIJ
  77 Quad* BiCylinder::getQuadIJ (int cyl, int nx, int ny, int nz)
  78 {
  79    int key = getKey (cyl, dir_z, nx, ny, nz);
  80    int nro = map_quad [key];
  81    return tab_quad [nro];
  82 }
  83 // ====================================================== getQuadJK
  84 Quad* BiCylinder::getQuadJK (int cyl, int nx, int ny, int nz)
  85 {
  86    int key = getKey (cyl, dir_x, nx, ny, nz);
  87    int nro = map_quad [key];
  88    return tab_quad [nro];
  89 }
  90 // ====================================================== getQuadIK
  91 Quad* BiCylinder::getQuadIK (int cyl, int nx, int ny, int nz)
  92 {
  93    int key = getKey (cyl, dir_y, nx, ny, nz);
  94    int nro = map_quad [key];
  95    return tab_quad [nro];
  96 }
  97 // ====================================================== getEdgeI
  98 Edge* BiCylinder::getEdgeI (int cyl, int nx, int ny, int nz)
  99 {
 100    int key = getKey (cyl, dir_x, nx, ny, nz);
 101    int nro = map_edge [key];
 102    return tab_edge [nro];
 103 }
 104 // ====================================================== getEdgeJ
 105 Edge* BiCylinder::getEdgeJ (int cyl, int nx, int ny, int nz)
 106 {
 107    int key = getKey (cyl, dir_y, nx, ny, nz);
 108    int nro = map_edge [key];
 109    return tab_edge [nro];
 110 }
 111 // ====================================================== getEdgeK
 112 Edge* BiCylinder::getEdgeK (int cyl, int nx, int ny, int nz)
 113 {
 114    int key = getKey (cyl, dir_z, nx, ny, nz);
 115    int nro = map_edge [key];
 116    return tab_edge [nro];
 117 }
 118 // ====================================================== getVertexIJK
 119 Vertex* BiCylinder::getVertexIJK (int cyl, int nx, int ny, int nz)
 120 {
 121    int key = getKey (cyl, nx, ny, nz);
 122    int nro = map_vertex [key];
 123    if (db)
 124       {
 125       printf ("getVertexIJK (%d, %d,%d,%d) = V%d = ", cyl, nx, ny, nz, nro);
 126       if (tab_vertex[nro] == NULL)
 127          printf ("NULL\n");
 128       else
 129          printf ("%s = (%g, %g, %g)\n", tab_vertex[nro]->getName(),
 130                  tab_vertex[nro]->getX(), tab_vertex[nro]->getY(),
 131                  tab_vertex[nro]->getZ());
 132       }
 133    return tab_vertex [nro];
 134 }
 135 // ------------------------------------------------------------------------
 136 // ====================================================== addVertex
 137 Vertex* BiCylinder::addVertex (double px, double py, double pz, int cyl,
 138                                int nx, int ny, int nz)
 139 {
 140    Vertex* node = el_root->addVertex (px, py, pz);
 141    int     key  = getKey (cyl, nx, ny, nz);
 142    tab_vertex.push_back (node);
 143    map_vertex [key] = nbr_vertex;
 144
 145    if (db)
 146       {
 147       printf (" tab_vertex [%d, %d,%d,%d] = V%d = ", cyl, nx, ny, nz,
 148                                                      nbr_vertex);
 149       if (node == NULL)
 150          printf ("NULL\n");
 151       else
 152          printf ("%s = (%g, %g, %g)\n", node->getName(), node->getX(),
 153                                         node->getY(),    node->getZ());
 154       }
 155    nbr_vertex ++;
 156    return  node;
 157 }
 158 // ====================================================== addEdge
 159 Edge* BiCylinder::addEdge (Vertex* v1, Vertex* v2, int cyl, int dir, int nx,
 160                            int ny, int nz)
 161 {
 162    int key = getKey  (cyl, dir, nx, ny, nz);
 163    int nro = map_edge [key];
 164    if (nro>0)
 165       {
 166       if (db)
 167          {
 168          Edge* edge = tab_edge [nro];
 169
 170          printf ("pres_edge [%d,%d, %d,%d,%d] = E%d = ", cyl, dir, nx, ny, nz,
 171                                                         nbr_edges);
 172          printf ("%s  ((%s, %s)) = (%s,%s)\n", edge->getName(),
 173                  edge->getVertex(0)->getName(), edge->getVertex(1)->getName(),
 174                  v1->getName(), v2->getName());
 175          if (NOT edge->definedBy (v1,v2))
 176             printf (" **** Incoherence !!\n");
 177          }
 178       return tab_edge [nro];
 179       }
 180
 181    if (v1==NULL || v2==NULL)
 182       return NULL;
 183
 184    Edge* edge = findEdge (v1, v2);
 185    if (edge==NULL)
 186        edge = newEdge (v1, v2);
 187
 188    map_edge [key] = nbr_edges;
 189    tab_edge.push_back (edge);
 190
 191    if (db)
 192       {
 193       printf (" tab_edge [%d,%d, %d,%d,%d] = E%d = ", cyl, dir, nx, ny, nz,
 194                                                      nbr_edges);
 195       if (edge == NULL)
 196          printf ("NULL\n");
 197       else
 198          printf ("%s = (%s, %s)\n", edge->getName(),
 199                                     edge->getVertex(0)->getName(),
 200                                     edge->getVertex(1)->getName());
 201       }
 202    nbr_edges ++;
 203    return edge;
 204 }
 205 // ====================================================== addQuad
 206 Quad* BiCylinder::addQuad (Edge* e1, Edge* e2, Edge* e3, Edge* e4, int cyl,
 207                            int dir, int nx, int ny, int nz)
 208 {
 209    int key = getKey (cyl, dir, nx, ny, nz);
 210    int nro = map_quad [key];
 211    if (nro>0)
 212       return tab_quad [nro];
 213
 214    Quad* quad = newQuad (e1, e2, e3, e4);
 215    map_quad [key] = nbr_quads;
 216    tab_quad.push_back (quad);
 217
 218    if (db)
 219       {
 220       printf (" tab_quad [%d,%d, %d,%d,%d] = Q%d = ", cyl, dir, nx, ny, nz,
 221                                                      nbr_quads);
 222       if (quad == NULL)
 223          printf ("NULL\n");
 224       else
 225          printf ("%s = (%s, %s, %s, %s)\n",  quad->getName(),
 226                  quad->getEdge(0)->getName(), quad->getEdge(1)->getName(),
 227                  quad->getEdge(2)->getName(), quad->getEdge(3)->getName());
 228       }
 229    nbr_quads ++;
 230    return quad;
 231 }
 232 // ====================================================== addHexa
 233 Hexa* BiCylinder::addHexa (Quad* q1, Quad* q2, Quad* q3, Quad* q4, Quad* q5,
 234                            Quad* q6, int cyl, int nx, int ny, int nz)
 235 {
 236    int key = getKey (cyl, nx, ny, nz);
 237    int nro = map_hexa [key];
 238    if (nro>0)
 239       {
 240       printf (" tab_hexa [%d, %d,%d,%d] = H%d est deja la\n ",
 241                           cyl, nx, ny, nz, nbr_hexas);
 242       return tab_hexa [nro];
 243       }
 244
 245    Hexa* hexa = newHexa (q1, q2, q3, q4, q5, q6);
 246    map_hexa [key] = nbr_hexas;
 247    tab_hexa.push_back (hexa);
 248
 249    if (db)
 250       {
 251       printf (" tab_hexa [%d, %d,%d,%d] = H%d = ", cyl, nx, ny, nz, nbr_hexas);
 252       if (hexa == NULL)
 253          printf ("NULL\n");
 254       else
 255          printf ("%s = (%s, %s, %s, %s, %s, %s)\n",  hexa->getName(),
 256                  hexa->getQuad(0)->getName(), hexa->getQuad(1)->getName(),
 257                  hexa->getQuad(2)->getName(), hexa->getQuad(3)->getName(),
 258                  hexa->getQuad(4)->getName(), hexa->getQuad(5)->getName());
 259       }
 260    nbr_hexas ++;
 261    return   hexa;
 262 }
 263 // ------------------------------------------------------------------------
 264 // ====================================================== findVertex
 265 Vertex* BiCylinder::findVertex (double px, double py, double pz,
 266                                int nx, int ny, int nz)
 267 {
 268    for (it_map=map_vertex.begin() ; it_map!=map_vertex.end() ; ++it_map)
 269        {
 270        int     nro = it_map->second;
 271        Vertex* elt = tab_vertex[nro];
 272        if (elt != NULL && elt->definedBy (px, py, pz))
 273           {
 274           int key = getKey (CylBig, nx, ny, nz);
 275           map_vertex [key] = nro;
 276           if (db)
 277              printf ("findVertex [Big,%d,%d,%d] = V%d = '%d' = %s\n",
 278                        nx, ny, nz, nro, it_map->first,
 279                        tab_vertex[nro]->getName());
 280           return elt;
 281           }
 282        }
 283    printf ("**** Recherche du vertex (%g, %g, %g)\n", px, py, pz);
 284    fatal_error ("HexBiCylinder : Vertex non trouve");
 285    return NULL;
 286 }
 287 // ====================================================== findEdge
 288 Edge* BiCylinder::findEdge (Vertex* v1, Vertex* v2)
 289 {
 290    int nbedges = tab_edge.size();
 291    for (int nro = 0 ; nro<nbedges ; nro++)
 292        {
 293        Edge* elt = tab_edge[nro];
 294        if (elt != NULL && elt->definedBy (v1, v2))
 295           return elt;
 296        }
 297    return NULL;
 298 }
 299 // ====================================================== findEdge
 300 Edge* BiCylinder::findEdge (Vertex* v1, Vertex* v2, int dir, int nx,
 301                            int ny, int nz)
 302 {
 303    for (it_map=map_edge.begin() ; it_map!=map_edge.end() ; ++it_map)
 304        {
 305        int   nro = it_map->second;
 306        Edge* elt = tab_edge[nro];
 307        if (elt != NULL && elt->definedBy (v1, v2))
 308           {
 309           int key = getKey (CylBig, dir, nx, ny, nz);
 310           map_edge [key] = nro;
 311           if (db)
 312              printf ("findEdge [%d, %d,%d,%d] = E%d = '%d' = %s\n",
 313                        dir, nx, ny, nz, nro, it_map->first,
 314                        elt->getName());
 315           return elt;
 316           }
 317        }
 318    fatal_error ("HexBiCylinder : Edge non trouve");
 319    return NULL;
 320 }
 321 // ====================================================== findQuad
 322 Quad* BiCylinder::findQuad (Vertex* v1, Vertex* v2,
 323                            int dir, int nx, int ny, int nz)
 324 {
 325    for (it_map=map_quad.begin() ; it_map!=map_quad.end() ; ++it_map)
 326        {
 327        int   nro = it_map->second;
 328        Quad* elt = tab_quad[nro];
 329        if (elt != NULL && elt->definedBy (v1, v2))
 330           {
 331           int key = getKey (CylBig, dir, nx, ny, nz);
 332           map_quad [key] = nro;
 333           if (db)
 334              printf ("findQuad [%d, %d,%d,%d] = E%d = '%d' = %s\n",
 335                        dir, nx, ny, nz, nro, it_map->first,
 336                        elt->getName());
 337           return elt;
 338           }
 339        }
 340    fatal_error ("HexBiCylinder : Quad non trouve");
 341    return NULL;
 342 }
 343 // ====================================================== findQuad
 344 Quad* BiCylinder::findQuad (Edge* e1, Edge* e2, int dir, int nx, int ny, int nz)
 345 {
 346    for (it_map=map_quad.begin() ; it_map!=map_quad.end() ; ++it_map)
 347        {
 348        int   nro = it_map->second;
 349        Quad* elt = tab_quad[nro];
 350        if (elt != NULL && elt->definedBy (e1, e2))
 351           {
 352           int key = getKey (CylBig, dir, nx, ny, nz);
 353           map_quad [key] = nro;
 354           if (db)
 355              printf ("findQuad [%d, %d,%d,%d] = E%d = '%d' = %s\n",
 356                        dir, nx, ny, nz, nro, it_map->first, elt->getName());
 357           return elt;
 358           }
 359        }
 360    fatal_error ("HexBiCylinder : Quad non trouve");
 361    return NULL;
 362 }
 363 // ====================================================== findHexa
 364 Hexa* BiCylinder::findHexa (Quad* q1, Quad* q2, int nx, int ny, int nz)
 365 {
 366    for (it_map=map_hexa.begin() ; it_map!=map_hexa.end() ; ++it_map)
 367        {
 368        int   nro = it_map->second;
 369        Hexa* elt = tab_hexa[nro];
 370        if (elt != NULL && elt->definedBy (q1, q2))
 371           {
 372           int key = getKey (CylBig, nx, ny, nz);
 373           map_hexa [key] = nro;
 374           if (db)
 375              printf ("findHexa [%d,%d,%d] = H%d = '%d'\n",
 376                        nx, ny, nz, nro, it_map->first);
 377           return elt;
 378           }
 379        }
 380    fatal_error ("HexBiCylinder : Hexa non trouve");
 381    return NULL;
 382 }
 383 // ====================================================== crossCylinders
 384 int BiCylinder::crossCylinders (Cylinder* lun, Cylinder* lautre)
 385 {
 386    if (lun->getRadius() < lautre->getRadius())
 387       {
 388       cross_cyl1 = lun;
 389       cross_cyl2 = lautre;
 390       }
 391    else
 392       {
 393       cross_cyl1 = lautre;
 394       cross_cyl2 = lun;
 395       }
 396
 397    int ier = cross_cyl2->interCylinder (cross_cyl1, at_left, at_right,
 398                                         cross_center);
 399    if (ier != HOK)
 400       return ier;
 401
 402    cyl_fill = false;
 403    cross_rayext  [CylSmall] = cross_cyl1->getRadius ();
 404    cross_rayext  [CylBig]   = cross_cyl2->getRadius ();
 405    cross_rayint  [CylSmall] = cross_rayext  [CylSmall] / 2;
 406    cross_rayint  [CylBig  ] = cross_rayext  [CylBig  ] / 2;
 407    cross_hauteur [CylSmall] = cross_cyl1->getHeight ();
 408    cross_hauteur [CylBig]   = cross_cyl2->getHeight ();
 409
 410    createLittleCyl ();
 411    createBigCyl    ();
 412    adjustLittleSlice (1, 1);
 413    adjustLittleSlice (0, 2);
 414    adjustLittleSlice (0, 3);
 415    adjustLittleSlice (1, 4);
 416
 417    Vector* iprim = new Vector (cross_cyl1->getDirection());
 418    Vector* kprim = new Vector (cross_cyl2->getDirection());
 419    Vector* jprim = new Vector (kprim);
 420
 421    iprim->renormer ();
 422    kprim->renormer ();
 423    jprim->vectoriel (kprim, iprim);
 424
 425    // transfoVertices (cross_center, iprim, jprim, kprim);
 426
 427    // Real3 snorm, bnorm;
 428    // iprim->getCoord (snorm);
 429    // kprim->getCoord (bnorm);
 430
 431    Real3 bnorm  = {0, 0, 1};
 432    Real3 snorm  = {1, 0, 0};
 433    assoCylinders (snorm, bnorm);
 434
 435 /*********************************************
 436    if (at_left)
 437       {
 438       assoIntersection (NxExt, 1, snorm, bnorm);
 439       if (grid_type == GR_BIPIPE)
 440          {
 441          assoIntersection (NxInt, 2, snorm, bnorm);
 442          }
 443       }
 444
 445    if (at_right)
 446       {
 447       assoIntersection (NxExt, NbrSlices1-1, snorm, bnorm);
 448       if (grid_type == GR_BIPIPE)
 449          {
 450          assoIntersection (NxInt, NbrSlices1-2, snorm, bnorm);
 451          }
 452       }
 453
 454   ******************************************* */
 455    assoResiduelle ();
 456    return HOK;
 457 }
 458 // ====================================================== createLittleCyl
 459 void BiCylinder::createLittleCyl ()
 460 {
 461    Real3   base;
 462    Vertex* vbase = cross_cyl1->getBase();
 463
 464    Real lg   =  cross_hauteur[CylSmall];
 465    Real h1   =  calc_distance (cross_center, vbase->getPoint (base));
 466    Real h2   =  cross_rayext[CylBig]*cos45;
 467    Real h3   = (cross_rayext[CylBig] - cross_rayint[CylBig])*cos45;
 468
 469    double t_haut [NbrVslices] = { -h1, -h2, -h2+h3, h2-h3, h2, lg-h1 };
 470
 471    for (int nk=0; nk<NbrVslices ; nk++)
 472        {
 473        switch (nk)
 474           {
 475           case 0 : case 5 :
 476               addSlice (CylSmall, 0, nk, t_haut[nk], cross_rayint [CylSmall]);
 477               addSlice (CylSmall, 1, nk, t_haut[nk], cross_rayext [CylSmall]);
 478               break;
 479           case 1 : case 4 :
 480               addSlice (CylSmall, 1, nk, t_haut[nk], cross_rayext [CylSmall]);
 481               addSlice (CylSmall, 2, nk, t_haut[nk], cross_rayext [CylBig]);
 482               break;
 483           case 2 : case 3 :
 484               addSlice (CylSmall, 0, nk, t_haut[nk], cross_rayint [CylSmall]);
 485               addSlice (CylSmall, 1, nk, t_haut[nk], cross_rayint [CylBig]);
 486               break;
 487           }
 488        }
 489                     //    ka kb kc kd
 490    fillSlice (CylSmall, 0,0, 0,2, 1,1, 1,0);
 491    fillSlice (CylSmall, 0,2, 1,2, 2,1, 1,1);  // OK
 492    // fillSlice (CylSmall, 1,1, 0,2, 1,2, 2,1);     // Test
 493    fillSlice (CylSmall, 1,2, 1,3, 2,4, 2,1, true);
 494    fillSlice (CylSmall, 1,3, 0,3, 1,4, 2,4);
 495    fillSlice (CylSmall, 0,3, 0,5, 1,5, 1,4);
 496    // fillSmallCyl ();
 497 }
 498 // ====================================================== createBigCyl
 499 void BiCylinder::createBigCyl ()
 500 {
 501    const int cyl = CylBig;
 502    Real3   base;
 503    Vertex* vbase = cross_cyl2->getBase();
 504    Real lg   = cross_hauteur[cyl];
 505    Real rext = cross_rayext [cyl];
 506    Real rint = cross_rayint [cyl];
 507    Real h1   = calc_distance (cross_center, vbase->getPoint (base));
 508    Real h2   = rext * cos45;
 509    Real h3   = lg - h1;
 510    Real dh   = (rext - rint)*cos45;
 511
 512    addSlice (CylBig, 0, 0, -h1,    rint);
 513    addSlice (CylBig, 1, 0, -h1,    rext);
 514    addSlice (CylBig, 0, 1, -h2+dh, rint, true);
 515    addSlice (CylBig, 1, 1, -h2,    rext, true);
 516    addSlice (CylBig, 0, 2,  h2-dh, rint, true);
 517    addSlice (CylBig, 1, 2,  h2,    rext, true);
 518    addSlice (CylBig, 0, 3,  h3,    rint);
 519    addSlice (CylBig, 1, 3,  h3,    rext);
 520
 521                   //   A    B    C    D
 522    fillSlice (CylBig, 0,0, 0,1, 1,1, 1,0);
 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 // ====================================================== fillSlice
 574 /* *****************************************************************
 575        H=bed  +----bd-----+ bdf=G
 576              /|          /|
 577            be |   B    bf |
 578            /  |        /  |
 579     E=bce +----bc-----+...|...bcf=F
 580           |  de     D |   df
 581           | E |       | F |             J
 582          ce   | C     cf  |             ^
 583   D=ade...|...+----ad-|---+ adf=C       |   I
 584           |  /        |  /              |  /
 585           | ae    A   | af              | /
 586           |/          |/                |/
 587     A=ace +----ac-----+ acf=B           +----->  K
 588
 589    ***************************************************************** */
 590 void BiCylinder::fillSlice (int cyl, int nia, int nka, int nib, int nkb,
 591                                      int nic, int nkc, int nid, int nkd,
 592                                      bool meddle)
 593 {
 594    for (int nj0=0; nj0<NbrCotes ; nj0++)
 595        {
 596        if (meddle) nj0++;
 597        int nj1 = nj0+1;
 598        if (nj1>=NbrCotes) nj1=0;
 599        Vertex* vace = getVertexIJK (cyl, nia, nj0, nka);
 600        Vertex* vacf = getVertexIJK (cyl, nib, nj0, nkb);
 601        Vertex* vadf = getVertexIJK (cyl, nic, nj0, nkc);
 602        Vertex* vade = getVertexIJK (cyl, nid, nj0, nkd);
 603
 604        Vertex* vbce = getVertexIJK (cyl, nia, nj1, nka);
 605        Vertex* vbcf = getVertexIJK (cyl, nib, nj1, nkb);
 606        Vertex* vbdf = getVertexIJK (cyl, nic, nj1, nkc);
 607        Vertex* vbde = getVertexIJK (cyl, nid, nj1, nkd);
 608
 609 /* *******************
 610        PutName (vace);
 611        PutName (vacf);
 612        PutName (vadf);
 613        PutName (vade);
 614        PutName (vbce);
 615        PutName (vbcf);
 616        PutName (vbdf);
 617        PutName (vbde);
 618    ******************* */
 619
 620        Edge* eac = addEdge (vace, vacf, cyl, dir_z, nia, nj0, nka);
 621        Edge* ead = addEdge (vade, vadf, cyl, dir_z, nid, nj0, nkd);
 622        Edge* eae = addEdge (vace, vade, cyl, dir_x, nia, nj0, nka);
 623        Edge* eaf = addEdge (vacf, vadf, cyl, dir_x, nib, nj0, nkb);
 624
 625        Edge* ebc = addEdge (vbce, vbcf, cyl, dir_z, nia, nj1, nka);
 626        Edge* ebd = addEdge (vbde, vbdf, cyl, dir_z, nid, nj1, nkd);
 627        Edge* ebe = addEdge (vbce, vbde, cyl, dir_x, nia, nj1, nka);
 628        Edge* ebf = addEdge (vbcf, vbdf, cyl, dir_x, nib, nj1, nkb);
 629
 630        Edge* ece = addEdge (vace, vbce, cyl, dir_y, nia, nj0, nka);
 631        Edge* ecf = addEdge (vacf, vbcf, cyl, dir_y, nib, nj0, nkb);
 632        Edge* edf = addEdge (vadf, vbdf, cyl, dir_y, nic, nj0, nkc);
 633        Edge* ede = addEdge (vade, vbde, cyl, dir_y, nid, nj0, nkd);
 634
 635        Quad* qa  = addQuad (eac, eaf, ead, eae, cyl, dir_y, nia , nj0, nka);
 636        Quad* qb  = addQuad (ebc, ebf, ebd, ebe, cyl, dir_y, nia , nj1, nka);
 637        Quad* qc  = addQuad (eac, ecf, ebc, ece, cyl, dir_x, nia , nj0, nka);
 638        Quad* qd  = addQuad (ead, edf, ebd, ede, cyl, dir_x, nid , nj0, nkd);
 639        Quad* qe  = addQuad (eae, ede, ebe, ece, cyl, dir_z, nia , nj0, nka);
 640        Quad* qf  = addQuad (eaf, edf, ebf, ecf, cyl, dir_z, nib , nj0, nkb);
 641
 642        addHexa (qa, qb, qc, qd, qe, qf, cyl, nia, nj0, nka);
 643        }
 644 }
 645 // ====================================================== assoCylinders
 646 void BiCylinder::assoCylinders (double* snormal, double* gnormal)
 647 {
 648    assoSlice (CylSmall, 0, 0, snormal);
 649    assoSlice (CylSmall, 1, 0, snormal);
 650    assoSlice (CylSmall, 0, 5, snormal);
 651    assoSlice (CylSmall, 1, 5, snormal);
 652
 653    for (int nk=0 ; nk<4 ; nk++)
 654        for (int ni=0 ; ni<2 ; ni++)
 655             assoSlice (CylBig, ni, nk, gnormal);
 656
 657    assoIntersection (1, 1, snormal, gnormal);
 658    assoIntersection (0, 2, snormal, gnormal);
 659    assoIntersection (0, 3, snormal, gnormal);
 660    assoIntersection (1, 4, snormal, gnormal);
 661
 662 }
 663 // ====================================================== assoSlice
 664 void BiCylinder::assoSlice (int cyl, int nx, int nzs, double* normal)
 665 {
 666    Real3  center, pnt1, pnt2, vbase;
 667    string brep;
 668    int ny0 = 0;
 669    int nyd = NbrCotes/2;
 670
 671    Vertex* v0 = getVertexIJK (cyl, nx, ny0 , nzs);
 672    Vertex* vd = getVertexIJK (cyl, nx, nyd , nzs);
 673
 674    if (vd==NULL || v0==NULL)
 675       return;
 676
 677    v0->getPoint (pnt1);
 678    vd->getPoint (pnt2);
 679
 680    double rayon = 0;
 681    for (int nro=0 ; nro<DIM3 ; nro++)
 682        {
 683        center[nro] = (pnt1[nro] + pnt2[nro])/2;
 684        vbase [nro] = (pnt1[nro] - center[nro]);
 685        rayon += vbase [nro]*vbase[nro];
 686        }
 687    rayon = sqrt (rayon);
 688
 689    PutCoord (pnt1);
 690    PutCoord (pnt2);
 691    PutCoord (vbase);
 692    PutCoord (center);
 693    PutData  (rayon);
 694
 695    geom_create_circle (center, rayon, normal, vbase, brep);
 696    geom_define_line   (brep);   // pour geom_asso_point
 697
 698    for (int ny=0 ; ny<NbrCotes ; ny++)
 699        {
 700        assoArc (cyl, nx, ny, nzs, brep, rayon);
 701        }
 702 }
 703 // ===================================================== assoArc
 704 void BiCylinder::assoArc (int cyl, int nx, int ny, int nz, string& brep,
 705                              double rayon)
 706 {
 707     static const double Alpha = 1.0 / NbrCotes;
 708     static const double Theta = 2*M_PI/ NbrCotes;
 709     int ny1 = (ny+1) MODULO NbrCotes;
 710
 711     Edge*   edge  = getEdgeJ  (cyl, nx, ny, nz);
 712     Vertex* node0 = getVertexIJK (cyl, nx, ny, nz);
 713     Vertex* node1 = getVertexIJK (cyl, nx, ny1, nz);
 714
 715     if (db)
 716        printf ("AssoArc : Edge = %s  = (%s,%s) -> (%s,%s)\n", edge->getName(),
 717                edge->getVertex(0)->getName(), edge->getVertex(1)->getName(),
 718                node0->getName(), node1->getName());
 719
 720     // Shape*  shape = new Shape (brep);
 721
 722     // shape->setBounds (ny*Alpha, (ny+1)*Alpha);
 723     // edge ->addAssociation (shape);
 724
 725     geom_asso_point ( ny   *Theta*rayon, node0);
 726     geom_asso_point ((ny+1)*Theta*rayon, node1);
 727 }
 728 // ===================================================== assoIntersection
 729 int BiCylinder::assoIntersection (int nxs, int nzs, double* snorm,
 730                                                     double* bnorm)
 731 {
 732    Real3  X_center = {0, 0, 0};  // provisoire, a la place de cross_center
 733    Real3  pse, psw, sorig, sbase;
 734    Real3  pbe, pbw, borig, bbase;
 735    string brep;
 736    int ny0 = 0;
 737    int nyd = NbrCotes/2;
 738    int MiddleSlice1 = 3;
 739
 740    int nz = nzs < MiddleSlice1 ? 0 : 5;
 741
 742    getVertexIJK (CylSmall, nxs, ny0 , nz)->getPoint (pse);
 743    getVertexIJK (CylSmall, nxs, nyd , nz)->getPoint (psw);
 744    getVertexIJK (CylBig,   nxs, ny0 , 0) ->getPoint (pbe);
 745    getVertexIJK (CylBig,   nxs, nyd , 0) ->getPoint (pbw);
 746
 747    double srayon = calc_distance (pse, psw)/2;
 748    double brayon = calc_distance (pbe, pbw)/2;
 749
 750    calc_milieu  (psw, pse, sorig);
 751    calc_milieu  (pbw, pbe, borig);
 752    calc_vecteur (psw, pse, sbase);
 753    calc_vecteur (pbw, pbe, bbase);
 754
 755    double shaut = calc_distance (X_center, sorig);
 756    double bhaut = calc_distance (X_center, borig)*2;
 757    double* orig = nzs < MiddleSlice1 ? sorig : X_center; // Pb orientation
 758
 759    if (db)
 760       {
 761       PutCoord  (borig);
 762       PutCoord  (sorig);
 763       PutCoord  (orig);
 764       PutData (nz);
 765       PutData  (srayon);
 766       PutData  (brayon);
 767       }
 768
 769    int ier = geom_create_cylcyl (borig, bnorm, bbase, brayon, bhaut,
 770                                   orig, snorm, sbase, srayon, shaut);
 771    if (ier != HOK)
 772       return ier;
 773
 774    for (int ny=0 ; ny<NbrCotes ; ny++)
 775        {
 776        Vertex* node = getVertexIJK (CylSmall, nxs, ny, nzs);
 777        if (node!=NULL)
 778            node->clearAssociation ();
 779        }
 780
 781    for (int ny=0 ; ny<NbrCotes ; ny++)
 782        {
 783        Edge* edge = getEdgeJ (CylSmall, nxs, ny, nzs);
 784        geom_asso_cylcyl (edge);
 785        }
 786
 787    return HOK;
 788 }
 789 // ======================================================== test_bicylinder
 790 BiCylinder* test_bicylinder (Document* docu, int option)
 791 {
 792    Vertex* ori1 = docu->addVertex ( 0,0,0);
 793    Vertex* ori2 = docu->addVertex (-5,0,5);
 794    Vector* vz   = docu->addVector ( 0,0,1);
 795    Vector* vx   = docu->addVector ( 1,0,0);
 796
 797    double r1 = 2*sqrt (2.0);
 798    double r2 = 3*sqrt (2.0);
 799    double l2 = 10;
 800    double l1 = 10;
 801
 802    Cylinder* cyl1 = docu->addCylinder (ori1, vz, r1, l1);
 803    Cylinder* cyl2 = docu->addCylinder (ori2, vx, r2, l2);
 804
 805    BiCylinder* grid = new BiCylinder (docu);
 806    grid->crossCylinders (cyl1, cyl2);
 807    return grid;
 808 }
 809 END_NAMESPACE_HEXA