src/HEXABLOCK/HexElements_bis.cxx

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