src/HEXABLOCK/HexElements_ter.cxx

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