src/HEXABLOCK/HexElements_piq.cxx

   1
   2 // C++ : Table d'hexaedres (Evol Versions 3)
   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 "HexVector.hxx"
  26 #include "HexVertex.hxx"
  27 #include "HexEdge.hxx"
  28 #include "HexDiagnostics.hxx"
  29 #include "HexDocument.hxx"
  30 #include "HexHexa.hxx"
  31 #include "HexGlobale.hxx"
  32 #include "HexQpattern.hxx"
  33
  34 #include <cmath>
  35 #include <map>
  36
  37 BEGIN_NAMESPACE_HEXA
  38
  39 static bool db = false;
  40
  41 // ====================================================== opposed_edge
  42 Edge* opposed_edge (Edge* edge, Hexas& thexas, Quads& tquads)
  43 {
  44    int nbquads = tquads.size();
  45    for (int nq = 0 ; nq<nbquads ; ++nq)
  46        {
  47        Quad* quad = tquads[nq];
  48        int nro = quad->indexEdge (edge);
  49        if (nro>=0)
  50           {
  51           Quad*  perp  = thexas[nro]->perpendicularQuad (quad, edge);
  52           Edge*  oppos = perp->getOpposEdge (edge);
  53           return oppos;
  54           }
  55        }
  56    return NULL;
  57 }
  58 // ====================================================== replaceHexas
  59 int Elements::replaceHexas (Quads& motif, Quads& cible, Vertex* p1,
  60                             Vertex* c1, Vertex* p2, Vertex* c2)
  61 {
  62     Vertices ed_motif, ed_cible;
  63     checkContour (motif, p1, p2, false, ed_motif);
  64     checkContour (cible, c1, c2, true , ed_cible);
  65     if (el_status != HOK)
  66         return el_status;
  67
  68     int nbedp   = ed_motif.size();
  69     int nbedges = ed_cible.size();
  70     if (nbedp != nbedges)
  71        {
  72        setError (HERR);
  73        Mess << " Contour of pattern ("  << nbedp << " edges) and target ("
  74             << nbedges << " edges) must have the same size ";
  75        return el_status;
  76        }
  77                         // Constitution du pattern
  78     Qpattern pattern;
  79     pattern.setPattern (ed_motif, motif);
  80     pattern.setTarget  (ed_cible, cible);
  81     int nbq, nbed, nbv, prof;
  82     pattern.getSize (nbq, nbed, nbv, prof);
  83
  84     resize (GR_REPLACE, nbq, prof+1, nbv, nbed);
  85
  86     replaceQuad (1, &pattern);
  87     for (int nh=1 ; nh<=prof ; ++nh)
  88         {
  89         pattern.stepDown();
  90         replaceQuad (nh+1, &pattern);
  91         replaceHexa (nh,   &pattern);
  92         }
  93
  94     extrudeQuad (&pattern);
  95     replaceHexa (0, &pattern);
  96     return HOK;
  97 }
  98 // -------------------------------------------------------------------
  99 // -------------------------------------------------------------------
 100 // -------------------------------------------------------------------
 101 // -------------------------------------------------------------------
 102 // ====================================================== repVertex
 103 void Elements::repVertex (int nh, int nro, Vertex* elt)
 104 {
 105    int addr = nh * pat_nbvertex + nro;
 106    if (tab_vertex[addr] == NULL)
 107        tab_vertex[addr]  = elt;
 108 }
 109 // ====================================================== repVertex
 110 Vertex* Elements::repVertex (int nh, int nro, double px, double py, double pz)
 111 {
 112    int addr = nh * pat_nbvertex + nro;
 113
 114    if (tab_vertex[addr] == NULL)
 115        tab_vertex[addr] = el_root->addVertex (px, py, pz);
 116
 117    return tab_vertex[addr];
 118 }
 119 // ====================================================== repVertex
 120 Vertex* Elements::repVertex (int nh, int nro)
 121 {
 122    int addr = nh * pat_nbvertex + nro;
 123    return tab_vertex[addr];
 124 }
 125 // ====================================================== repEdgeH
 126 void Elements::repEdgeH (int nh, int nro, Edge* elt)
 127 {
 128    int addr = nh * (pat_nbedges + pat_nbvertex) + nro;
 129    if (tab_edge[addr] == NULL)
 130        tab_edge[addr]  = elt;
 131 }
 132 // ====================================================== repEdgeH
 133 Edge* Elements::repEdgeH (int nh, int nro, Vertex* v1, Vertex* v2)
 134 {
 135    int addr = nh * (pat_nbedges + pat_nbvertex) + nro;
 136    if (tab_edge[addr] == NULL)
 137        tab_edge[addr]  = el_root->addEdge (v1, v2);
 138    return tab_edge[addr];
 139 }
 140 // ====================================================== repEdgeH
 141 Edge* Elements::repEdgeH (int nh, int nro)
 142 {
 143    int addr = nh * (pat_nbedges + pat_nbvertex) + nro;
 144    return tab_edge[addr];
 145 }
 146 // ====================================================== repEdgeV
 147 void Elements::repEdgeV (int nh, int nro, Edge* elt)
 148 {
 149    int addr = nh * (pat_nbedges + pat_nbvertex) + pat_nbedges + nro;
 150    if (tab_edge[addr] == NULL)
 151        tab_edge[addr]  = elt;
 152 }
 153 // ====================================================== repEdgeV
 154 Edge* Elements::repEdgeV (int nh, int nro, Vertex* v1, Vertex* v2)
 155 {
 156    int addr = nh * (pat_nbedges + pat_nbvertex) + pat_nbedges + nro;
 157    if (tab_edge[addr] == NULL)
 158        tab_edge[addr]  = el_root->addEdge (v1, v2);
 159    return tab_edge[addr];
 160 }
 161 // ====================================================== repEdgeV
 162 Edge* Elements::repEdgeV (int nh, int nro)
 163 {
 164    int addr = nh * (pat_nbedges + pat_nbvertex) + pat_nbedges + nro;
 165    return tab_edge[addr];
 166 }
 167 // ====================================================== repQuadH
 168 void Elements::repQuadH (int nh, int nro, Quad* elt)
 169 {
 170    int addr = nh * (nbr_orig + pat_nbedges) + nro;
 171    if (tab_quad[addr] == NULL)
 172        tab_quad[addr]  = elt;
 173 }
 174 // ====================================================== repQuadH
 175 Quad* Elements::repQuadH (int nh, int nro)
 176 {
 177    int addr = nh * (nbr_orig + pat_nbedges) + nro;
 178    return tab_quad [addr];
 179 }
 180 // ====================================================== repQuadH
 181 Quad* Elements::repQuadH (int nh, int nro, Edge* e1, Edge* e2, Edge* e3,
 182                                            Edge* e4)
 183 {
 184    int addr = nh * (nbr_orig + pat_nbedges) + nro;
 185    if (tab_quad[addr] == NULL)
 186        tab_quad[addr]  = el_root->addQuad (e1, e2, e3, e4);
 187    return tab_quad [addr];
 188 }
 189 // ====================================================== repQuadV
 190 void Elements::repQuadV (int nh, int nro, Quad* elt)
 191 {
 192    int addr = nh * (nbr_orig + pat_nbedges) + nbr_orig + nro;
 193    if (tab_quad[addr] == NULL)
 194        tab_quad[addr]  = elt;
 195 }
 196 // ====================================================== repQuadV
 197 Quad* Elements::repQuadV (int nh, int nro)
 198 {
 199    int addr = nh * (nbr_orig + pat_nbedges) + nbr_orig + nro;
 200    return tab_quad [addr];
 201 }
 202 // ====================================================== repQuadV
 203 Quad* Elements::repQuadV (int nh, int nro, Edge* e1, Edge* e2, Edge* e3,
 204                                            Edge* e4)
 205 {
 206    int addr = nh * (nbr_orig + pat_nbedges) + nbr_orig + nro;
 207    if (tab_quad[addr] == NULL)
 208        tab_quad[addr]  = el_root->addQuad (e1, e2, e3, e4);
 209    return tab_quad [addr];
 210 }
 211 // ====================================================== repHexa
 212 Hexa* Elements::repHexa (int nh, int nro, Quad* qa, Quad* qb, Quad* qc,
 213                                           Quad* qd, Quad* qe, Quad* qf)
 214 {
 215    int addr = nh * nbr_orig  + nro;
 216    if (tab_hexa[addr] == NULL)
 217        tab_hexa[addr]  = el_root->addHexa (qa, qb, qc, qd, qe, qf);
 218
 219    return tab_hexa [addr];
 220 }
 221 // ====================================================== replaceQuad
 222 // ==== Creation des quads horizontaux
 223 int Elements::replaceQuad (int nh, Qpattern* pat)
 224 {
 225                               // Enregistrement des vertex existant
 226    int size_contour = pat->cont_nbnodes;
 227    for (int nro=0 ; nro<size_contour ; nro++)
 228        {
 229        Vertex* vh = pat->getTargetVertex (nro);
 230        repVertex (nh,  nro, vh);
 231        }
 232
 233    Real3 orig, pmin, pmax, cg, delta, ib, jb, kb;
 234    pat->anaTarget (pmin, pmax, cg);
 235    for (int nc=0 ; nc<DIM3 ; nc++)
 236        delta[nc] = pmax[nc] - pmin[nc];
 237
 238                               // Creation des vertex
 239    pat->getTargetVertex (0)->getPoint (orig);
 240    pat->getTargetVertex (1)->getPoint (kb);
 241    calc_vecteur (orig, kb, ib);
 242    calc_vecteur (orig, cg, jb);
 243    normer_vecteur (ib);
 244    normer_vecteur (jb);
 245    prod_vectoriel (ib, jb, kb);
 246    prod_vectoriel (kb, ib, jb);
 247
 248    for (int nro=0 ; nro<pat->nbr_vertex ; nro++)
 249        {
 250        double lambda = pat->pat_vertex [nro].v_x;
 251        double mu     = pat->pat_vertex [nro].v_y;
 252        double px = pmin[dir_x] + (lambda*ib[dir_x] + mu*jb[dir_x])*delta[dir_x];
 253        double py = pmin[dir_y] + (lambda*ib[dir_y] + mu*jb[dir_y])*delta[dir_y];
 254        double pz = pmin[dir_z] + (lambda*ib[dir_z] + mu*jb[dir_z])*delta[dir_z];
 255        repVertex (nh, nro, px, py, pz);
 256        }
 257                               // Creation des edges horizontaux
 258    for (int nro=0 ; nro<pat->nbr_edges ; nro++)
 259        {
 260        int nv1 = pat->pat_edge [nro].v_amont;
 261        int nv2 = pat->pat_edge [nro].v_aval;
 262        Vertex* v1 = repVertex (nh, nv1);
 263        Vertex* v2 = repVertex (nh, nv2);
 264        repEdgeH (nh, nro, v1, v2);
 265        }
 266                               // Creation des quads horizontaux
 267    for (int nro=0 ; nro<pat->nbr_quads ; nro++)
 268        {
 269        Edge* eda = repEdgeH (nh, pat->pat_quad [nro].q_edge [0]);
 270        Edge* edb = repEdgeH (nh, pat->pat_quad [nro].q_edge [1]);
 271        Edge* edc = repEdgeH (nh, pat->pat_quad [nro].q_edge [2]);
 272        Edge* edd = repEdgeH (nh, pat->pat_quad [nro].q_edge [3]);
 273
 274        repQuadH (nh,  nro, eda, edb, edc, edd);
 275        }
 276    return HOK;
 277 }
 278 // ====================================================== extrudeQuad
 279 // ==== Creation des quads horizontaux
 280 int Elements::extrudeQuad (Qpattern* pat)
 281 {
 282                               // Creation des vertex de niveau 0
 283    for (int nro=0 ; nro<pat->nbr_vertex ; nro++)
 284        {
 285        Vertex* v1 = repVertex (1, nro);
 286        Vertex* v2 = repVertex (2, nro);
 287        double px = 2*v1->getX() - v2->getX();
 288        double py = 2*v1->getY() - v2->getY();
 289        double pz = 2*v1->getZ() - v2->getZ();
 290        repVertex (0, nro, px, py, pz);
 291        }
 292                               // Creation des edges horizontaux
 293    for (int nro=0 ; nro<pat->nbr_edges ; nro++)
 294        {
 295        int nv1 = pat->pat_edge [nro].v_amont;
 296        int nv2 = pat->pat_edge [nro].v_aval;
 297        Vertex* v1 = repVertex (0, nv1);
 298        Vertex* v2 = repVertex (0, nv2);
 299        repEdgeH (0, nro, v1, v2);
 300        }
 301                               // Creation des quads horizontaux
 302    for (int nro=0 ; nro<pat->nbr_quads ; nro++)
 303        {
 304        Edge* eda = repEdgeH (0, pat->pat_quad [nro].q_edge [0]);
 305        Edge* edb = repEdgeH (0, pat->pat_quad [nro].q_edge [1]);
 306        Edge* edc = repEdgeH (0, pat->pat_quad [nro].q_edge [2]);
 307        Edge* edd = repEdgeH (0, pat->pat_quad [nro].q_edge [3]);
 308
 309        repQuadH (0,  nro, eda, edb, edc, edd);
 310        }
 311    return HOK;
 312 }
 313 // ====================================================== replaceHexa
 314 int Elements::replaceHexa (int nh, Qpattern* pat)
 315 {
 316  /* *************************************************
 317    int vnro [4] = { 0, 1, 2, pat->pos_vertex4};
 318    int ednro[4] = { 0, 1, pat->pos_edge3, pat->pos_edge4 };
 319
 320                               // Enregistrement des edges & quads existants
 321    if (nh!=0)
 322        {
 323        for (int nro=0 ; nro<QUAD4 ; nro++)
 324            {
 325            Vertex* vh  = repVertex (nh,   vnro [nro]);
 326            Vertex* vh1 = repVertex (nh+1, vnro [nro]);
 327            Edge*   edv = hexa->findEdge (vh, vh1);
 328            repEdgeV (nh, vnro [nro], edv);
 329
 330            Edge*   edh   = repEdgeH (nh,   ednro [nro]);
 331            Edge*   edh1  = repEdgeH (nh+1, ednro [nro]);
 332            Quad*   quadv = hexa->findQuad (edh, edh1);
 333            repQuadV (nh, ednro [nro], quadv);
 334            }
 335        }
 336     ************************************************* */
 337                               // Creation des edges verticaux
 338    for (int nro=0 ; nro<pat->nbr_vertex ; nro++)
 339        {
 340        Vertex* v1 = repVertex (nh,   nro);
 341        Vertex* v2 = repVertex (nh+1, nro);
 342        repEdgeV (nh, nro, v1, v2);
 343        }
 344                               // Creation des quads verticaux
 345    for (int nro=0 ; nro<pat->nbr_edges ; nro++)
 346        {
 347        int nv1 = pat->pat_edge [nro].v_amont;
 348        int nv2 = pat->pat_edge [nro].v_aval;
 349
 350        Edge* eda = repEdgeH (nh,   nro);
 351        Edge* edb = repEdgeV (nh,   nv1);
 352        Edge* edc = repEdgeH (nh+1, nro);
 353        Edge* edd = repEdgeV (nh,   nv2);
 354        repQuadV (nh,  nro, eda, edb, edc, edd);
 355        }
 356                               // Creation des hexaedres
 357    for (int nro=0 ; nro<pat->nbr_quads ; nro++)
 358        {
 359        Quad* qa = repQuadH (nh,   nro);
 360        Quad* qb = repQuadH (nh+1, nro);
 361        Quad* qc = repQuadV (nh, pat->pat_quad [nro].q_edge [0]);
 362        Quad* qd = repQuadV (nh, pat->pat_quad [nro].q_edge [2]);
 363        Quad* qe = repQuadV (nh, pat->pat_quad [nro].q_edge [1]);
 364        Quad* qf = repQuadV (nh, pat->pat_quad [nro].q_edge [3]);
 365        repHexa (nh, nro, qa, qb, qc, qd, qe, qf);
 366        }
 367
 368    return HOK;
 369 }
 370 END_NAMESPACE_HEXA