src/HEXABLOCK/HexElements_grid.cxx

   1 // Copyright (C) 2009-2023  CEA, EDF
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 // C++ : Table des noeuds
  21
  22 #include "HexElements.hxx"
  23
  24 #include "HexDocument.hxx"
  25 #include "HexVector.hxx"
  26 #include "HexVertex.hxx"
  27 #include "HexHexa.hxx"
  28 #include "HexEdge.hxx"
  29
  30 #include "HexGlobale.hxx"
  31
  32 #include <cmath>
  33
  34 BEGIN_NAMESPACE_HEXA
  35 // ====================================================== fillGrid
  36 int Elements::fillGrid ()
  37 {
  38    if (cyl_closed)
  39       for (int nx=0 ; nx<size_vx ; nx++)
  40            for (int nz=0 ; nz<size_vz ; nz++)
  41                setVertex (getVertexIJK (nx, 0, nz), nx, size_vy-1, nz);
  42
  43    for (int nz=0 ; nz<size_vz ; nz++)
  44        for (int ny=0 ; ny<size_vy ; ny++)
  45            for (int nx=0 ; nx<size_vx ; nx++)
  46                {
  47                Vertex* v0 = getVertexIJK (nx, ny, nz  );
  48                Vertex* vx = getVertexIJK (nx+1, ny, nz);
  49                Vertex* vy = getVertexIJK (nx, ny+1, nz);
  50                Vertex* vz = getVertexIJK (nx, ny, nz+1);
  51
  52                Edge* e1 = NULL;
  53                Edge* e2 = newEdge (v0, vy);
  54                Edge* e3 = NULL;
  55
  56                if (cyl_closed && ny==size_vy-1)
  57                   {
  58                   e1 = getEdgeI (nx, 0, nz);
  59                   e3 = getEdgeK (nx, 0, nz);
  60                   }
  61                else
  62                   {
  63                   e3 = newEdge (v0, vz);
  64                   e1 = newEdge (v0, vx);
  65                   }
  66
  67                setEdge (e1, dir_x, nx, ny, nz);
  68                setEdge (e2, dir_y, nx, ny, nz);
  69                setEdge (e3, dir_z, nx, ny, nz);
  70                }
  71
  72    for (int nz=0 ; nz<size_vz ; nz++)
  73        for (int ny=0 ; ny<size_vy ; ny++)
  74            for (int nx=0 ; nx<size_vx ; nx++)
  75                {
  76                Edge* eae = getEdgeI (nx, ny,   nz);
  77                Edge* ebe = getEdgeI (nx, ny,   nz+1);
  78                Edge* eaf = getEdgeI (nx, ny+1, nz);
  79
  80                Edge* eac = getEdgeJ (nx,   ny, nz);
  81                Edge* ead = getEdgeJ (nx+1, ny, nz);
  82                Edge* ebc = getEdgeJ (nx,   ny, nz+1);
  83
  84                Edge* ece = getEdgeK (nx,   ny,   nz);
  85                Edge* ede = getEdgeK (nx+1, ny,   nz);
  86                Edge* ecf = getEdgeK (nx,   ny+1, nz);
  87
  88                Quad* q1 = newQuad (eac, eaf, ead, eae);
  89                Quad* q2 = newQuad (eac, ecf, ebc, ece);
  90                Quad* q3 = NULL;
  91                Quad* q30 = getQuadIK (nx, 0 ,nz);
  92
  93                if (cyl_closed && ny==size_vy-1 && q30!= NULL)
  94                   {
  95                   q3 = q30;
  96                   // Display(q3);
  97                   }
  98                else
  99                   {
 100                   q3 = newQuad (eae, ede, ebe, ece);
 101                   }
 102
 103                setQuad (q1, dir_z, nx,ny,nz);
 104                setQuad (q2, dir_x, nx,ny,nz);
 105                setQuad (q3, dir_y, nx,ny,nz);
 106                }
 107
 108    for (int nz=0 ; nz<size_hz ; nz++)
 109        for (int ny=0 ; ny<size_hy ; ny++)
 110            for (int nx=0 ; nx<size_hx ; nx++)
 111                {
 112                Quad* qa = getQuadIJ (nx, ny, nz);
 113                Quad* qb = getQuadIJ (nx, ny, nz+1);
 114
 115                Quad* qc = getQuadIK (nx, ny,   nz);
 116                Quad* qd = getQuadIK (nx, ny+1, nz);
 117
 118                Quad* qe = getQuadJK (nx,   ny, nz);
 119                Quad* qf = getQuadJK (nx+1, ny, nz);
 120
 121                setHexa (newHexa (qa, qb, qc, qd, qe, qf), nx, ny, nz);
 122                }
 123
 124    switch (cyl_dispo)
 125       {
 126       case CYL_CLOSED :
 127       case CYL_PEER   :  fillCenter ();
 128            break ;
 129       case CYL_CL4 :     fillCenter4 ();
 130            break ;
 131       case CYL_CL6 :     fillCenter6 ();
 132            break ;
 133       case CYL_ODD :     fillCenterOdd ();
 134            break ;
 135       case CYL_NOFILL  :
 136       default : ;
 137       }
 138    return HOK;
 139 }
 140 // ====================================================== fillCenter
 141 // === Remplissage radial
 142 #define IndElt(nc,nz)   (nbsecteurs*(nz) + nc)
 143 #define IndRedge(nc,nz) (nbrayons  *(nz) + nc)
 144 #define IndVquad(nc,nz) (nbrayons  *(nz-1) + nc)
 145 void Elements::fillCenter ()
 146 {
 147    int nx0 = 0;
 148    int nbsecteurs = size_hy / 2;
 149    int nbrayons   = cyl_closed ? nbsecteurs : nbsecteurs + 1;
 150
 151    size_hplus  = nbsecteurs * size_hz;
 152    size_qhplus = nbsecteurs * size_vz;
 153    size_qvplus = nbrayons   * size_hz;
 154    size_ehplus = nbrayons   * size_vz;
 155    size_evplus = size_hz;
 156
 157    ker_hexa .resize (size_hplus,  NULL);
 158    ker_hquad.resize (size_qhplus, NULL);
 159    ker_vquad.resize (size_qvplus, NULL);
 160    ker_hedge.resize (size_ehplus, NULL);
 161    ker_vedge.resize (size_evplus, NULL);
 162
 163    Vertex* pcenter = NULL;
 164
 165    for (int nz=0 ; nz<size_vz ; nz++)
 166        {
 167                                  //   Vertex central
 168        Vertex* center = getVertex (ker_vertex+nz);
 169                                  //   Edges horizontaux radiaux
 170        for (int nc=0 ; nc<nbrayons ; nc++)
 171            {
 172            Vertex* vv = getVertexIJK (nx0, 2*nc, nz);
 173            Edge*   edge = newEdge (center, vv);
 174            ker_hedge [IndRedge(nc,nz)] = edge;
 175            }
 176                                  //   Quads horizontaux
 177        for (int nc=0; nc<nbsecteurs ; nc++)
 178            {
 179            int nc1  = (nc + 1) MODULO nbrayons;
 180            Edge* e1 = ker_hedge [IndRedge (nc, nz)];
 181            Edge* e2 = getEdgeJ (nx0, 2*nc,  nz);
 182            Edge* e3 = getEdgeJ (nx0, 2*nc+1,nz);
 183            Edge* e4 = ker_hedge [IndRedge (nc1, nz)];
 184
 185            ker_hquad [IndElt (nc,nz)] = newQuad (e1, e2, e3, e4);
 186            if (debug())
 187               {
 188               printf ("hquad (%d,%d) = ", nc, nz);
 189               ker_hquad [IndElt (nc,nz)]->dumpPlus();
 190               }
 191            }
 192
 193        if (nz>0)
 194           {
 195                                  //   Edges verticaux + cloisons interieures
 196           Edge* pilier = ker_vedge [nz-1] = newEdge (pcenter, center);
 197
 198           for (int nc=0 ; nc<nbrayons ; nc++)
 199               {
 200               Edge* e1 = ker_hedge [IndRedge (nc, nz)];
 201               Edge* e2 = getEdgeK (nx0, 2*nc,  nz-1);
 202               Edge* e3 = ker_hedge [IndRedge (nc, nz-1)];
 203               ker_vquad [IndVquad (nc,nz)] = newQuad (e1, e2, e3, pilier);
 204               if (debug())
 205                  {
 206                  printf ("vquad (%d,%d) = ", nc, nz);
 207                  ker_vquad [IndVquad (nc,nz)]->dumpPlus();
 208                  }
 209               }
 210                                  //   Hexas
 211           for (int nc=0 ; nc < nbsecteurs ; nc++)
 212               {
 213               int nc1  = nc + 1;
 214               if (cyl_closed)
 215                   nc1  = nc1 MODULO nbsecteurs;
 216               Quad* qa = ker_hquad [IndElt (nc, nz-1)];
 217               Quad* qb = ker_hquad [IndElt (nc, nz)];
 218
 219               Quad* qc = ker_vquad [IndVquad (nc, nz)];
 220               Quad* qd = getQuadJK (nx0, 2*nc+1, nz-1);
 221
 222               Quad* qe = getQuadJK (nx0, 2*nc,   nz-1);
 223               Quad* qf = ker_vquad [IndVquad (nc1,  nz)];
 224
 225               if (debug())
 226                  {
 227                  printf (" --------------- Hexa : nc=%d, nz=%d\n", nc, nz);
 228                  HexDump (qa);
 229                  HexDump (qb);
 230                  HexDump (qc);
 231                  HexDump (qd);
 232                  HexDump (qe);
 233                  HexDump (qf);
 234                  }
 235               Hexa* cell = newHexa (qa, qb, qc, qd, qe, qf);
 236               tab_hexa.push_back (cell);
 237               ker_hexa [IndElt (nc,  nz-1)] = cell;
 238               }
 239           }
 240        pcenter = center;
 241        }
 242 }
 243 // ====================================================== fillCenter4
 244 // === Remplissage radial
 245 void Elements::fillCenter4 ()
 246 {
 247    int   nx0  = 0;
 248    Quad* sol = NULL;
 249
 250    for (int nz=0 ; nz<size_vz ; nz++)
 251        {
 252                                  //   Quad horizontal
 253
 254        Quad* plafond = newQuad (getEdgeJ (nx0, 0,  nz), getEdgeJ (nx0, 1, nz),
 255                                 getEdgeJ (nx0, 2,  nz), getEdgeJ (nx0, 3, nz));
 256
 257        if (nz>0)
 258           {
 259           Quad* qc = getQuadJK (nx0, 0, nz-1);
 260           Quad* qd = getQuadJK (nx0, 2, nz-1);
 261           Quad* qe = getQuadJK (nx0, 1, nz-1);
 262           Quad* qf = getQuadJK (nx0, 3, nz-1);
 263
 264           if (debug())
 265              {
 266              printf (" --------------- Hexa grille4 : nz=%d\n", nz);
 267              HexDump (plafond);
 268              HexDump (sol);
 269              HexDump (qc);
 270              HexDump (qd);
 271              HexDump (qe);
 272              HexDump (qf);
 273              }
 274           Hexa* cell = newHexa (plafond, sol, qc, qd, qe, qf);
 275           tab_hexa.push_back (cell);
 276           }
 277        sol = plafond;
 278        }
 279 }
 280 // ====================================================== fillCenter6
 281 void Elements::fillCenter6 ()
 282 {
 283    int nx0 = 0;
 284    int nydemi = size_hy / 2;
 285
 286    Edge* s_barre = NULL;
 287    Quad* sr_quad = NULL;
 288    Quad* sl_quad = NULL;
 289
 290    for (int nz=0 ; nz<size_vz ; nz++)
 291        {
 292                                  //   Edges horizontal radial
 293        Edge* p_barre = newEdge (getVertexIJK (nx0, 0,      nz),
 294                                 getVertexIJK (nx0, nydemi, nz));
 295                                  //   Quads horizontaux
 296        Edge* e0 = getEdgeJ (nx0, 0,  nz);
 297        Edge* e1 = getEdgeJ (nx0, 1,  nz);
 298        Edge* e2 = getEdgeJ (nx0, 2,  nz);
 299        Edge* e3 = getEdgeJ (nx0, 3,  nz);
 300        Edge* e4 = getEdgeJ (nx0, 4,  nz);
 301        Edge* e5 = getEdgeJ (nx0, 5,  nz);
 302
 303        Quad* pl_quad = newQuad (e0, e1, e2, p_barre);
 304        Quad* pr_quad = newQuad (e3, e4, e5, p_barre);
 305
 306        if (nz>0)
 307           {
 308                                  //   Cloison interieure
 309           Quad* cloison = newQuad (p_barre, getEdgeK (nx0, 0,       nz-1),
 310                                    s_barre, getEdgeK (nx0, nydemi,  nz-1));
 311                                  //   Hexas
 312           Quad* q0 = getQuadJK (nx0, 0, nz-1);
 313           Quad* q1 = getQuadJK (nx0, 1, nz-1);
 314           Quad* q2 = getQuadJK (nx0, 2, nz-1);
 315           Quad* q3 = getQuadJK (nx0, 3, nz-1);
 316           Quad* q4 = getQuadJK (nx0, 4, nz-1);
 317           Quad* q5 = getQuadJK (nx0, 5, nz-1);
 318
 319           Hexa* left  = newHexa (sl_quad, pl_quad, q0, q2,  q1, cloison);
 320           Hexa* right = newHexa (sr_quad, pr_quad, q3, q5,  q4, cloison);
 321           tab_hexa.push_back (left);
 322           tab_hexa.push_back (right);
 323           }
 324        s_barre = p_barre;
 325        sr_quad = pr_quad;
 326        sl_quad = pl_quad;
 327        }
 328 }
 329 // ====================================================== fillCenterOdd
 330 #undef  IndElt
 331 #define IndElt(nc,nz) (nbsecteurs*(nz) + nc)
 332 void Elements::fillCenterOdd ()
 333 {
 334    int nx0 = 0;
 335    int nbsecteurs = size_hy / 2;
 336
 337    std::vector <Edge*> ker_hedge (nbsecteurs*size_vz);
 338    std::vector <Quad*> ker_hquad (nbsecteurs*size_vz);
 339    std::vector <Quad*> ker_vquad (nbsecteurs*size_vz);
 340
 341    for (int nz=0 ; nz<size_vz ; nz++)
 342        {
 343                                  //   Edges horizontaux radiaux
 344        for (int nc=0 ; nc<nbsecteurs ; nc++)
 345            {
 346            int nc1 = size_hy - nc;
 347            Edge* edge = newEdge (getVertexIJK (nx0, nc,  nz),
 348                                  getVertexIJK (nx0, nc1, nz));
 349            ker_hedge [IndElt(nc,nz)] = edge;
 350            }
 351                                  //   Quads horizontaux
 352        for (int nc=0; nc<nbsecteurs ; nc++)
 353            {
 354            Edge* e1 = ker_hedge [IndElt (nc, nz)];
 355            Edge* e2 = getEdgeJ (nx0, nc,           nz);
 356            Edge* e4 = getEdgeJ (nx0, size_hy-nc-1, nz);
 357
 358            Edge* e3 = nc<nbsecteurs-1 ? ker_hedge [IndElt (nc+1, nz)]
 359                                       : getEdgeJ (nx0, nbsecteurs, nz);
 360
 361            ker_hquad [IndElt (nc,nz)] = newQuad (e1, e2, e3, e4);
 362            if (debug())
 363               {
 364               printf ("hquad (%d,%d) = ", nc, nz);
 365               ker_hquad [IndElt (nc,nz)]->dumpPlus();
 366               }
 367            }
 368
 369        if (nz>0)
 370           {
 371                                  //   Edges verticaux + cloisons interieures
 372           for (int nc=0 ; nc<nbsecteurs ; nc++)
 373               {
 374               int nc1 = size_hy - nc;
 375               Edge* e1 = ker_hedge [IndElt (nc, nz)];
 376               Edge* e3 = ker_hedge [IndElt (nc, nz-1)];
 377               Edge* e2 = getEdgeK (nx0, nc,   nz-1);
 378               Edge* e4 = getEdgeK (nx0, nc1,  nz-1);
 379               ker_vquad [IndElt (nc,nz)] = newQuad (e1, e2, e3, e4);
 380               if (debug())
 381                  {
 382                  printf ("vquad (%d,%d) = ", nc, nz);
 383                  ker_vquad [IndElt (nc,nz)]->dumpPlus();
 384                  }
 385               }
 386                                  //   Hexas
 387           for (int nc=0 ; nc < nbsecteurs ; nc++)
 388               {
 389               int nc1 = size_hy - nc-1;
 390               Quad* qa = ker_hquad [IndElt (nc, nz-1)];
 391               Quad* qb = ker_hquad [IndElt (nc, nz)];
 392
 393               Quad* qc = getQuadJK (nx0, nc,  nz-1);
 394               Quad* qd = getQuadJK (nx0, nc1, nz-1);
 395
 396               Quad* qe = ker_vquad [IndElt (nc, nz)];
 397               Quad* qf = nc<nbsecteurs-1 ? ker_vquad [IndElt (nc+1,  nz)]
 398                                          : getQuadJK (nx0, nbsecteurs, nz-1);
 399               if (debug())
 400                  {
 401                  printf (" --------------- Hexa : nc=%d, nz=%d\n", nc, nz);
 402                  HexDump (qa);
 403                  HexDump (qb);
 404                  HexDump (qc);
 405                  HexDump (qd);
 406                  HexDump (qe);
 407                  HexDump (qf);
 408                  }
 409               Hexa* cell = newHexa (qa, qb, qc, qd, qe, qf);
 410               tab_hexa.push_back (cell);
 411               }
 412           }
 413        }
 414 }
 415 // --------------------------------------------------------------------------
 416 // ----------------------------------------- Evols Hexa 3
 417 // --------------------------------------------------------------------------
 418
 419 // ====================================================== makeBasicCylinder
 420 // ==== Version avec vecteurs
 421 int Elements::makeBasicCylinder (RealVector& tdr, RealVector& tda,
 422                                  RealVector& tdh, bool fill)
 423 {
 424    int na    = tda.size()-1;
 425    cyl_dispo = CYL_NOFILL;
 426    if (fill && na > 3)
 427       {
 428       if (cyl_closed)
 429          {
 430          if (na==4)
 431             cyl_dispo = CYL_CL4;
 432          else if (na==6)
 433             cyl_dispo = CYL_CL6;
 434          else if (na MODULO 2 == 0)
 435             cyl_dispo = CYL_CLOSED;
 436          }
 437       else if ((na MODULO 2)==0)
 438          cyl_dispo = CYL_PEER;
 439       else
 440          cyl_dispo = CYL_ODD;
 441       }
 442
 443    cyl_fill = cyl_dispo != CYL_NOFILL;
 444
 445    double alpha  = 0;
 446    int    nb_secteurs = cyl_closed ? size_vy-1 : size_vy;
 447
 448    for (int ny=0 ; ny<nb_secteurs ; ny++)
 449        {
 450        alpha  = tda[ny];
 451        double theta     = M_PI*alpha/180;
 452        double cos_theta = cos (theta);
 453        double sin_theta = sin (theta);
 454
 455        for (int nx=0 ; nx<size_vx ; nx++)
 456            {
 457            double rayon = tdr [nx];
 458            double px = rayon*cos_theta;
 459            double py = rayon*sin_theta;
 460
 461            for (int nz=0 ; nz<size_vz ; nz++)
 462                {
 463                double pz = tdh [nz];
 464                Vertex* node = el_root->addVertex (px, py, pz);
 465                setVertex (node, nx, ny, nz);
 466                }
 467            }
 468        }
 469
 470    if (cyl_closed)
 471       {
 472       for (int nx=0 ; nx<size_vx ; nx++)
 473           for (int nz=0 ; nz<size_vz ; nz++)
 474               {
 475               Vertex* node = getVertexIJK (nx, 0, nz);
 476               setVertex (node, nx, size_vy-1, nz);
 477               }
 478       }
 479                       // Les vertex centraux
 480    if (cyl_fill)
 481       {
 482       ker_vertex = nbr_vertex;
 483       for (int nz=0 ; nz<size_vz ; nz++)
 484           {
 485           double pz = tdh [nz];
 486           Vertex* node = el_root->addVertex (0, 0, pz);
 487           tab_vertex.push_back (node);
 488           nbr_vertex ++;
 489           }
 490       }
 491
 492    return HOK;
 493 }
 494 END_NAMESPACE_HEXA
 495