src/HEXABLOCK/HexElements_ter.cxx

   1
   2 // C++ : Table d'hexaedres (Evol Versions 3)
   3
   4 // Copyright (C) 2009-2012  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 "HexVector.hxx"
  26 #include "HexVertex.hxx"
  27 #include "HexEdge.hxx"
  28 #include "HexHexa.hxx"
  29 #include "HexMatrix.hxx"
  30 #include "HexShape.hxx"
  31 #include "HexGlobale.hxx"
  32
  33 #include <cmath>
  34
  35 BEGIN_NAMESPACE_HEXA
  36
  37 void geom_create_circle (double* milieu, double rayon, double* normale,
  38                          double* base, string& brep);
  39 void geom_create_sphere (double* milieu, double radius, string& brep);
  40
  41 void translate_brep  (string& brep, double vdir[], string& trep);
  42 void transfo_brep  (string& brep, Matrix* matrice, string& trep);
  43 void geom_asso_point (Vertex* node);
  44
  45 static bool db=false;
  46
  47 // ====================================================== makeRind
  48 int Elements::makeRind (EnumGrid type, Vertex* center, Vector* vx, Vector* vz,
  49                         double radext, double radint, double radhole,
  50                         Vertex* plorig, double angle, int nr, int na, int nl)
  51 {
  52    double phi1;
  53    int ier = controlRind (type, center, vx, vz, radext, radint, radhole,
  54                           plorig, angle, nr, na, nl, cyl_phi0, phi1);
  55    if (ier!=HOK)
  56       return ier;
  57
  58    resize (type, nr, na, nl);
  59
  60    cyl_radext  = radext;
  61    cyl_radint  = radint;
  62    cyl_radhole = radhole;
  63    cyl_closed  = type==GR_HEMISPHERIC || type==GR_RIND;
  64
  65    double theta  = cyl_closed ? 2*M_PI : M_PI*angle/180;
  66    cyl_dphi      = (phi1-cyl_phi0)/ size_hz;
  67    cyl_dtheta    = theta / size_hy;
  68
  69    int nb_secteurs = cyl_closed ? size_vy-1 : size_vy;
  70
  71    for (int ny=0 ; ny<nb_secteurs ; ny++)
  72        for (int nx=0 ; nx<size_vx ; nx++)
  73            for (int nz=0 ; nz<size_vz ; nz++)
  74                {
  75                double px, py, pz;
  76                getCylPoint (nx, ny, nz, px, py, pz);
  77                Vertex* node = el_root->addVertex (px, py, pz);
  78                setVertex (node, nx, ny, nz);
  79                }
  80    if (cyl_closed)
  81       for (int nx=0 ; nx<size_vx ; nx++)
  82           for (int nz=0 ; nz<size_vz ; nz++)
  83               {
  84               Vertex* node = getVertexIJK (nx, 0, nz);
  85               setVertex (node, nx, size_vy-1, nz);
  86               }
  87
  88    transfoVertices (center, vx, vz);
  89    fillGrid ();
  90    assoCylinder (center, vz, angle);
  91    return HOK;
  92 }
  93 // ====================================================== getCylPoint
  94 int Elements::getCylPoint (int nr, int na, int nh, double& px, double& py,
  95                            double& pz)
  96 {
  97    if (grid_type == GR_CYLINDRIC)
  98       {
  99       px = cyl_radext * cos (na*cyl_dtheta);
 100       py = cyl_radext * sin (na*cyl_dtheta);
 101       pz = cyl_length * nh;
 102       return HOK;
 103       }
 104
 105    bool   rind   = (grid_type == GR_RIND || grid_type == GR_PART_RIND);
 106    double sinphi = sin (cyl_phi0 + nh*cyl_dphi);
 107    double cosphi = cos (cyl_phi0 + nh*cyl_dphi);
 108
 109    double rayon = 0;
 110    if (rind)
 111       {
 112       rayon = cyl_radint + nr*(cyl_radext-cyl_radint)/size_hx;
 113       pz    = rayon*sinphi;
 114       rayon = rayon*cosphi;
 115       }
 116    else
 117       {
 118       pz    = cyl_radext*sinphi;
 119       rayon = cyl_radhole + nr*(cyl_radext*cosphi - cyl_radhole)/size_hx;
 120       rayon = std::max (cyl_radhole, rayon);
 121       }
 122
 123    px = rayon * cos (na*cyl_dtheta);
 124    py = rayon * sin (na*cyl_dtheta);
 125
 126    return HOK;
 127 }
 128 // ====================================================== controlRind
 129 int Elements::controlRind (EnumGrid type, Vertex* cx, Vector* vx, Vector* vz,
 130                            double rext, double rint, double rhole,
 131                            Vertex* px, double angle,
 132                            int nrad, int nang, int nhaut,
 133                            double &phi0, double &phi1)
 134 {
 135    const double Epsil1 = 1e-6;
 136    phi0  = phi1 = 0;
 137
 138    if (cx == NULL || vx == NULL || vz == NULL)
 139       return HERR;
 140
 141    if (nrad<=0 || nang<=0 || nhaut<=0)
 142       return HERR;
 143
 144    if (rext <= 0.0)
 145       return HERR;
 146
 147    if (rint >= rext)
 148       return HERR;
 149
 150    if (rhole > rint)
 151       return HERR;
 152
 153    double nvx = vx->getNorm();
 154    double nvz = vz->getNorm();
 155
 156    if (nvx < Epsil1 || nvz <  Epsil1)
 157       return HERR;
 158
 159    double alpha = asin (rhole/rext);
 160    double beta  = -M_PI*DEMI;
 161    if (type==GR_HEMISPHERIC || type==GR_PART_SPHERIC)
 162        alpha = 2*alpha;
 163
 164    if (px!=NULL)
 165       {
 166           // oh = oa.n/|n|
 167       double oh = ((px->getX() - cx->getX()) * vz->getDx()
 168                 +  (px->getY() - cx->getY()) * vz->getDy()
 169                 +  (px->getZ() - cx->getZ()) * vz->getDz()) / nvz;
 170       if (oh > rext)
 171          return HERR;
 172       else if (oh > -rext)
 173          beta  = asin (oh/rext);
 174       }
 175
 176    phi0 = std::max (alpha - M_PI*DEMI, beta);
 177    phi1 = M_PI*DEMI - alpha;
 178    return HOK;
 179 }
 180 // ====================================================== getHexas
 181 int Elements::getHexas (Hexas& liste)
 182 {
 183    liste.clear ();
 184    for (int nro = 0 ; nro<nbr_hexas ; nro++)
 185        {
 186        Hexa* cell = tab_hexa [nro];
 187        if (cell!=NULL && cell->isValid())
 188           liste.push_back (cell);
 189        }
 190    return HOK;
 191 }
 192 // ====================================================== assoCylinder
 193 void Elements::assoCylinder (Vertex* ori, Vector* normal, double angle)
 194 {
 195    RealVector tangles;
 196    assoCylinders (ori, normal, angle, tangles);
 197 }
 198 // ====================================================== assoCylinders
 199 void Elements::assoCylinders (Vertex* ori, Vector* normal, double angle,
 200                               RealVector& t_angles)
 201 {
 202    int    na      = t_angles.size();
 203    bool   regul   = na == 0;
 204    double alpha   = angle/size_hy;
 205
 206    string brep;
 207    Real3 vk = { normal->getDx(), normal->getDy(), normal->getDz() };
 208    normer_vecteur (vk);
 209
 210    for (int nz=0 ; nz<size_vz ; nz++)
 211        {
 212        for (int nx=0 ; nx<size_vx ; nx++)
 213            {
 214            Vertex* pm = getVertexIJK (nx, 0, nz);
 215            Real3   om = { pm->getX() - ori->getX(),
 216                           pm->getY() - ori->getY(),
 217                           pm->getZ() - ori->getZ() };
 218
 219            double oh     = prod_scalaire (om, vk);
 220            double rayon  = 0;
 221            Real3  ph, hm;
 222            for (int dd=dir_x; dd<=dir_z ; dd++)
 223                {
 224                ph [dd] = ori->getCoord(dd) + oh*vk[dd];
 225                hm [dd] = pm ->getCoord(dd) - ph[dd];
 226                rayon  += hm[dd] * hm[dd];
 227                }
 228
 229            rayon = sqrt (rayon);
 230            geom_create_circle (ph, rayon, vk, hm, brep);
 231
 232            double  pmax = 0;
 233            for (int ny=0 ; ny<size_hy ; ny++)
 234                {
 235                double beta  = regul ? alpha : t_angles [ny];
 236                double pmin = pmax;
 237                pmax  += beta/360;
 238                Edge* edge   =  getEdgeJ  (nx, ny, nz);
 239                Shape* shape = new Shape (brep);
 240                shape-> setBounds (pmin, pmax);
 241                edge->addAssociation (shape);
 242                }
 243            }
 244        }
 245
 246    // Association automatique des vertex non associes -> bph
 247    // Traitement des faces spheriques
 248
 249    Real3 vi = { -vk[dir_x],  -vk[dir_y],  -vk[dir_z] };
 250    Real3 po = { ori->getX(), ori->getY(), ori->getZ() };
 251
 252    switch (grid_type)
 253       {
 254       case GR_HEMISPHERIC  :    // Pour l'exterieur
 255       case GR_PART_SPHERIC :
 256            assoRind (po, vi, size_vx-1);
 257            break;
 258       case GR_PART_RIND    :    // Exterieur + interieur
 259       case GR_RIND         :
 260            assoRind (po, vi, 0);
 261            assoRind (po, vi, size_vx-1);
 262            break;
 263       default :
 264            break;
 265       }
 266    assoResiduelle ();  // Association des sommets residuels
 267 }
 268 // ====================================================== assoRind
 269 // Association des meridiennes
 270 // Creation sphere geometrique + association faces
 271 void Elements::assoRind (double* ori, double* vi, int nx)
 272 {
 273    Real3  vk;
 274    Edges  contour;
 275    string c_brep, s_brep;
 276    Shape  shape (c_brep);
 277    Shapes t_shape;
 278    t_shape.push_back (&shape);
 279
 280    double radius = nx==0 ? cyl_radint : cyl_radext;
 281    double paramin = (cyl_phi0 + M_PI/2) / (2*M_PI);
 282    double paramax = paramin + size_hz*cyl_dphi/(2*M_PI);
 283
 284    paramin = std::max (paramin, 0.0);
 285    paramax = std::min (paramax, 1.0);
 286    geom_create_sphere (ori, radius, s_brep);
 287
 288    int nz1 = size_vz/2;
 289    int nb_secteurs = cyl_closed ? size_vy-1 : size_vy;
 290
 291    for (int ny=0 ; ny<nb_secteurs ; ny++)
 292        {
 293        Vertex* pm = getVertexIJK (nx, ny, nz1);
 294        Real3   vj = { pm->getX(), pm->getY(), pm->getZ() };
 295        prod_vectoriel (vi, vj, vk);
 296        double rayon = cyl_radint + nx*(cyl_radext-cyl_radint)/size_hx;
 297        geom_create_circle (ori, rayon, vk, vi, c_brep);
 298        shape.setBrep   (c_brep);
 299        shape.setBounds (paramin, paramax);
 300        contour.clear ();
 301
 302        for (int nz=0 ; nz<size_hz ; nz++)
 303            {
 304            contour.push_back (getEdgeK  (nx, ny, nz));
 305            Quad* quad = getQuadJK (nx, ny, nz);
 306            if (quad != NULL)
 307               {
 308               Shape* sshape = new Shape (s_brep);
 309               quad->addAssociation (sshape);
 310               }
 311            }
 312        cutAssociation (t_shape, contour, false);
 313        }
 314 }
 315 // ====================================================== assoCircle
 316 // ==== utilise pour les spheres carrees
 317 void Elements::assoCircle (double* center, Edge* ed1, Edge* ed2)
 318 {
 319    Real3 oa,  ob, normal;
 320    Real3 pta, ptb, ptc, ptd;
 321    string brep;
 322
 323 //  Les 2 edges dont les petits cotes d'un rectangle de rapport L/l=sqrt(2)
 324 //  Soit le cercle circonscrit a ce rectangle.
 325 //    * la largeur est balayee par l'angle alpha
 326 //    * la longueur par l'angle beta = pi -alpha
 327
 328    ed1->getVertex(V_AMONT)->getPoint (pta);
 329    ed1->getVertex(V_AVAL )->getPoint (ptb);
 330    ed2->getVertex(V_AMONT)->getPoint (ptc);
 331    ed2->getVertex(V_AVAL )->getPoint (ptd);
 332
 333    double d1 = calc_distance (pta, ptc);
 334    double d2 = calc_distance (pta, ptd);
 335
 336    if (d1 < d2)
 337       {
 338       ed2->getVertex(V_AMONT)->getPoint (ptd);
 339       ed2->getVertex(V_AVAL )->getPoint (ptc);
 340       }
 341
 342    calc_vecteur   (center, pta, oa);
 343    calc_vecteur   (center, ptb, ob);
 344    prod_vectoriel (oa, ob, normal);
 345    double rayon = calc_norme (oa);
 346
 347    geom_create_circle (center, rayon, normal, oa, brep);
 348
 349    Shape* asso1 = new Shape (brep);
 350    Shape* asso2 = new Shape (brep);
 351
 352    const double alpha = atan (sqrt(2)/2)/M_PI; //  angle proche de 70.528 degres
 353    asso1->setBounds (0,   alpha);
 354    asso2->setBounds (0.5, alpha + 0.5);
 355
 356    ed1->addAssociation (asso1);
 357    ed2->addAssociation (asso2);
 358 }
 359 // ====================================================== assoSphere
 360 void Elements::assoSphere (Vertex* ori, Edge* t_edge[], Quad* t_quad[])
 361 {
 362    Real3 center, sommet;
 363    ori->getPoint(center);
 364
 365    assoCircle (center, t_edge [E_AC], t_edge [E_BD]);
 366    assoCircle (center, t_edge [E_AD], t_edge [E_BC]);
 367    assoCircle (center, t_edge [E_AE], t_edge [E_BF]);
 368    assoCircle (center, t_edge [E_AF], t_edge [E_BE]);
 369    assoCircle (center, t_edge [E_CE], t_edge [E_DF]);
 370    assoCircle (center, t_edge [E_CF], t_edge [E_DE]);
 371
 372    t_edge[E_AC]->getVertex(V_AMONT)->getPoint (sommet);
 373    double radius = calc_distance (center, sommet);;
 374
 375    string brep;
 376    geom_create_sphere (center, radius, brep);
 377
 378    for (int nf=0 ; nf < HQ_MAXI ; nf++)
 379        {
 380        Shape* shape = new Shape (brep);
 381        t_quad [nf]->addAssociation (shape);
 382        }
 383
 384    assoResiduelle ();  // Association des sommets residuels
 385 }
 386 // ====================================================== makeSphericalGrid
 387 int Elements::makeSphericalGrid (Vertex* c, double rayon, int nb, double  k)
 388 {
 389    if (nb<=0 || rayon <=Epsil || k <= Epsil || BadElement (c))
 390       {
 391       setError ();
 392       return HERR;
 393       }
 394
 395    resize (GR_SPHERIC, nb);
 396
 397    Vertex* i_node [HV_MAXI];    // Les noeuds de l'hexa englobant
 398    Edge*   i_edge [HE_MAXI];    // Les noeuds de l'hexa englobant
 399    Quad*   i_quad [HQ_MAXI];    // Les noeuds de l'hexa englobant
 400
 401    for (int nro=0 ; nro<HV_MAXI; nro++)
 402        {
 403        double dx = glob->CoordVertex (nro, dir_x) * rayon;
 404        double dy = glob->CoordVertex (nro, dir_y) * rayon;
 405        double dz = glob->CoordVertex (nro, dir_z) * rayon;
 406
 407        i_node [nro] = el_root->addVertex (c->getX ()+dx, c->getY ()+dy,
 408                                           c->getZ ()+dz);
 409        }
 410
 411    for (int nro=0 ; nro<HE_MAXI; nro++)
 412        {
 413        int v1 = glob->EdgeVertex (nro, V_AMONT);
 414        int v2 = glob->EdgeVertex (nro, V_AVAL);
 415        i_edge[nro] = newEdge (i_node[v1], i_node[v2]);
 416
 417        if (db)
 418           {
 419           char nm0[8], nm1 [8], nm2 [8];
 420           printf (" %2d : %s = %s = [%s, %s] = [%d,%d] = [%s,%s]\n", nro,
 421                  glob->namofHexaEdge(nro), i_edge[nro]->getName(nm0),
 422                  glob->namofHexaVertex(v1), glob->namofHexaVertex(v2), v1, v2,
 423                  i_node[v1]->getName(nm1), i_node[v2]->getName(nm2));
 424           }
 425        }
 426
 427    for (int nro=0 ; nro<HQ_MAXI; nro++)
 428        i_quad[nro] = newQuad (i_edge[glob->QuadEdge (nro, E_A)],
 429                               i_edge[glob->QuadEdge (nro, E_B)],
 430                               i_edge[glob->QuadEdge (nro, E_C)],
 431                               i_edge[glob->QuadEdge (nro, E_D)]);
 432
 433    tab_hexa.push_back (newHexa (i_quad[Q_A], i_quad[Q_B], i_quad[Q_C],
 434                                 i_quad[Q_D], i_quad[Q_E], i_quad[Q_F]));
 435    double lambda  = 1;
 436    double dcell   = 1;
 437    for (int niv=0; niv<gr_rayon ; niv++)
 438        {
 439        double lambda0 = lambda;
 440        dcell  *= k;
 441        lambda += dcell;
 442        addStrate (i_quad, i_edge, i_node, c,  lambda/lambda0);
 443        }
 444
 445    assoSphere (c, i_edge, i_quad);
 446    return HOK;
 447 }
 448 // ==================================================== propagateAssociation
 449 int Elements::propagateAssociation (Edge* orig, Edge* dest, Edge* dir)
 450 {
 451    if (revo_lution || orig==NULL || dest==NULL || dir==NULL)
 452       return HERR;
 453
 454    const Shapes& tab_shapes = orig->getAssociations ();
 455    const Shapes& tab_dest   = dest->getAssociations ();
 456    int   nbdest             = tab_dest.size();
 457    int   nbshapes           = tab_shapes.size();
 458    bool  on_edge            = nbshapes!=0 && nbdest==0;
 459
 460    Vertex* vo1 = orig->commonVertex  (dir);
 461    Vertex* vd1 = dest->commonVertex  (dir);
 462    Vertex* vo2 = orig->opposedVertex (vo1);
 463    Vertex* vd2 = dest->opposedVertex (vd1);
 464
 465    if (vo1==NULL || vd1==NULL)
 466       return HERR;
 467
 468    string  trep;
 469    Real3   pa, pb, vdir1, vdir2;
 470    calc_vecteur (vo1->getPoint (pa), vd1->getPoint (pb), vdir1);
 471    calc_vecteur (vo2->getPoint (pa), vd2->getPoint (pb), vdir2);
 472
 473    double dd = calc_distance (vdir1, vdir2);
 474    bool para = dd < 1.0e-3;
 475
 476    if (para && on_edge)
 477       {
 478       for (int nro=0 ; nro<nbshapes ; nro++)
 479           {
 480           Shape* shape  = tab_shapes[nro];
 481           if (shape!=NULL)
 482              {
 483              string brep   = shape->getBrep();
 484              translate_brep (brep, vdir1, trep);
 485              Shape* tshape = new Shape (trep);
 486              tshape->setBounds (shape->debut, shape->fin);
 487              dest->addAssociation (tshape);
 488              }
 489           }
 490       }
 491
 492    double* vdir = vdir1;
 493    for (int nro=V_AMONT ; nro<=V_AVAL ; nro++)
 494        {
 495        Shape* shape = vo1->getAssociation ();
 496        if (shape!=NULL && vd1->getAssociation ()==NULL)
 497           {
 498           string brep   = shape->getBrep();
 499           translate_brep (brep, vdir, trep);
 500           Shape* tshape = new Shape (trep);
 501           vd1->setAssociation (tshape);
 502           }
 503        vo1  = vo2;
 504        vd1  = vd2;
 505        vdir = vdir2;
 506        }
 507
 508    return HOK;
 509 }
 510 // ==================================================== prismAssociation
 511 int Elements::prismAssociation (Edge* orig, Edge* dest, int nh, Edge* dir)
 512 {
 513    if (orig==NULL || dest==NULL || dir==NULL)
 514       return HERR;
 515
 516    updateMatrix (nh);
 517
 518    const Shapes& tab_shapes = orig->getAssociations ();
 519    const Shapes& tab_dest   = dest->getAssociations ();
 520    int   nbdest             = tab_dest.size();
 521    int   nbshapes           = tab_shapes.size();
 522    bool  on_edge            = nbshapes>0 && nbdest==0;
 523
 524    Vertex* vo1 = orig->commonVertex (dir);
 525    Vertex* vd1 = dest->commonVertex (dir);
 526
 527    if (vo1==NULL || vd1==NULL)
 528       return HERR;
 529
 530    string  trep;
 531    Real3   porig, pdest, vdir;
 532    vo1->getPoint (porig);
 533    vd1->getPoint (pdest);
 534    calc_vecteur  (porig, pdest, vdir);
 535
 536    if (on_edge)
 537       {
 538       for (int nro=0 ; nro<nbshapes ; nro++)
 539           {
 540           Shape* shape  = tab_shapes[nro];
 541           if (shape!=NULL)
 542              {
 543              string brep   = shape->getBrep();
 544              //   translate_brep (brep, vdir, trep);
 545              transfo_brep (brep, &gen_matrix, trep);
 546              Shape* tshape = new Shape (trep);
 547              tshape->setBounds (shape->debut, shape->fin);
 548              dest->addAssociation (tshape);
 549              }
 550           }
 551       }
 552
 553    for (int nro=V_AMONT ; nro<=V_AVAL ; nro++)
 554        {
 555        Shape* shape = vo1->getAssociation ();
 556        if (shape!=NULL && vd1->getAssociation ()==NULL)
 557           {
 558           string brep   = shape->getBrep();
 559           //  translate_brep (brep, vdir, trep);
 560           transfo_brep (brep, &gen_matrix, trep);
 561           Shape* tshape = new Shape (trep);
 562           vd1->setAssociation (tshape);
 563           }
 564        vo1 = orig->opposedVertex (vo1);
 565        vd1 = dest->opposedVertex (vd1);
 566        }
 567    return HOK;
 568 }
 569 // ====================================================== assoResiduelle
 570 void Elements::assoResiduelle ()
 571 {
 572    int nbre = tab_vertex.size();
 573    for (int nv=0 ; nv<nbre ; nv++)
 574        {
 575        geom_asso_point (tab_vertex [nv]);
 576        }
 577 }
 578 // ====================================================== moveDisco
 579 void Elements::moveDisco (Hexa* hexa)
 580 {
 581    Real3  center;
 582    Matrix matrix;
 583    hexa->getCenter (center);
 584    matrix.defScale (center, 0.55);
 585
 586    int nbre = tab_vertex.size();
 587    for (int nv=0 ; nv<nbre ; nv++)
 588        {
 589        matrix.perform (tab_vertex [nv]);
 590        }
 591 }
 592 END_NAMESPACE_HEXA