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