src/HEXABLOCK/HexElements_bis.cxx

   1
   2 // C++ : Table d'hexaedres
   3
   4 // Copyright (C) 2009-2023  CEA, EDF
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include "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 "HexNewShape.hxx"
  33
  34 #include <map>
  35
  36 BEGIN_NAMESPACE_HEXA
  37
  38 // ====================================================== getHexaIJK
  39 Hexa* Elements::getHexaIJK (int nx, int ny, int nz)
  40 {
  41    if (isBad())
  42       return NULL;
  43    if (nx<0 || nx>=size_hx ||  ny<0 || ny>=size_hy || nz<0 || nz>=size_hz)
  44       return NULL;
  45    else if (grid_nocart)
  46       return NULL;
  47
  48    int nro = nx + size_hx*ny + size_hx*size_hy*nz;
  49
  50    DumpStart  ("getHexaIJK", nx << ny << nz);
  51    DumpReturn (tab_hexa [nro]);
  52    return      tab_hexa [nro];
  53 }
  54 // ====================================================== getQuadIJ
  55 Quad* Elements::getQuadIJ (int nx, int ny, int nz)
  56 {
  57    if (isBad())
  58       return NULL;
  59    if (nx<0 || nx>=size_qx ||  ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
  60       return NULL;
  61    else if (grid_nocart)
  62       return NULL;
  63
  64    int nro = nx + size_qx*ny + size_qx*size_qy*nz
  65                 + size_qx*size_qy*size_qz*dir_z;
  66
  67    DumpStart  ("getQuadIJ", nx << ny << nz);
  68    DumpReturn (tab_quad [nro]);
  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    DumpStart  ("getQuadJK", nx << ny << nz);
  84    DumpReturn (tab_quad [nro]);
  85    return tab_quad [nro];
  86 }
  87 // ====================================================== getQuadIK
  88 Quad* Elements::getQuadIK (int nx, int ny, int nz)
  89 {
  90    if (isBad())
  91       return NULL;
  92    if (nx<0 || nx>=size_qx ||  ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
  93       return NULL;
  94    else if (grid_nocart)
  95       return NULL;
  96
  97    int nro = nx + size_qx*ny + size_qx*size_qy*nz + size_qx*size_qy*size_qz;
  98
  99    DumpStart  ("getQuadIK", nx << ny << nz);
 100    DumpReturn (tab_quad [nro]);
 101    return tab_quad [nro];
 102 }
 103 // ====================================================== getEdgeI
 104 Edge* Elements::getEdgeI (int nx, int ny, int nz)
 105 {
 106    if (isBad())
 107       return NULL;
 108    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
 109       return NULL;
 110    else if (grid_nocart)
 111       return NULL;
 112
 113    int nro = nx + size_ex*ny + size_ex*size_ey*nz;
 114
 115    DumpStart  ("getEdgeI", nx << ny << nz);
 116    DumpReturn (tab_edge [nro]);
 117    return tab_edge [nro];
 118 }
 119 // ====================================================== getEdgeJ
 120 Edge* Elements::getEdgeJ (int nx, int ny, int nz)
 121 {
 122    if (isBad())
 123       return NULL;
 124    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
 125       return NULL;
 126    else if (grid_nocart)
 127       return NULL;
 128
 129    int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez;
 130
 131    DumpStart  ("getEdgeJ", nx << ny << nz);
 132    DumpReturn (tab_edge [nro]);
 133    return tab_edge [nro];
 134 }
 135 // ====================================================== getEdgeK
 136 Edge* Elements::getEdgeK (int nx, int ny, int nz)
 137 {
 138    if (isBad())
 139       return NULL;
 140    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
 141       return NULL;
 142    else if (grid_nocart)
 143       return NULL;
 144
 145    int nro = nx + size_ex*ny + size_ex*size_ey*nz
 146                 + size_ex*size_ey*size_ez*dir_z;
 147
 148    DumpStart  ("getEdgeK", nx << ny << nz);
 149    DumpReturn (tab_edge [nro]);
 150    return tab_edge [nro];
 151 }
 152 // ====================================================== getVertexIJK
 153 Vertex* Elements::getVertexIJK (int nx, int ny, int nz)
 154 {
 155    if (isBad())
 156       return NULL;
 157    if (nx<0 || nx>=size_vx ||  ny<0 || ny>=size_vy || nz<0 || nz>=size_vz)
 158       return NULL;
 159    else if (grid_nocart)
 160       return NULL;
 161
 162    int nro = nx + size_vx*ny + size_vx*size_vy*nz;
 163
 164    DumpStart  ("getVertexIJK", nx << ny << nz);
 165    DumpReturn (tab_vertex [nro]);
 166    return tab_vertex [nro];
 167 }
 168 // ====================================================== setVertex
 169 void Elements::setVertex (Vertex* elt, int nx, int ny, int nz)
 170 {
 171    if (   nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy
 172        || nz < 0 || nz >= size_vz) return;
 173
 174    int nro = nx + size_vx*ny + size_vx*size_vy*nz;
 175    tab_vertex [nro] = elt;
 176 }
 177 // ====================================================== setVertex (2)
 178 void Elements::setVertex (int nx, int ny, int nz, double px, double py,
 179                                                   double pz)
 180 {
 181    if (   nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy
 182        || nz < 0 || nz >= size_vz) return;
 183
 184    Vertex*    node = el_root->addVertex (px, py, pz);
 185    setVertex (node, nx, ny, nz);
 186 }
 187 // ====================================================== setEdge
 188 void Elements::setEdge (Edge* elt, EnumCoord dir, int nx, int ny, int nz)
 189 {
 190    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez
 191             || dir < dir_x || dir > dir_z )
 192       return;
 193
 194    int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez*dir;
 195    tab_edge [nro] = elt;
 196 }
 197 // ====================================================== setQuad
 198 void Elements::setQuad (Quad* elt, EnumCoord dir, int nx, int ny, int nz)
 199 {
 200    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez
 201             || dir < dir_x || dir > dir_z )
 202       return;
 203
 204    int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez*dir;
 205    tab_quad [nro] = elt;
 206 }
 207 // ====================================================== setHexa
 208 void Elements::setHexa (Hexa* elt, int nx, int ny, int nz)
 209 {
 210    if (   nx < 0 || nx >= size_hx || ny < 0 || ny >= size_hy
 211        || nz < 0 || nz >= size_hz) return;
 212
 213    int nro = nx + size_hx*ny + size_hx*size_hy*nz;
 214    tab_hexa [nro] = elt;
 215 }
 216 // ====================================================== remove
 217 void Elements::remove ()
 218 {
 219    int nbre=tab_hexa.size ();
 220    nbre = nbr_hexas;
 221    for (int nh=0 ; nh<nbre ; nh++)
 222        if (tab_hexa[nh] != NULL)
 223            tab_hexa[nh]->remove();
 224 }
 225 // ====================================================== prismHexas
 226 int  Elements::prismHexas (int nro, Quad* qbase, int hauteur)
 227 {
 228    int   ind_node [QUAD4];
 229    int   subid = 0;
 230    Real3 koord, transfo;
 231
 232            // ----------------------------- Vertex + aretes verticales
 233    for (int ns=0 ; ns<QUAD4 ; ns++)
 234        {
 235        Vertex* vbase = qbase->getVertex (ns);
 236        int     indx  = vbase->getMark ();
 237        if (indx<0)
 238           {
 239           bool asso = vbase->isAssociated();
 240           vbase->getAssoCoord (koord);
 241
 242           indx = nbr_vertex++;
 243           vbase->setMark (indx);
 244           Vertex* nd0 = vbase;
 245           Vertex* nd1 = NULL;
 246           double beta = 0;
 247           if (revo_lution)
 248              {
 249              Real3 centre, point, om;
 250              revo_center->getPoint (centre);
 251              vbase      ->getPoint (point);
 252
 253              calc_vecteur  (centre, point, om);
 254              double oh    = prod_scalaire (om, revo_axe);
 255              double rayon = 0;
 256              Real3  ph, hm;
 257              for (int dd=dir_x; dd<=dir_z ; dd++)
 258                  {
 259                  ph [dd] = centre [dd] + oh*revo_axe[dd];
 260                  hm [dd] = point  [dd] - ph[dd];
 261                  rayon  += hm[dd] * hm[dd];
 262                  }
 263              subid = grid_geom->addCircle (ph, sqrt(rayon), revo_axe, hm);
 264              }
 265
 266           for (int nh=0 ; nh<hauteur ; nh++)
 267               {
 268               nd1 = el_root->addVertex (nd0->getX(), nd0->getY(), nd0->getZ());
 269               updateMatrix (nh);
 270               gen_matrix.perform   (nd1);
 271               tab_vertex.push_back (nd1);
 272               Edge* pilier = newEdge (nd0, nd1);
 273               tab_pilier.push_back (pilier);
 274               if (revo_lution)
 275                  {
 276                  double alpha = beta;
 277                  beta = alpha + gen_values[nh];
 278                  grid_geom->addAssociation (pilier, subid, alpha/360, beta/360);
 279                  }
 280               if (asso)
 281                  {
 282                  cum_matrix.perform  (koord, transfo);
 283                  nd1->setAssociation (transfo);
 284                  }
 285               nd0 = nd1;
 286               }
 287           }
 288        ind_node [ns] = indx;
 289        }
 290            // ----------------------------- Aretes horizontales
 291            // ----------------------------- + face verticales
 292    int ind_poutre [QUAD4];
 293    for (int ns=0 ; ns<QUAD4 ; ns++)
 294        {
 295        Edge* ebase = qbase->getEdge (ns);
 296        int   indx  = ebase->getMark ();
 297        if (indx<0)
 298           {
 299           indx = nbr_edges ++;
 300           ebase->setMark (indx);
 301           int   nd1 = ind_node [ns];
 302           int   nd2 = ind_node [(ns+1) MODULO QUAD4];
 303           Edge* ed0 = ebase;
 304           Edge *ed1, *ed2, *ed3;
 305           for (int nh=0 ; nh<hauteur ; nh++)
 306               {
 307               ed2 = ed0;
 308               ed0 = newEdge (tab_vertex [nd1*hauteur + nh],
 309                              tab_vertex [nd2*hauteur + nh]);
 310               ed1 = tab_pilier [nd1*hauteur + nh];
 311               ed3 = tab_pilier [nd2*hauteur + nh];
 312
 313               Quad* mur = newQuad (ed0, ed1, ed2, ed3);
 314               tab_edge.push_back (ed0);
 315               tab_quad.push_back (mur);
 316               if (ebase->isAssociated ())
 317                   prismAssociation (ebase, ed0, nh);
 318               }
 319           }
 320        ind_poutre [ns] = indx;
 321        }
 322            // ----------------------------- Faces horizontales
 323            // ----------------------------- + Hexas
 324    Quad* qa = qbase;
 325    Quad *qb, *qc, *qd, *qe, *qf;
 326    int nv0 = hauteur*ind_poutre [0];
 327    int nv1 = hauteur*ind_poutre [1];
 328    int nv2 = hauteur*ind_poutre [2];
 329    int nv3 = hauteur*ind_poutre [3];
 330    for (int nh=0 ; nh<hauteur ; nh++)
 331        {
 332        qb = newQuad (tab_edge [nh+nv0], tab_edge [nh+nv1],
 333                      tab_edge [nh+nv2], tab_edge [nh+nv3]);
 334        qc = tab_quad [nh + nv0];
 335        qd = tab_quad [nh + nv2];
 336        qe = tab_quad [nh + nv1];
 337        qf = tab_quad [nh + nv3];
 338
 339 // *** tab_hexa [nh*hauteur + nro] = newHexa (qa, qb, qc, qd, qe, qf); Abu
 340        tab_hexa [nro*hauteur + nh] = newHexa (qa, qb, qc, qd, qe, qf);
 341        ker_hquad.push_back (qb);
 342        qa = qb;
 343        }
 344    return HOK;
 345 }
 346 // ====================================================== updateMatrix
 347 void Elements::updateMatrix (int nh)
 348 {
 349    if (revo_lution)
 350       {
 351       gen_matrix.defRotation (revo_center, revo_axe, gen_values[nh]);
 352       cum_matrix.defRotation (revo_center, revo_axe, cum_values[nh]);
 353       }
 354    else
 355       {
 356       double hauteur = (nh+1)*prism_len;
 357       double dh      = prism_len;
 358       if (prism_vec)
 359          {
 360          if (under_v6)
 361             {
 362             hauteur = gen_values[nh];
 363             dh      = nh>0 ? hauteur-gen_values[nh-1] : hauteur;
 364             }
 365          else
 366             {
 367             hauteur = cum_values[nh];
 368             dh      = gen_values[nh];
 369             }
 370          }
 371       Real3 trans, decal;
 372       for (int nc=dir_x ; nc<=dir_z ; nc++)
 373           {
 374           decal [nc] = prism_dir [nc]*dh;
 375           trans [nc] = prism_dir [nc]*hauteur;
 376           }
 377
 378       gen_matrix.defTranslation (decal);
 379       cum_matrix.defTranslation (trans);
 380       }
 381 }
 382 // ====================================================== endPrism
 383 void Elements::endPrism ()
 384 {
 385    closeShape();
 386
 387    int nbelts = ker_hquad.size();
 388    for (int nro=0 ; nro<nbelts ; nro++)
 389        tab_quad.push_back (ker_hquad[nro]);
 390
 391    nbelts = tab_pilier.size();
 392    for (int nro=0 ; nro<nbelts ; nro++)
 393        tab_edge.push_back (tab_pilier[nro]);
 394
 395
 396    nbr_hexas  = tab_hexa.size ();
 397    nbr_edges  = tab_edge.size ();
 398    nbr_quads  = tab_quad.size ();
 399    nbr_vertex = tab_vertex.size ();
 400 }
 401 // ====================================================== getShape
 402 NewShape* Elements::getShape()
 403 {
 404    if (grid_geom==NULL)
 405       {
 406       std::string name = "extrud_" + el_name;
 407       grid_geom   = el_root->addShape (name.c_str(), SH_EXTRUD);
 408       grid_geom -> openShape();
 409       }
 410
 411    return grid_geom;
 412 }
 413 // ====================================================== closeShape
 414 void Elements::closeShape()
 415 {
 416    if (grid_geom==NULL)
 417        return;
 418
 419    grid_geom -> closeShape();
 420    grid_geom = NULL;
 421 }
 422 END_NAMESPACE_HEXA