src/HEXABLOCK/HexElements_bis.cxx

   1
   2 // C++ : Table d'hexaedres
   3
   4 // Copyright (C) 2009-2012  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 "HexElements.hxx"
  24
  25 #include "HexDocument.hxx"
  26 #include "HexVector.hxx"
  27 #include "HexVertex.hxx"
  28 #include "HexHexa.hxx"
  29 #include "HexEdge.hxx"
  30
  31 #include "HexGlobale.hxx"
  32 #include "HexCylinder.hxx"
  33 #include "HexShape.hxx"
  34
  35 #include <map>
  36
  37 BEGIN_NAMESPACE_HEXA
  38
  39 void geom_dump_asso     (Edge* edge);
  40 void geom_create_circle (double* milieu, double rayon, double* normale,
  41                          double* base, string& brep);
  42
  43 // ====================================================== getHexaIJK
  44 Hexa* Elements::getHexaIJK (int nx, int ny, int nz)
  45 {
  46    if (isBad())
  47       return NULL;
  48    if (nx<0 || nx>=size_hx ||  ny<0 || ny>=size_hy || nz<0 || nz>=size_hz)
  49       return NULL;
  50    else if (grid_nocart)
  51       return NULL;
  52
  53    int nro = nx + size_hx*ny + size_hx*size_hy*nz;
  54
  55    return tab_hexa [nro];
  56 }
  57 // ====================================================== getQuadIJ
  58 Quad* Elements::getQuadIJ (int nx, int ny, int nz)
  59 {
  60    if (isBad())
  61       return NULL;
  62    if (nx<0 || nx>=size_qx ||  ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
  63       return NULL;
  64    else if (grid_nocart)
  65       return NULL;
  66
  67    int nro = nx + size_qx*ny + size_qx*size_qy*nz
  68                 + size_qx*size_qy*size_qz*dir_z;
  69    return tab_quad [nro];
  70 }
  71 // ====================================================== getQuadJK
  72 Quad* Elements::getQuadJK (int nx, int ny, int nz)
  73 {
  74    if (isBad())
  75       return NULL;
  76    if (nx<0 || nx>=size_qx ||  ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
  77       return NULL;
  78    else if (grid_nocart)
  79       return NULL;
  80
  81    int nro = nx + size_qx*ny + size_qx*size_qy*nz; // + dir_x*...
  82
  83    return tab_quad [nro];
  84 }
  85 // ====================================================== getQuadIK
  86 Quad* Elements::getQuadIK (int nx, int ny, int nz)
  87 {
  88    if (isBad())
  89       return NULL;
  90    if (nx<0 || nx>=size_qx ||  ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
  91       return NULL;
  92    else if (grid_nocart)
  93       return NULL;
  94
  95    int nro = nx + size_qx*ny + size_qx*size_qy*nz + size_qx*size_qy*size_qz;
  96
  97    return tab_quad [nro];
  98 }
  99 // ====================================================== getEdgeI
 100 Edge* Elements::getEdgeI (int nx, int ny, int nz)
 101 {
 102    if (isBad())
 103       return NULL;
 104    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
 105       return NULL;
 106    else if (grid_nocart)
 107       return NULL;
 108
 109    int nro = nx + size_ex*ny + size_ex*size_ey*nz;
 110
 111    return tab_edge [nro];
 112 }
 113 // ====================================================== getEdgeJ
 114 Edge* Elements::getEdgeJ (int nx, int ny, int nz)
 115 {
 116    if (isBad())
 117       return NULL;
 118    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
 119       return NULL;
 120    else if (grid_nocart)
 121       return NULL;
 122
 123    int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez;
 124
 125    return tab_edge [nro];
 126 }
 127 // ====================================================== getEdgeK
 128 Edge* Elements::getEdgeK (int nx, int ny, int nz)
 129 {
 130    if (isBad())
 131       return NULL;
 132    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
 133       return NULL;
 134    else if (grid_nocart)
 135       return NULL;
 136
 137    int nro = nx + size_ex*ny + size_ex*size_ey*nz
 138                 + size_ex*size_ey*size_ez*dir_z;
 139    return tab_edge [nro];
 140 }
 141 // ====================================================== getVertexIJK
 142 Vertex* Elements::getVertexIJK (int nx, int ny, int nz)
 143 {
 144    if (isBad())
 145       return NULL;
 146    if (nx<0 || nx>=size_vx ||  ny<0 || ny>=size_vy || nz<0 || nz>=size_vz)
 147       return NULL;
 148    else if (grid_nocart)
 149       return NULL;
 150
 151    int nro = nx + size_vx*ny + size_vx*size_vy*nz;
 152
 153    return tab_vertex [nro];
 154 }
 155 // ====================================================== setVertex
 156 void Elements::setVertex (Vertex* elt, int nx, int ny, int nz)
 157 {
 158    if (   nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy
 159        || nz < 0 || nz >= size_vz) return;
 160
 161    int nro = nx + size_vx*ny + size_vx*size_vy*nz;
 162    tab_vertex [nro] = elt;
 163 }
 164 // ====================================================== setVertex (2)
 165 void Elements::setVertex (int nx, int ny, int nz, double px, double py,
 166                                                   double pz)
 167 {
 168    if (   nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy
 169        || nz < 0 || nz >= size_vz) return;
 170
 171    Vertex*    node = el_root->addVertex (px, py, pz);
 172    setVertex (node, nx, ny, nz);
 173 }
 174 // ====================================================== setEdge
 175 void Elements::setEdge (Edge* elt, EnumCoord dir, int nx, int ny, int nz)
 176 {
 177    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez
 178             || dir < dir_x || dir > dir_z )
 179       return;
 180
 181    int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez*dir;
 182    tab_edge [nro] = elt;
 183 }
 184 // ====================================================== setQuad
 185 void Elements::setQuad (Quad* elt, EnumCoord dir, int nx, int ny, int nz)
 186 {
 187    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez
 188             || dir < dir_x || dir > dir_z )
 189       return;
 190
 191    int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez*dir;
 192    tab_quad [nro] = elt;
 193 }
 194 // ====================================================== setHexa
 195 void Elements::setHexa (Hexa* elt, int nx, int ny, int nz)
 196 {
 197    if (   nx < 0 || nx >= size_hx || ny < 0 || ny >= size_hy
 198        || nz < 0 || nz >= size_hz) return;
 199
 200    int nro = nx + size_hx*ny + size_hx*size_hy*nz;
 201    tab_hexa [nro] = elt;
 202 }
 203 // ====================================================== remove
 204 void Elements::remove ()
 205 {
 206    int nbre=tab_hexa.size ();
 207    nbre = nbr_hexas;
 208    for (int nh=0 ; nh<nbre ; nh++)
 209        if (tab_hexa[nh] != NULL)
 210            tab_hexa[nh]->remove();
 211 }
 212 // ====================================================== makeCylinder
 213 int Elements::makeCylinder (Cylinder* cyl, Vector* vx, int nr, int na, int nl)
 214 {
 215    if (BadElement (cyl) || BadElement (vx) ||  nr<=0 || na <=3 || nl <=0
 216                         || vx->getNorm () <= Epsil)
 217       {
 218       setError ();
 219       return HERR;
 220       }
 221
 222    Vertex* orig = cyl->getBase ();
 223    Vector* dir  = cyl->getDirection ();
 224    double  ray  = cyl->getRadius ();
 225    double  haut = cyl->getHeight ();
 226
 227    resize (GR_CYLINDRIC, nr, na, nl);
 228    cyl_closed = true;
 229    makeCylindricalNodes (orig, vx, dir, ray/(nr+1), 360, haut/nl,
 230                          nr, na, nl, true);
 231    fillGrid ();
 232    assoCylinder (orig, dir, 360);
 233    return HOK;
 234 }
 235 // ====================================================== makePipe
 236 int Elements::makePipe (Cylinder* cyl, Vector* vx, int nr, int na, int nl)
 237 {
 238    if (BadElement (cyl) || BadElement (vx) ||  nr<=0 || na <=3 || nl <=0
 239                         || vx->getNorm () <= Epsil)
 240       {
 241       setError ();
 242       return HERR;
 243       }
 244
 245    Vertex* orig = cyl->getBase ();
 246    Vector* dir  = cyl->getDirection ();
 247    double  ray  = cyl->getRadius ();
 248    double  haut = cyl->getHeight ();
 249
 250    resize (GR_CYLINDRIC, nr, na, nl);
 251    cyl_closed = true;
 252    makeCylindricalNodes (orig, vx, dir, ray, 360, haut, nr, na, nl, false);
 253    fillGrid ();
 254    assoCylinder (orig, dir, 360);
 255    return HOK;
 256 }
 257 //
 258 // ---------------------------------------- prism Quads
 259 //
 260 // ====================================================== prismQuads
 261 int Elements::prismQuads (Quads& tstart, Vector* dir, int nbiter)
 262 {
 263    if (BadElement (dir) || dir->getNorm () <= Epsil || nbiter <= 0)
 264       {
 265       setError ();
 266       return HERR;
 267       }
 268
 269    el_root->markAll (NO_USED);
 270    int nbcells   = tstart.size ();
 271    nbr_vertex    = 0;
 272    nbr_edges     = 0;
 273
 274    nbr_hexas = nbcells*nbiter;
 275
 276    tab_hexa.resize (nbr_hexas);
 277    tab_quad.clear ();          // verticaux
 278    ker_hquad.clear ();         // Horizontaux
 279    tab_edge.clear ();
 280    tab_pilier.clear ();
 281    tab_vertex.clear ();
 282
 283    revo_lution = false;
 284    prism_vec   = false;
 285    gen_matrix.defTranslation (dir);
 286
 287    for (int nro=0 ; nro<nbcells ; nro++)
 288        {
 289        prismHexas (nro, tstart[nro], nbiter);
 290        }
 291
 292    endPrism ();
 293    return HOK;
 294 }
 295 // ====================================================== prismQuadsVec
 296 int Elements::prismQuadsVec (Quads& tstart, Vector* dir, RealVector& tlen,
 297                              int mode)
 298 {
 299    int nbiter = tlen.size();
 300    if (BadElement (dir) || dir->getNorm () <= Epsil || nbiter <= 0)
 301       {
 302       setError ();
 303       return HERR;
 304       }
 305
 306    el_root->markAll (NO_USED);
 307    int nbcells   = tstart.size ();
 308    nbr_vertex    = 0;
 309    nbr_edges     = 0;
 310
 311    nbr_hexas = nbcells*nbiter;
 312
 313    tab_hexa.resize (nbr_hexas);
 314    tab_quad.clear ();          // verticaux
 315    ker_hquad.clear ();         // Horizontaux
 316    tab_edge.clear ();
 317    tab_pilier.clear ();
 318    tab_vertex.clear ();
 319
 320    revo_lution = false;
 321    prism_vec   = true;
 322    dir->getCoord  (prism_dir);
 323    normer_vecteur (prism_dir);
 324    gen_values = tlen;
 325
 326    for (int nro=0 ; nro<nbcells ; nro++)
 327        {
 328        prismHexas (nro, tstart[nro], nbiter);
 329        }
 330
 331    endPrism ();
 332    return HOK;
 333 }
 334 // ======================================================== revolutionQuads
 335 int Elements::revolutionQuads (Quads& start, Vertex* center, Vector* axis,
 336                                RealVector &angles)
 337 {
 338    int nbiter  = angles.size();
 339    int nbcells = start.size ();
 340    if (BadElement (center)  || BadElement(axis) || nbiter==0 || nbcells==0
 341                             || axis->getNorm () <= Epsil)
 342       {
 343       setError ();
 344       return HERR;
 345       }
 346
 347    el_root->markAll (NO_USED);
 348    nbr_vertex    = 0;
 349    nbr_edges     = 0;
 350
 351    nbr_hexas   = nbcells*nbiter;
 352
 353    tab_hexa.resize (nbr_hexas);
 354    tab_quad.clear ();          // verticaux
 355    ker_hquad.clear ();         // Horizontaux
 356    tab_edge.clear ();
 357    tab_pilier.clear ();
 358    tab_vertex.clear ();
 359
 360    revo_lution  = true;
 361    prism_vec    = false;
 362    revo_axis    = axis;
 363    revo_center  = center;
 364    gen_values = angles;
 365
 366    for (int nro=0 ; nro<nbcells ; nro++)
 367        {
 368        prismHexas (nro, start[nro], nbiter);
 369        }
 370
 371    endPrism ();
 372    return HOK;
 373 }
 374 // ====================================================== prismHexas
 375 int  Elements::prismHexas (int nro, Quad* qbase, int hauteur)
 376 {
 377    int ind_node [QUAD4];
 378    string c_rep;
 379
 380            // ----------------------------- Vertex + aretes verticales
 381    for (int ns=0 ; ns<QUAD4 ; ns++)
 382        {
 383        Vertex* vbase = qbase ->getVertex (ns);
 384        int     indx  = vbase->getMark ();
 385        if (indx<0)
 386           {
 387           indx = nbr_vertex++;
 388           vbase->setMark (indx);
 389           Vertex* nd0 = vbase;
 390           Vertex* nd1 = NULL;
 391           double beta = 0;
 392           if (revo_lution)
 393              {
 394              Real3 centre, vk, point, om;
 395              revo_center->getPoint (centre);
 396              vbase      ->getPoint (point);
 397              revo_axis  ->getCoord (vk);
 398              normer_vecteur (vk);
 399
 400              calc_vecteur   (centre, point, om);
 401              double oh     = prod_scalaire (om, vk);
 402              double rayon  = 0;
 403              Real3  ph, hm;
 404              for (int dd=dir_x; dd<=dir_z ; dd++)
 405                  {
 406                  ph [dd] = centre [dd] + oh*vk[dd];
 407                  hm [dd] = point  [dd] - ph[dd];
 408                  rayon  += hm[dd] * hm[dd];
 409                  }
 410              rayon = sqrt (rayon);
 411 /********************************
 412              PutCoord (centre);
 413              PutCoord (point);
 414              PutData  (oh);
 415              PutCoord (ph);
 416              PutData  (rayon);
 417              PutCoord (vk);
 418              PutCoord (hm);
 419 ********************************/
 420              geom_create_circle (ph, rayon, vk, hm, c_rep);
 421              }
 422
 423           for (int nh=0 ; nh<hauteur ; nh++)
 424               {
 425               nd1 = el_root->addVertex (nd0->getX(), nd0->getY(), nd0->getZ());
 426               updateMatrix (nh);
 427               gen_matrix.perform   (nd1);
 428               tab_vertex.push_back (nd1);
 429               Edge* pilier = newEdge (nd0, nd1);
 430               tab_pilier.push_back (pilier);
 431               if (revo_lution)
 432                  {
 433                  double alpha = beta;
 434                  beta = alpha + gen_values[nh];
 435                  Shape* shape = new Shape (c_rep);
 436                  shape->setBounds (alpha/360, beta/360);
 437                  pilier->addAssociation (shape);
 438                  //      geom_dump_asso (pilier);
 439                  }
 440               nd0 = nd1;
 441               }
 442           }
 443        ind_node [ns] = indx;
 444        }
 445            // ----------------------------- Aretes horizontales
 446            // ----------------------------- + face verticales
 447    int ind_poutre [QUAD4];
 448    for (int ns=0 ; ns<QUAD4 ; ns++)
 449        {
 450        Edge* ebase = qbase->getEdge (ns);
 451        int   indx  = ebase->getMark ();
 452        if (indx<0)
 453           {
 454           indx = nbr_edges ++;
 455           ebase->setMark (indx);
 456           int   nd1 = ind_node [ns];
 457           int   nd2 = ind_node [(ns+1) MODULO QUAD4];
 458           Edge* ed0 = ebase;
 459           Edge *ed1, *ed2, *ed3;
 460           for (int nh=0 ; nh<hauteur ; nh++)
 461               {
 462               ed2 = ed0;
 463               ed0 = newEdge (tab_vertex [nd1*hauteur + nh],
 464                              tab_vertex [nd2*hauteur + nh]);
 465               ed1 = tab_pilier [nd1*hauteur + nh];
 466               ed3 = tab_pilier [nd2*hauteur + nh];
 467
 468               Quad* mur = newQuad (ed0, ed1, ed2, ed3);
 469               tab_edge.push_back (ed0);
 470               tab_quad.push_back (mur);
 471               prismAssociation (ed2, ed0, nh, ed1);
 472               }
 473           }
 474        ind_poutre [ns] = indx;
 475        }
 476            // ----------------------------- Faces horizontales
 477            // ----------------------------- + Hexas
 478    Quad* qa = qbase;
 479    Quad *qb, *qc, *qd, *qe, *qf;
 480    int nv0 = hauteur*ind_poutre [0];
 481    int nv1 = hauteur*ind_poutre [1];
 482    int nv2 = hauteur*ind_poutre [2];
 483    int nv3 = hauteur*ind_poutre [3];
 484    for (int nh=0 ; nh<hauteur ; nh++)
 485        {
 486        qb = newQuad (tab_edge [nh+nv0], tab_edge [nh+nv1],
 487                      tab_edge [nh+nv2], tab_edge [nh+nv3]);
 488        qc = tab_quad [nh + nv0];
 489        qd = tab_quad [nh + nv2];
 490        qe = tab_quad [nh + nv1];
 491        qf = tab_quad [nh + nv3];
 492
 493 // *** tab_hexa [nh*hauteur + nro] = newHexa (qa, qb, qc, qd, qe, qf); Abu
 494        tab_hexa [nro*hauteur + nh] = newHexa (qa, qb, qc, qd, qe, qf);
 495        ker_hquad.push_back (qb);
 496        qa = qb;
 497        }
 498    return HOK;
 499 }
 500 // ====================================================== updateMatrix
 501 void Elements::updateMatrix (int hauteur)
 502 {
 503    if (revo_lution)
 504       {
 505       gen_matrix.defRotation (revo_center, revo_axis, gen_values[hauteur]);
 506       }
 507    else if (prism_vec)
 508       {
 509       double h0 = hauteur>0 ?  gen_values[hauteur-1] : 0;
 510       double dh = gen_values[hauteur] - h0;
 511       Real3 decal;
 512       for (int nc=dir_x ; nc<=dir_z ; nc++)
 513           decal [nc] = prism_dir [nc]*dh;
 514       gen_matrix.defTranslation (decal);
 515       }
 516 }
 517 // ====================================================== endPrism
 518 void Elements::endPrism ()
 519 {
 520    int nbelts = ker_hquad.size();
 521    for (int nro=0 ; nro<nbelts ; nro++)
 522        tab_quad.push_back (ker_hquad[nro]);
 523
 524    nbelts = tab_pilier.size();
 525    for (int nro=0 ; nro<nbelts ; nro++)
 526        tab_edge.push_back (tab_pilier[nro]);
 527
 528
 529    nbr_hexas  = tab_hexa.size ();
 530    nbr_edges  = tab_edge.size ();
 531    nbr_quads  = tab_quad.size ();
 532    nbr_vertex = tab_vertex.size ();
 533 }
 534 END_NAMESPACE_HEXA