src/HEXABLOCK/HexElements_bis.cxx

   1
   2 // C++ : Table d'hexaedres
   3
   4 #include "HexElements.hxx"
   5
   6 #include "HexDocument.hxx"
   7 #include "HexVector.hxx"
   8 #include "HexVertex.hxx"
   9 #include "HexHexa.hxx"
  10 #include "HexEdge.hxx"
  11
  12 #include "HexGlobale.hxx"
  13 #include "HexCylinder.hxx"
  14
  15 #include <map>
  16
  17 static bool db=false;
  18
  19 BEGIN_NAMESPACE_HEXA
  20
  21 // ====================================================== getHexaIJK
  22 Hexa* Elements::getHexaIJK (int nx, int ny, int nz)
  23 {
  24    if (nx<0 || nx>=size_hx ||  ny<0 || ny>=size_hy || nz<0 || nz>=size_hz)
  25       return NULL;
  26    else if (grid_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
  27       return NULL;
  28
  29    int nro = nx + size_hx*ny + size_hx*size_hy*nz;
  30
  31    return tab_hexa [nro];
  32 }
  33 // ====================================================== getQuadIJ
  34 Quad* Elements::getQuadIJ (int nx, int ny, int nz)
  35 {
  36    if (nx<0 || nx>=size_qx ||  ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
  37       return NULL;
  38    else if (grid_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
  39       return NULL;
  40
  41    int nro = nx + size_qx*ny + size_qx*size_qy*nz
  42                 + size_qx*size_qy*size_qz*dir_z;
  43    return tab_quad [nro];
  44 }
  45 // ====================================================== getQuadJK
  46 Quad* Elements::getQuadJK (int nx, int ny, int nz)
  47 {
  48    if (nx<0 || nx>=size_qx ||  ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
  49       return NULL;
  50    else if (grid_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
  51       return NULL;
  52
  53    int nro = nx + size_qx*ny + size_qx*size_qy*nz; // + dir_x*...
  54
  55    return tab_quad [nro];
  56 }
  57 // ====================================================== getQuadIK
  58 Quad* Elements::getQuadIK (int nx, int ny, int nz)
  59 {
  60    if (nx<0 || nx>=size_qx ||  ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
  61       return NULL;
  62    else if (grid_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
  63       return NULL;
  64
  65    int nro = nx + size_qx*ny + size_qx*size_qy*nz + size_qx*size_qy*size_qz;
  66
  67    return tab_quad [nro];
  68 }
  69 // ====================================================== getEdgeI
  70 Edge* Elements::getEdgeI (int nx, int ny, int nz)
  71 {
  72    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
  73       return NULL;
  74    else if (grid_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
  75       return NULL;
  76
  77    int nro = nx + size_ex*ny + size_ex*size_ey*nz;
  78
  79    return tab_edge [nro];
  80 }
  81 // ====================================================== getEdgeJ
  82 Edge* Elements::getEdgeJ (int nx, int ny, int nz)
  83 {
  84    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
  85       return NULL;
  86    else if (grid_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
  87       return NULL;
  88
  89    int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez;
  90
  91    return tab_edge [nro];
  92 }
  93 // ====================================================== getEdgeK
  94 Edge* Elements::getEdgeK (int nx, int ny, int nz)
  95 {
  96    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
  97       return NULL;
  98    else if (grid_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
  99       return NULL;
 100
 101    int nro = nx + size_ex*ny + size_ex*size_ey*nz
 102                 + size_ex*size_ey*size_ez*dir_z;
 103    return tab_edge [nro];
 104 }
 105 // ====================================================== getVertexIJK
 106 Vertex* Elements::getVertexIJK (int nx, int ny, int nz)
 107 {
 108    if (nx<0 || nx>=size_vx ||  ny<0 || ny>=size_vy || nz<0 || nz>=size_vz)
 109       return NULL;
 110    else if (grid_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
 111       return NULL;
 112
 113    int nro = nx + size_vx*ny + size_vx*size_vy*nz;
 114
 115    return tab_vertex [nro];
 116 }
 117 // ====================================================== setVertex
 118 void Elements::setVertex (Vertex* elt, int nx, int ny, int nz)
 119 {
 120    if (   nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy
 121        || nz < 0 || nz >= size_vz) return;
 122
 123    int nro = nx + size_vx*ny + size_vx*size_vy*nz;
 124    tab_vertex [nro] = elt;
 125 }
 126 // ====================================================== setVertex (2)
 127 void Elements::setVertex (int nx, int ny, int nz, double px, double py,
 128                                                   double pz)
 129 {
 130    if (   nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy
 131        || nz < 0 || nz >= size_vz) return;
 132
 133    Vertex*    node = el_root->addVertex (px, py, pz);
 134    setVertex (node, nx, ny, nz);
 135 }
 136 // ====================================================== setEdge
 137 void Elements::setEdge (Edge* elt, EnumCoord dir, int nx, int ny, int nz)
 138 {
 139    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez
 140             || dir < dir_x || dir > dir_z )
 141       return;
 142
 143    int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez*dir;
 144    tab_edge [nro] = elt;
 145 }
 146 // ====================================================== setQuad
 147 void Elements::setQuad (Quad* elt, EnumCoord dir, int nx, int ny, int nz)
 148 {
 149    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez
 150             || dir < dir_x || dir > dir_z )
 151       return;
 152
 153    int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez*dir;
 154    tab_quad [nro] = elt;
 155 }
 156 // ====================================================== setHexa
 157 void Elements::setHexa (Hexa* elt, int nx, int ny, int nz)
 158 {
 159    if (   nx < 0 || nx >= size_hx || ny < 0 || ny >= size_hy
 160        || nz < 0 || nz >= size_hz) return;
 161
 162    int nro = nx + size_hx*ny + size_hx*size_hy*nz;
 163    tab_hexa [nro] = elt;
 164 }
 165 // ====================================================== remove
 166 void Elements::remove ()
 167 {
 168    int nbre=tab_hexa.size ();
 169    nbre = nbr_hexas;
 170    for (int nh=0 ; nh<nbre ; nh++)
 171        if (tab_hexa[nh] != NULL)
 172            tab_hexa[nh]->remove();
 173 }
 174 // ====================================================== prismQuads
 175 int Elements::prismQuads (Quads& tstart, Vector* dir, int hauteur)
 176 {
 177    el_root->markAll (NO_USED);
 178    int nbcells   = tstart.size ();
 179    nbr_vertex    = 0;
 180    nbr_edges     = 0;
 181
 182    nbr_hexas = nbcells*hauteur;
 183
 184    tab_hexa.resize (nbr_hexas);
 185    tab_quad.clear ();          // verticaux
 186    tab_edge.clear ();
 187    tab_pilier.clear ();
 188    tab_vertex.clear ();
 189
 190    for (int nro=0 ; nro<nbcells ; nro++)
 191        {
 192        pushHexas (nro, tstart[nro], dir->getDx(), dir->getDy(), dir->getDz(),
 193                   hauteur);
 194        }
 195    return HOK;
 196 }
 197 // ====================================================== pushHexas
 198 int  Elements::pushHexas (int nro, Quad* qbase, double dx, double dy,
 199                           double dz, int hauteur, Quad* toit, int tcorr[])
 200 {
 201    int ind_node [QUAD4];
 202
 203            // ----------------------------- Vertex + aretes verticales
 204    for (int ns=0 ; ns<QUAD4 ; ns++)
 205        {
 206        Vertex* vbase = qbase ->getVertex (ns);
 207        int     indx  = vbase->getMark ();
 208        if (indx<0)
 209           {
 210           indx = nbr_vertex++;
 211           vbase->setMark (indx);
 212           Vertex* nd0 = vbase;
 213           Vertex* nd1 = NULL;
 214           for (int nh=0 ; nh<hauteur ; nh++)
 215               {
 216               nd1 = el_root->addVertex (nd0->getX() + dx, nd0->getY() + dy,
 217                                         nd0->getZ() + dz);
 218               tab_vertex.push_back (nd1);
 219               tab_pilier.push_back (newEdge (nd0, nd1));
 220               nd0 = nd1;
 221               }
 222           }
 223        ind_node [ns] = indx;
 224        }
 225            // ----------------------------- Aretes horizontales
 226            // ----------------------------- + face verticales
 227    int ind_poutre [QUAD4];
 228    for (int ns=0 ; ns<QUAD4 ; ns++)
 229        {
 230        Edge* ebase = qbase ->getEdge (ns);
 231        int   indx  = ebase->getMark ();
 232        if (indx<0)
 233           {
 234           indx = nbr_edges ++;
 235           ebase->setMark (indx);
 236           int   nd1 = ind_node [ns];
 237           int   nd2 = ind_node [(ns+1) MODULO QUAD4];
 238           Edge* ed0 = ebase;
 239           Edge *ed1, *ed2, *ed3;
 240           for (int nh=0 ; nh<hauteur ; nh++)
 241               {
 242               ed2 = ed0;
 243               ed0 = newEdge (tab_vertex [nd1*hauteur + nh],
 244                              tab_vertex [nd2*hauteur + nh]);
 245               ed1 = tab_pilier [nd1*hauteur + nh];
 246               ed3 = tab_pilier [nd2*hauteur + nh];
 247
 248               Quad* mur = newQuad (ed0, ed1, ed2, ed3);
 249               tab_edge.push_back (ed0);
 250               tab_quad.push_back (mur);
 251               }
 252           }
 253        ind_poutre [ns] = indx;
 254        }
 255            // ----------------------------- Faces horizontales
 256            // ----------------------------- + Hexas
 257    Quad* qa = qbase;
 258    Quad *qb, *qc, *qd, *qe, *qf;
 259    int nv0 = hauteur*ind_poutre [0];
 260    int nv1 = hauteur*ind_poutre [1];
 261    int nv2 = hauteur*ind_poutre [2];
 262    int nv3 = hauteur*ind_poutre [3];
 263    for (int nh=0 ; nh<hauteur ; nh++)
 264        {
 265        qb = newQuad (tab_edge [nh+nv0], tab_edge [nh+nv1],
 266                      tab_edge [nh+nv2], tab_edge [nh+nv3]);
 267        qc = tab_quad [nh + nv0];
 268        qd = tab_quad [nh + nv2];
 269        qe = tab_quad [nh + nv1];
 270        qf = tab_quad [nh + nv3];
 271
 272 // *** tab_hexa [nh*hauteur + nro] = newHexa (qa, qb, qc, qd, qe, qf); Abu
 273        tab_hexa [nro*hauteur + nh] = newHexa (qa, qb, qc, qd, qe, qf);
 274        qa = qb;
 275        }
 276    return HOK;
 277 }
 278 // ====================================================== cutHexas
 279 int Elements::cutHexas0 (const Edges& t_edges, int nbcuts)
 280 {
 281                                        // 1) marquage des hexas
 282    el_root->markAll (NO_USED);
 283                                        // 2) Memo noeuds
 284    vector <Quad*> q_amont;
 285    vector <Quad*> q_aval;
 286
 287    int nbnodes = t_edges.size();
 288    vector <Vertex*> v_amont (nbnodes);
 289    vector <Vertex*> v_aval  (nbnodes);
 290
 291    std::map <Vertex*, Vertex*> vertex_aval;
 292
 293    for (int nro=0; nro<nbnodes ; nro++)
 294        {
 295        Edge*   arete = t_edges [nro];
 296        Vertex* amont = arete->getAmont ();
 297        v_amont [nro] = amont;
 298        vertex_aval [amont] = arete->getAval  ();
 299        int nbcells   = arete->getNbrParents ();
 300
 301        for (int nq=0 ; nq<nbcells ; nq++)
 302            {
 303            Quad* quad = arete->getParent (nq);
 304            if (quad->getMark () != IS_USED)
 305               {
 306               quad->setMark (IS_USED);
 307               int nbcubes = quad->getNbrParents ();
 308               for (int nh=0 ; nh<nbcubes ; nh++)
 309                   {
 310                   Hexa* hexa = quad->getParent (nh);
 311                   if (hexa->getMark () != IS_USED)
 312                      {
 313                      hexa->setMark (IS_USED);
 314                      int namont = hexa->getBase (v_amont[nro], arete);
 315                      int naval  = glob->getOpposedQuad (namont);
 316                      q_amont.push_back (hexa->getQuad  (namont));
 317                      q_aval .push_back (hexa->getQuad  (naval ));
 318
 319                      if (db)
 320                         {
 321                         printf (" Quad = ");
 322                         hexa->printName (", ");
 323                         printf (" Faces = (");
 324                         hexa->getQuad (namont)->printName (", ");
 325                         hexa->getQuad (naval )->printName (")\n");
 326                         }
 327                      }
 328                   }
 329               }
 330            }
 331        }
 332                    // ------------------- Dimensionnement
 333    int hauteur = nbcuts + 1;
 334    int nbcells = q_amont.size ();
 335    nbr_hexas   = nbcells * hauteur;
 336
 337    tab_hexa.resize   (nbr_hexas);
 338    tab_quad.clear ();          // verticaux
 339    tab_edge.clear ();
 340    tab_pilier.clear ();
 341    tab_vertex.clear ();
 342
 343    for (int ned=0; ned<nbnodes ; ned++)
 344        t_edges [ned]->remove ();
 345
 346
 347    for (int nro=0; nro<nbcells ; nro++)
 348        {
 349        int tcorr [QUAD4];
 350        Quad* sol  = q_amont[nro];
 351        Quad* toit = q_aval [nro];
 352
 353        for (int nv=0; nv<QUAD4 ; nv++)
 354            {
 355            Vertex* dest = vertex_aval [sol->getVertex(nv)];
 356            tcorr [nv] = toit->indexVertex (dest);
 357            }
 358
 359        pushHexas (nro, sol, 0,0,0, hauteur, toit, tcorr);
 360        }
 361    return HOK;
 362 }
 363 // ====================================================== makeCylinder
 364 int Elements::makeCylinder (Cylinder* cyl, Vector* vx, int nr, int na, int nl)
 365 {
 366    Vertex* orig = cyl->getBase ();
 367    Vector* dir  = cyl->getDirection ();
 368    double  ray  = cyl->getRadius ();
 369    double  haut = cyl->getHeight ();
 370
 371    resize (GR_CYLINDRIC, nr, na, nl);
 372    cyl_closed = true;
 373    makeCylindricalNodes (orig, vx, dir, ray, 360, haut, nr, na, nl, true);
 374    fillGrid ();
 375    return HOK;
 376 }
 377 // ====================================================== makePipe
 378 int Elements::makePipe (Cylinder* cyl, Vector* vx, int nr, int na, int nl)
 379 {
 380    Vertex* orig = cyl->getBase ();
 381    Vector* dir  = cyl->getDirection ();
 382    double  ray  = cyl->getRadius ();
 383    double  haut = cyl->getHeight ();
 384
 385    resize (GR_CYLINDRIC, nr, na, nl);
 386    cyl_closed = true;
 387    makeCylindricalNodes (orig, vx, dir, ray, 360, haut, nr, na, nl, false);
 388    fillGrid ();
 389    return HOK;
 390 }
 391
 392 END_NAMESPACE_HEXA