src/HEXABLOCK/HexElements_v6.cxx

   1
   2 // C++ : Grilles hexa 6
   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 #include "HexGlobale.hxx"
  25 #include "HexDocument.hxx"
  26 #include "HexVector.hxx"
  27 #include "HexVertex.hxx"
  28 #include "HexHexa.hxx"
  29 #include "HexEdge.hxx"
  30 #include "HexPropagation.hxx"
  31
  32 #include <cmath>
  33 #include <map>
  34
  35 BEGIN_NAMESPACE_HEXA
  36
  37 static bool db=false;
  38 double calcul_phimax (double radhole, double radext, bool sphere);
  39
  40 // ====================================================== makeCartesian
  41 int Elements::makeCartesian (Vertex* orig, double* ux, double* uy, double* uz,
  42                              RealVector& tx, RealVector& ty, RealVector& tz)
  43 {
  44    resize (GR_CARTESIAN, tx.size(), ty.size(), tz.size());
  45
  46    double dx, dy, dz;
  47    double kx, ky, kz;
  48    double px0 = orig->getX ();
  49    double py0 = orig->getY ();
  50    double pz0 = orig->getZ ();
  51
  52    setVertex (orig, 0, 0, 0);
  53    for (int nz=0 ; nz<size_vz ; nz++)
  54        {
  55        kz = nz==0 ? 0 : tz [nz-1];
  56        for (int ny=0 ; ny<size_vy ; ny++)
  57            {
  58            ky = ny==0 ? 0 : ty [ny-1];
  59            for (int nx=0 ; nx<size_vx ; nx++)
  60                {
  61                kx = nx==0 ? 0 : tx [nx-1];
  62                dx = px0 + kx*ux[dir_x] + ky*uy[dir_x] + kz*uz[dir_x];
  63                dy = py0 + kx*ux[dir_y] + ky*uy[dir_y] + kz*uz[dir_y];
  64                dz = pz0 + kx*ux[dir_z] + ky*uy[dir_z] + kz*uz[dir_z];
  65
  66                Vertex* node = orig;
  67                if (nx+ny+nz !=0)
  68                   node = el_root->addVertex (dx, dy, dz);
  69                setVertex (node, nx, ny, nz);
  70                }
  71            }
  72        }
  73    fillGrid ();
  74    return HOK;
  75 }
  76 // ====================================================== makeCylinder
  77 int Elements::makeCylinder (Vertex* vorig, double* base, double* haut,
  78                          RealVector& tray, RealVector& tang, RealVector& thaut,
  79                             bool fill)
  80 {
  81    Real3 orig = { 0, 0, 0 };
  82    if (vorig!=NULL)
  83        vorig->getPoint (orig);
  84
  85    int nr = tray .size() - 1;
  86    int na = tang .size() - 1;
  87    int nl = thaut.size() - 1;
  88    double angle = tang[na];
  89
  90    resize (GR_CYLINDRIC, nr, na, nl);
  91    cyl_closed = angle >= 359.9;
  92
  93    val_absolues = true;
  94    int ier = makeBasicCylinder (tray, tang, thaut, fill);
  95    if (ier!=HOK)
  96        return ier;
  97
  98    transfoVertices  (orig,  base, haut);
  99
 100    fillGrid ();
 101    assoCylinders (orig, haut, angle, tang);
 102    return HOK;
 103 }
 104 // ====================================================== transfoVertices
 105 void Elements::transfoVertices (double* omega, double* base, double* haut)
 106 {
 107    Matrix matrice;
 108    Real3  iprim, jprim;
 109
 110    prod_vectoriel (haut,  base, jprim);
 111    prod_vectoriel (jprim, haut, iprim);
 112
 113    matrice.define (omega, iprim, jprim, haut);
 114    int nbre = tab_vertex.size ();
 115    for (int nro=0 ; nro<nbre ; nro++)
 116        {
 117        if (tab_vertex[nro] != NULL)
 118            tab_vertex[nro]->setMark (NO_USED);
 119        }
 120
 121    for (int nro=0 ; nro<nbre ; nro++)
 122        {
 123        Vertex* node =  tab_vertex[nro];
 124        if (node != NULL && node->getMark() == NO_USED)
 125           {
 126           matrice.perform (node);
 127           node->setMark (IS_USED);
 128           }
 129        }
 130 }
 131 // ====================================================== makeSpherical
 132 int Elements::makeSpherical (double* center, double* base, double* haut,
 133                              RealVector& trayon, int critere)
 134 {
 135    Real3  iprim, jprim;
 136    prod_vectoriel (haut,  base, jprim);
 137    prod_vectoriel (jprim, haut, iprim);
 138
 139    int nbre = trayon.size();
 140    resize (GR_SPHERIC, nbre);
 141
 142    Vertex* i_node [HV_MAXI];    // Les noeuds de l'hexa englobant
 143    Edge*   i_edge [HE_MAXI];    // Les noeuds de l'hexa englobant
 144    Quad*   i_quad [HQ_MAXI];    // Les noeuds de l'hexa englobant
 145
 146    double rayon = trayon [0];
 147    for (int nro=0 ; nro<HV_MAXI; nro++)
 148        {
 149        double dx = glob->CoordVertex (nro, dir_x) * rayon;
 150        double dy = glob->CoordVertex (nro, dir_y) * rayon;
 151        double dz = glob->CoordVertex (nro, dir_z) * rayon;
 152
 153        double px = center[0] + dx*iprim[0] + dy*jprim[0] + dz*haut[0] ;
 154        double py = center[1] + dx*iprim[1] + dy*jprim[1] + dz*haut[1] ;
 155        double pz = center[2] + dx*iprim[2] + dy*jprim[2] + dz*haut[2] ;
 156
 157        i_node [nro] = el_root->addVertex (px, py, pz);
 158        }
 159
 160    for (int nro=0 ; nro<HE_MAXI; nro++)
 161        {
 162        int v1 = glob->EdgeVertex (nro, V_AMONT);
 163        int v2 = glob->EdgeVertex (nro, V_AVAL);
 164        i_edge[nro] = newEdge (i_node[v1], i_node[v2]);
 165
 166        if (db)
 167           {
 168           char nm0[8], nm1 [8], nm2 [8];
 169           printf (" %2d : %s = %s = [%s, %s] = [%d,%d] = [%s,%s]\n", nro,
 170                  glob->namofHexaEdge(nro), i_edge[nro]->getName(nm0),
 171                  glob->namofHexaVertex(v1), glob->namofHexaVertex(v2), v1, v2,
 172                  i_node[v1]->getName(nm1), i_node[v2]->getName(nm2));
 173           }
 174        }
 175
 176    for (int nro=0 ; nro<HQ_MAXI; nro++)
 177        i_quad[nro] = newQuad (i_edge[glob->QuadEdge (nro, E_A)],
 178                               i_edge[glob->QuadEdge (nro, E_B)],
 179                               i_edge[glob->QuadEdge (nro, E_C)],
 180                               i_edge[glob->QuadEdge (nro, E_D)]);
 181
 182    tab_hexa.push_back (newHexa (i_quad[Q_A], i_quad[Q_B], i_quad[Q_C],
 183                                 i_quad[Q_D], i_quad[Q_E], i_quad[Q_F]));
 184
 185    for (int niv=1; niv<gr_rayon ; niv++)
 186        {
 187        double rayon0  = rayon;
 188        rayon  = trayon [niv];
 189        addStrate (i_quad, i_edge, i_node, center, rayon/rayon0);
 190        }
 191
 192    nbr_hexas = tab_hexa.size();
 193    assoSphere (center, i_edge, i_quad);
 194    return HOK;
 195 }
 196 // ====================================================== addStrate
 197 int Elements::addStrate (Quad* i_quad[], Edge* i_edge[], Vertex* i_node[],
 198                          Real* center,  double lambda)
 199 {
 200    Vertex* e_node [HV_MAXI];    // Les noeuds de l'hexa englobant
 201    Edge*   e_edge [HE_MAXI];    // Les noeuds de l'hexa englobant
 202    Quad*   e_quad [HQ_MAXI];    // Les noeuds de l'hexa englobant
 203
 204    Edge*   d_edge [HV_MAXI];    // Les aretes diagonales (1 par sommet)
 205    Quad*   d_quad [HE_MAXI];    // Les faces  diagonales (1 par arete)
 206
 207                                           // Les sommets
 208                                           //  + les aretes diagonales
 209    double px0 = center [dir_x];
 210    double py0 = center [dir_y];
 211    double pz0 = center [dir_z];
 212
 213    for (int nv=0 ; nv<HV_MAXI ; nv++)
 214        {
 215        e_node[nv] = el_root->addVertex (px0+lambda*(i_node[nv]->getX()-px0),
 216                                         py0+lambda*(i_node[nv]->getY()-py0),
 217                                         pz0+lambda*(i_node[nv]->getZ()-pz0));
 218
 219        d_edge[nv] = newEdge (i_node[nv], e_node[nv]);
 220        }
 221                                           // Les aretes exterieures
 222                                           //  + les faces diagonales
 223    for (int nro=0 ; nro<HE_MAXI ; nro++)
 224        {
 225        int nv0  = glob->EdgeVertex (nro, V_AMONT);
 226        int nv1  = glob->EdgeVertex (nro, V_AVAL );
 227        e_edge[nro] = newEdge (e_node [nv0], e_node [nv1]);
 228        d_quad[nro] = newQuad (i_edge [nro], d_edge [nv0],
 229                               e_edge [nro], d_edge [nv1]);
 230        }
 231                                           // Les faces exterieures
 232                                           //  + les hexas
 233    Hexa* strate = NULL;
 234    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 235        {
 236        int ne0 = glob->QuadEdge (nro, E_A);
 237        int ne1 = glob->QuadEdge (nro, E_B);
 238        int ne2 = glob->QuadEdge (nro, E_C);
 239        int ne3 = glob->QuadEdge (nro, E_D);
 240
 241        e_quad[nro] = newQuad (e_edge[ne0], e_edge[ne1],
 242                               e_edge[ne2], e_edge[ne3]);
 243        strate  = newHexa (i_quad[nro], e_quad[nro], d_quad[ne0],
 244                           d_quad[ne2], d_quad[ne1], d_quad[ne3]);
 245        tab_hexa.push_back (strate);
 246        }
 247
 248    for (int nv=0 ; nv<HV_MAXI ; nv++) i_node  [nv] = e_node [nv];
 249    for (int ns=0 ; ns<HE_MAXI ; ns++) i_edge  [ns] = e_edge [ns];
 250    for (int nq=0 ; nq<HQ_MAXI ; nq++) i_quad  [nq] = e_quad [nq];
 251    return HOK;
 252 }
 253 // ====================================================== makeRind
 254 int Elements::makeRind (EnumGrid type, double* center, double* vx, double* vz,
 255                         RealVector& tray, RealVector& tang, RealVector& tphi)
 256 {
 257    int na = tang.size() - 1;
 258    int nr = tray.size() - 1;
 259    int nh = tphi.size() - 1;
 260    int rad0 = 0;
 261    int nint = 0;
 262
 263    bool sphere = true;
 264    if (type==GR_PART_RIND || type==GR_RIND)
 265       {
 266       sphere = false;
 267       nr   --;
 268       rad0 = 1;
 269       nint = 1;
 270       }
 271
 272    resize (type, nr, na, nh);
 273
 274    cyl_radhole = tray[0];
 275    cyl_radint  = tray[nint];
 276    cyl_radext  = tray[rad0+nr];
 277    cyl_closed  = type==GR_HEMISPHERIC || type==GR_RIND;
 278    cyl_phimax  = calcul_phimax (cyl_radhole, cyl_radext, sphere);
 279    PutData (cyl_phimax);
 280    int nb_secteurs = cyl_closed ? size_vy-1 : size_vy;
 281
 282    for (int ny=0 ; ny<nb_secteurs ; ny++)
 283        {
 284        double theta = tang[ny];
 285        for (int nx=0 ; nx<size_vx ; nx++)
 286            {
 287            double rayon = tray[rad0+nx];
 288            if (sphere && ny==0 && nx==0)
 289                printf (" rayon = %g, rad_int=%g\n", rayon, cyl_radint);
 290
 291            for (int nz=0 ; nz<size_vz ; nz++)
 292                {
 293                double px, py, pz;
 294                double phi = tphi [nz];
 295                makeRindPoint (rayon, theta, phi, px, py, pz);
 296                Vertex* node = el_root->addVertex (px, py, pz);
 297                setVertex (node, nx, ny, nz);
 298                }
 299            }
 300        }
 301
 302    if (cyl_closed)
 303       {
 304       for (int nx=0 ; nx<size_vx ; nx++)
 305           for (int nz=0 ; nz<size_vz ; nz++)
 306               {
 307               Vertex* node = getVertexIJK (nx, 0, nz);
 308               setVertex (node, nx, size_vy-1, nz);
 309               }
 310       }
 311
 312    transfoVertices (center, vx, vz);
 313    fillGrid ();
 314    double angle = tang[na];
 315    assoCylinders (center, vz, angle, tang);
 316    return HOK;
 317 }
 318 // ====================================================== makeRindPoint
 319 int Elements::makeRindPoint (double ray, double theta, double phi,
 320                              double& px, double& py, double& pz)
 321 {
 322    bool   rind   = (grid_type == GR_RIND || grid_type == GR_PART_RIND);
 323    double rphi   = deg2radians (phi);
 324    rphi = std::min (cyl_phimax, std::max(-cyl_phimax, rphi));
 325    double rtheta = deg2radians (theta);
 326    double sinphi = sin (rphi);
 327    double cosphi = cos (rphi);
 328
 329    double rayon = 0;
 330    if (rind)
 331       {
 332       // rayon = cyl_radint + nr*(cyl_radext-cyl_radint)/size_hx;
 333       // pz    = rayon*sinphi;
 334       // rayon = rayon*cosphi;
 335       pz    = ray*sinphi;
 336       rayon = ray*cosphi;
 337       }
 338    else
 339       {
 340       // rayon = cyl_radhole + nr*(cyl_radext*cosphi - cyl_radhole)/size_hx;
 341       pz    = cyl_radext*sinphi;
 342       rayon = cyl_radhole + (ray-cyl_radhole)*(cyl_radext*cosphi-cyl_radhole)
 343                           / (cyl_radext-cyl_radhole);
 344       }
 345
 346    rayon = std::max (cyl_radhole, rayon);
 347    px = rayon * cos (rtheta);
 348    py = rayon * sin (rtheta);
 349    return HOK;
 350 }
 351 // ====================================================== extrudeQuads
 352 int Elements::extrudeQuads (Quads& tstart, double* axe, RealVector& tabval,
 353                             bool revol, Vertex* center)
 354 {
 355    under_v6    = false;
 356    int nbiter  = tabval.size ();
 357    int nbcells = tstart.size ();
 358
 359    el_root->markAll (NO_USED);
 360    nbr_vertex    = 0;
 361    nbr_edges     = 0;
 362
 363    nbr_hexas   = nbcells*nbiter;
 364
 365    tab_hexa.resize (nbr_hexas, NULL);
 366    tab_quad.clear ();          // verticaux
 367    ker_hquad.clear ();         // Horizontaux
 368    tab_edge.clear ();
 369    tab_pilier.clear ();
 370    tab_vertex.clear ();
 371    copy_vecteur (axe, prism_dir);
 372
 373    revo_lution  = revol;
 374    prism_vec    = NOT revol;
 375    copy_vecteur (axe, revo_axe);
 376    revo_center  = center;
 377    gen_values   = tabval;
 378    cum_values   = gen_values;
 379
 380    for (int nro=1 ; nro<nbiter ; nro++)
 381        gen_values [nro] = cum_values [nro] - cum_values [nro-1];
 382
 383    grid_geom   = NULL;
 384    getShape ();
 385
 386    for (int nro=0 ; nro<nbcells ; nro++)
 387        {
 388        prismHexas (nro, tstart[nro], nbiter);
 389        }
 390
 391    endPrism ();
 392    return HOK;
 393 }
 394 // ============================================================ joinQuads
 395 int Elements::joinQuads (Quads& tstart, Quad* cible, Vertex* va1, Vertex* vb1,
 396                          Vertex* va2, Vertex* vb2, RealVector& tlen)
 397 {
 398    under_v6    = false;
 399    int nbiter  = tlen  .size();
 400    int nbquads = tstart.size();
 401    resize (GR_JOINT, nbquads, nbiter);
 402
 403    el_root->markAll (IS_NONE);
 404    db = on_debug();
 405
 406    gr_hauteur = nbiter;
 407    nbr_orig   = nbquads;
 408    tab_orig   = tstart;
 409    cum_values = tlen;
 410
 411    for (int nq=0 ; nq<nbquads ; nq++)
 412          tab_orig [nq]->setMark (nq);
 413
 414    StrOrient orient (va1, va2, vb1, vb2);
 415    int ier =  this   ->coupler (0, cible, &orient);
 416    if (ier!=HOK)
 417       {
 418       setError ();
 419       return HERR;
 420       }
 421    ier = tstart[0]->coupler (cible, &orient, this);
 422    return ier;
 423 }
 424 // ============================================================ cutEdge
 425 int Elements::cutEdge (Edge* edge, RealVector& tlen)
 426 {
 427    under_v6    = false;
 428    gen_values  = tlen;
 429    cum_values  = gen_values;
 430
 431    Propagation* prop    = el_root->findPropagation (edge);
 432    const Edges& t_edges = prop->getEdges ();
 433
 434    int nbcuts = tlen.size ();
 435    cutHexas (t_edges, nbcuts);
 436
 437    // el_root->majPropagation ();
 438    return HOK;
 439 }
 440 // ============================================================ nearestVertex
 441 Vertex* Elements::nearestVertex (Vertex* other)
 442 {
 443    if (BadElement (other))
 444       return NULL;
 445
 446    Vertex* vsol = NULL;
 447    double  dmin = 1e+77;
 448    int nbre = countVertex ();
 449    for (int nro=0 ; nro<nbre ; nro++)
 450        {
 451        Vertex* vert = getVertex (nro);
 452        double  dist = other->dist2 (vert);
 453        if (dist < dmin)
 454           {
 455           dmin = dist;
 456           vsol = vert;
 457           }
 458        }
 459    return vsol;
 460 }
 461 END_NAMESPACE_HEXA