src/HEXABLOCK/HexDocument_v6.cxx

   1
   2 // C++ : Fonctions de creation HexaBlock v6
   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 "HexDocument.hxx"
  24
  25 #include "HexElements.hxx"
  26 #include "HexBiCylinder.hxx"
  27 #include "HexVector.hxx"
  28 #include "HexQuad.hxx"
  29 #include "HexHexa.hxx"
  30
  31 #include "HexGlobale.hxx"
  32
  33 BEGIN_NAMESPACE_HEXA
  34
  35
  36 // ======================================================== make_normale
  37 int make_normale (Quads& tquads, double axe[], double& largeur)
  38 {
  39    double lgmoy = 0;
  40    Real3  cf, cg, vn;
  41    int    nbre = tquads.size();
  42    if (nbre==0)
  43        return HERR;
  44
  45    double coeff = 0;
  46    for (int nro=0 ; nro<nbre ; ++nro)
  47        {
  48        Quad* quad = tquads[nro];
  49        int nbp    = quad->getNbrParents();
  50        if (nbp!=1)
  51           return HERR;
  52                         // Calcul des cotes
  53        for (int nc=0 ; nc<QUAD4 ; ++nc)
  54            lgmoy += quad->getEdge(nc)->getLength ();
  55                         // Calcul de la normale sortante
  56        Hexa* hexa = quad->getParent (0);
  57        hexa->getCenter (cg);
  58        quad->getCenter (cf);
  59        calc_vecteur (cg, cf, vn);
  60        for (int nc=0 ; nc<DIM3 ; ++nc)
  61            axe[nc] = coeff*axe[nc] + vn[nc];
  62        coeff = 1;
  63        }
  64
  65    normer_vecteur (axe);
  66    largeur = lgmoy / (4*nbre);
  67    return HOK;
  68 }
  69 // ======================================================== makeCartesianTop
  70 Elements* Document::makeCartesianTop (int nx, int ny, int nz)
  71 {
  72    DumpStart ("makeCartesianTop", nx << ny << nz );
  73    Elements* grid = new Elements (this);
  74
  75    RealVector tx, ty, tz;
  76    Real3 rx = { 1, 0, 0 };
  77    Real3 ry = { 0, 1, 0 };
  78    Real3 rz = { 0, 0, 1 };
  79
  80    grid->checkSize (1, nx, ny, nz, false);
  81    if (grid->isValid())
  82       {
  83       Vertex* orig = addVertex (0, 0, 0);
  84       for (int nro=1 ; nro<=nx; ++nro) tx.push_back (nro);
  85       for (int nro=1 ; nro<=ny; ++nro) ty.push_back (nro);
  86       for (int nro=1 ; nro<=nz; ++nro) tz.push_back (nro);
  87       grid->makeCartesian (orig, rx, ry, rz, tx, ty, tz);
  88       }
  89
  90    DumpReturn (grid);
  91    return grid;
  92 }
  93 // ======================================================== makeCartesianUni
  94 Elements* Document::makeCartesianUni (Vertex* orig,
  95                                       Vector* vx, Vector* vy, Vector* vz,
  96                                       double  lx, double  ly, double  lz,
  97                                       int     nx, int     ny, int     nz)
  98 {
  99    DumpStart ("makeCartesianUni", orig << vx << vy << vz
 100                                        << lx << ly << lz
 101                                        << nx << ny << nz);
 102
 103    Elements* grid = new Elements (this);
 104
 105    Real3 rx, ry, rz;
 106    grid->checkOrig   (1, orig);
 107    grid->checkSize   (8, nx, ny, nz, false);
 108    grid->checkSystem (2, vx, vy, vz, rx, ry, rz);
 109
 110    if (grid->isValid())
 111       {
 112       RealVector tx, ty, tz;
 113       double dx = lx / nx;
 114       double dy = ly / ny;
 115       double dz = lz / nz;
 116       for (int nro=1 ; nro<=nx; ++nro) tx.push_back (nro*dx);
 117       for (int nro=1 ; nro<=ny; ++nro) ty.push_back (nro*dy);
 118       for (int nro=1 ; nro<=nz; ++nro) tz.push_back (nro*dz);
 119
 120       grid->makeCartesian (orig, rx, ry, rz, tx, ty, tz);
 121       }
 122
 123    DumpReturn (grid);
 124    return grid;
 125 }
 126 // ======================================================== makeCartesian
 127 Elements* Document::makeCartesian (Vertex*    orig,
 128                                    Vector*    vx, Vector*    vy, Vector*    vz,
 129                                    RealVector tx, RealVector ty, RealVector tz)
 130 {
 131    DumpStart ("makeCartesian", orig << vx << vy << vz << tx << ty << tz);
 132
 133    Elements* grid = new Elements (this);
 134
 135    Real3 rx, ry, rz;
 136    grid->checkOrig   (1, orig);
 137    grid->checkSystem (2, vx, vy, vz, rx, ry, rz);
 138    grid->checkVector (5, tx);
 139    grid->checkVector (6, ty);
 140    grid->checkVector (7, tz);
 141
 142    if (grid->isValid())
 143        grid->makeCartesian (orig, rx, ry, rz, tx, ty, tz);
 144
 145    DumpReturn (grid);
 146    return grid;
 147 }
 148 // ======================================================== makeCylinderTop
 149 Elements* Document::makeCylinderTop (int nr, int na, int nh)
 150 {
 151    DumpStart ("makeCylinderTop", nr << na << nh);
 152
 153    Elements* grid = new Elements (this);
 154
 155    RealVector tray, tang, thaut;
 156    Real3 rx = { 1, 0, 0 };
 157    Real3 rz = { 0, 0, 1 };
 158
 159    grid->checkSize (1, nr, na, nh, true);
 160    if (grid->isValid())
 161       {
 162       double da = 360.0 / na;
 163       for (int nro=0 ; nro<=nr; ++nro) tray .push_back (nro+1);
 164       for (int nro=0 ; nro<=na; ++nro) tang .push_back (nro*da);
 165       for (int nro=0 ; nro<=nh; ++nro) thaut.push_back (nro);
 166       grid->makeCylinder (NULL, rx, rz, tray, tang, thaut, true);
 167       }
 168
 169    DumpReturn (grid);
 170    return grid;
 171 }
 172 // ======================================================== makeCylinderUni
 173 Elements* Document::makeCylinderUni (Vertex* orig, Vector* vx, Vector* vz,
 174                            double rint, double rext, double ang, double haut,
 175                            int nr, int na, int nh)
 176 {
 177    DumpStart ("makeCylinderUni", orig << vx << vz << rint << rext
 178                                       << ang << haut << nr << na << nh);
 179
 180    Elements* grid = new Elements (this);
 181
 182    Real3 rx, rz;
 183    grid->checkOrig   (1, orig);
 184    grid->checkSize   (8, nr, na, nh, true);
 185    grid->checkPipe   (4, rint, rext, ang, haut);
 186    grid->checkSystem (2, vx, vz, rx, rz);
 187
 188    if (grid->isValid())
 189       {
 190       double dray  = (rext-rint) / nr;
 191       double dang  = ang  / na;
 192       double dhaut = haut / nh;
 193       RealVector tang, tray, thaut;
 194
 195       for (int nro=0 ; nro<=nr; ++nro) tray .push_back (rint + nro*dray);
 196       for (int nro=0 ; nro<=na; ++nro) tang .push_back (nro*dang);
 197       for (int nro=0 ; nro<=nh; ++nro) thaut.push_back (nro*dhaut);
 198
 199       grid->makeCylinder (orig, rx, rz, tray, tang, thaut, true);
 200       }
 201
 202    DumpReturn (grid);
 203    return grid;
 204 }
 205 // ======================================================== makeCylinder
 206 Elements* Document::makeCylinder (Vertex* orig, Vector* vx, Vector* vz,
 207                         RealVector tray, RealVector tang, RealVector thaut)
 208 {
 209    DumpStart ("makeCylinder", orig << vx << vz << tray << tang << thaut);
 210
 211    Elements* grid = new Elements (this);
 212
 213    Real3 rx, rz;
 214    grid->checkOrig   (1, orig);
 215    grid->checkSystem (2, vx, vz, rx, rz);
 216    grid->checkVector (4, tray);
 217    grid->checkVector (5, tang,  3);
 218    grid->checkVector (6, thaut);
 219
 220    if (grid->isValid())
 221       grid->makeCylinder (orig, rx, rz, tray, tang, thaut, true);
 222
 223    DumpReturn (grid);
 224    return grid;
 225 }
 226 // ======================================================== makePipeTop
 227 Elements* Document::makePipeTop (int nr, int na, int nh)
 228 {
 229    DumpStart ("makePipeTop", nr << na << nh);
 230
 231    Elements* grid = new Elements (this);
 232
 233    RealVector tray, tang, thaut;
 234    Real3 rx = { 1, 0, 0 };
 235    Real3 rz = { 0, 0, 1 };
 236
 237    grid->checkSize (1,nr, na, nh, true);
 238    if (grid->isValid())
 239       {
 240       double da = 360.0 / na;
 241       for (int nro=0 ; nro<=nr; ++nro) tray .push_back (nro+1);
 242       for (int nro=0 ; nro<=na; ++nro) tang .push_back (nro*da);
 243       for (int nro=0 ; nro<=nh; ++nro) thaut.push_back (nro);
 244       grid->makeCylinder (NULL, rx, rz, tray, tang, thaut, false);
 245       }
 246
 247    DumpReturn (grid);
 248    return grid;
 249 }
 250 // ======================================================== makePipeUni
 251 Elements* Document::makePipeUni (Vertex* orig, Vector* vx, Vector* vz,
 252                        double rint, double rext, double angle, double hauteur,
 253                        int nr, int na, int nh)
 254 {
 255    DumpStart ("makePipeUni", orig << vx << vz << rint << rext
 256                                   << angle << hauteur << nr << na << nh);
 257
 258    Elements* grid = new Elements (this);
 259
 260    Real3 rx, rz;
 261    RealVector tray, tang, thaut;
 262    grid->checkOrig   (1, orig);
 263    grid->checkSize   (8, nr, na, nh, true);
 264    grid->checkPipe   (4, rint, rext, angle, hauteur);
 265    grid->checkSystem (1, vx, vz, rx, rz);
 266
 267    if (grid->isValid())
 268       {
 269       double dray  = (rext-rint) / nr;
 270       double dang  = angle   / na;
 271       double dhaut = hauteur / nh;
 272       for (int nro=0 ; nro<=nr; ++nro) tray .push_back (rint + nro*dray);
 273       for (int nro=0 ; nro<=na; ++nro) tang .push_back (nro  * dang);
 274       for (int nro=0 ; nro<=nh; ++nro) thaut.push_back (nro  * dhaut);
 275
 276       grid->makeCylinder (orig, rx, rz, tray, tang, thaut, false);
 277       }
 278    DumpReturn (grid);
 279    return grid;
 280 }
 281 // ======================================================== makePipe
 282 Elements* Document::makePipe (Vertex* orig, Vector* vx, Vector* vz,
 283                        RealVector tray, RealVector tang, RealVector thaut)
 284 {
 285    DumpStart ("makePipe", orig << vx << vz << tray << tang << thaut);
 286
 287    Elements* grid = new Elements (this);
 288
 289    Real3 rx, rz;
 290    grid->checkOrig   (1, orig);
 291    grid->checkSystem (2, vx, vz, rx, rz);
 292    grid->checkVector (4, tray,  2);
 293    grid->checkVector (5, tang,  3);
 294    grid->checkVector (6, thaut, 2);
 295
 296    if (grid->isValid())
 297       grid->makeCylinder (orig, rx, rz, tray, tang, thaut, false);
 298
 299    DumpReturn (grid);
 300    return grid;
 301 }
 302 // ======================================================== makeSphericalTop
 303 Elements* Document::makeSphericalTop (int nbre, int crit)
 304 {
 305    DumpStart ("makeSphericalTop", nbre << crit);
 306    Elements* grid = new Elements (this);
 307
 308    Real3 orig = { 0, 0, 0 };
 309    Real3 rx   = { 1, 0, 0 };
 310    Real3 rz   = { 0, 0, 1 };
 311
 312    if (nbre<2)
 313       {
 314       grid->setError (HERR);
 315       Mess << "Argument number 1 :";
 316       Mess << "Nbr slices must be greather or equal than 2";
 317       }
 318
 319    if (grid->isValid())
 320       {
 321       RealVector trayon;
 322       for (int nro=1 ; nro<=nbre; ++nro) trayon.push_back (nro);
 323       grid->makeSpherical (orig, rx, rz, trayon, crit);
 324       }
 325
 326    DumpReturn (grid);
 327    return grid;
 328 }
 329 // ======================================================== makeSphericalUni
 330 Elements* Document::makeSphericalUni (Vertex* center, Vector* vx, Vector* vz,
 331                                       double rayon, int nbre, int crit)
 332 {
 333    DumpStart ("makeSphericalUni", center << vx << vz << rayon << nbre << crit);
 334    Elements* grid = new Elements (this);
 335
 336    Real3 rx, rz;
 337    grid->checkOrig   (1, center);
 338    grid->checkSystem (2, vx, vz, rx, rz);
 339
 340    if (rayon<=1e-5)
 341       {
 342       grid->setError (HERR);
 343       Mess << "Argument number 4 :";
 344       Mess << "Radius must be positive";
 345       }
 346
 347    if (nbre<2)
 348       {
 349       grid->setError (HERR);
 350       Mess << "Argument number 5 :";
 351       Mess << "Nbr slices must be greather or equal than 2";
 352       }
 353
 354    if (grid->isValid())
 355       {
 356       Real3      orig;
 357       RealVector trayon;
 358       center->getPoint (orig);
 359       for (int nro=1 ; nro<=nbre; ++nro) trayon.push_back (nro*rayon);
 360       grid->makeSpherical (orig, rx, rz, trayon, crit);
 361       }
 362
 363    DumpReturn (grid);
 364    return grid;
 365 }
 366 // ======================================================== makeSpherical
 367 Elements* Document::makeSpherical (Vertex* center, Vector* vx, Vector* vz,
 368                                    RealVector trayon, int crit)
 369 {
 370    DumpStart ("makeSpherical", center << vx << vz << trayon << crit);
 371    Elements* grid = new Elements (this);
 372
 373    Real3 rx, rz;
 374    grid->checkOrig   (1, center);
 375    grid->checkSystem (2, vx, vz, rx, rz);
 376    grid->checkVector (4, trayon, 2);
 377
 378    if (grid->isValid())
 379       {
 380       Real3 orig;
 381       center->getPoint (orig);
 382       grid->makeSpherical (orig, rx, rz, trayon, crit);
 383       }
 384
 385
 386    DumpReturn (grid);
 387    return grid;
 388 }
 389 // ======================================================== makeSphereTop
 390 Elements* Document::makeSphereTop (int nr, int na, int nh)
 391 {
 392    DumpStart ("makeSphereTop", nr << na << nh);
 393    Elements* grid = new Elements (this);
 394
 395    Real3 orig = { 0, 0, 0 };
 396    Real3 rx   = { 1, 0, 0 };
 397    Real3 rz   = { 0, 0, 1 };
 398
 399    double rtrou = 1;
 400    double rext  = 10;
 401    double angle = 360;
 402    double phi0, phi1;
 403
 404    grid->checkPhi (true, orig, rz, rext, rtrou, NULL, phi0, phi1);
 405
 406    if (nr<2)
 407       {
 408       grid->setError (HERR);
 409       Mess << "Argument number 1 :";
 410       Mess << "Nbr slices must be greather or equal than 2";
 411       }
 412
 413    if (na<3)
 414       {
 415       grid->setError (HERR);
 416       Mess << "Argument number 2 :";
 417       Mess << "Nbr sectors must be greather or equal than 3";
 418       }
 419
 420    if (nh<2)
 421       {
 422       grid->setError (HERR);
 423       Mess << "Argument number 3 :";
 424       Mess << "Nbr slices must be greather or equal than 2";
 425       }
 426
 427    if (grid->isValid())
 428       {
 429       double dr = (rext-rtrou)/nr;
 430       double da = angle / na;
 431       double dh = (phi1-phi0) / nh;
 432       RealVector tray, tang, thaut;
 433       for (int nro=0 ; nro<=nr; ++nro) tray .push_back (rtrou + nro*dr);
 434       for (int nro=0 ; nro<=na; ++nro) tang .push_back (nro*da);
 435       for (int nro=0 ; nro<=nh; ++nro) thaut.push_back (phi0  + nro*dh);
 436       grid->makeRind (GR_HEMISPHERIC, orig, rx, rz, tray, tang, thaut);
 437       }
 438
 439    DumpReturn (grid);
 440    return grid;
 441 }
 442 // ======================================================== makeSphereUni
 443 Elements* Document::makeSphereUni (Vertex* center, Vector* vx, Vector* vz,
 444                             double rtrou, double rext, double ang,
 445                             Vertex* hplan, int nr, int na, int nh)
 446 {
 447    DumpStart ("makeSphereUni", center << vx << vz << rtrou << rext << ang
 448                                       << hplan << nr << na << nh);
 449    Elements* grid = new Elements (this);
 450
 451    Real3 rx, rz;
 452    grid->checkOrig   (1, center);
 453    grid->checkSystem (2, vx, vz, rx, rz);
 454
 455    Real3 orig = { 0, 0, 0 };
 456    if (center != NULL)
 457        center->getPoint (orig);
 458
 459    double phi0, phi1;
 460    grid->checkPhi (true, orig, rz, rext, rtrou, hplan, phi0, phi1);
 461
 462    if (rtrou<=1e-5)
 463       {
 464       grid->setError (HERR);
 465       Mess << "Argument number 4 :";
 466       Mess << "Hole radius must be positive";
 467       }
 468
 469    if (rext<rtrou)
 470       {
 471       grid->setError (HERR);
 472       Mess << "Argument number 5 :";
 473       Mess << "External radius must be greather than hole radius";
 474       }
 475
 476    if (ang<1.0)
 477       {
 478       grid->setError (HERR);
 479       Mess << "Argument number 6 :";
 480       Mess << "Angle must be positive";
 481       }
 482
 483    if (nr<2)
 484       {
 485       grid->setError (HERR);
 486       Mess << "Argument number 7 :";
 487       Mess << "Nbr slices must be greather or equal than 2";
 488       }
 489
 490    if (na<3)
 491       {
 492       grid->setError (HERR);
 493       Mess << "Argument number 8 :";
 494       Mess << "Nbr sectors must be greather or equal than 3";
 495       }
 496
 497    if (nh<2)
 498       {
 499       grid->setError (HERR);
 500       Mess << "Argument number 9 :";
 501       Mess << "Nbr slices must be greather or equal than 2";
 502       }
 503
 504    EnumGrid type = ang>= 359.9 ? GR_HEMISPHERIC : GR_PART_SPHERIC;
 505    if (grid->isValid())
 506       {
 507       double dr = (rext-rtrou)/nr;
 508       double dh = (phi1-phi0) /nh;
 509       double da = ang / na;
 510       RealVector tray, tang, thaut;
 511       for (int nro=0 ; nro<=nr; ++nro) tray .push_back (rtrou + nro*dr);
 512       for (int nro=0 ; nro<=na; ++nro) tang .push_back (nro*da);
 513       for (int nro=0 ; nro<=nh; ++nro) thaut.push_back (phi0  + nro*dh);
 514       grid->makeRind (type, orig, rx, rz, tray, tang, thaut);
 515       }
 516
 517    DumpReturn (grid);
 518    return grid;
 519 }
 520 // ======================================================== makeSphere
 521 Elements* Document::makeSphere (Vertex* center, Vector* vx, Vector* vz,
 522                             RealVector tray, RealVector tang, RealVector thaut)
 523 {
 524    DumpStart ("makeSphere", center << vx << vz << tray << tang << thaut );
 525    Elements* grid = new Elements (this);
 526
 527    Real3 rx, rz;
 528    grid->checkOrig   (1, center);
 529    grid->checkSystem (2, vx, vz, rx, rz);
 530    grid->checkVector (4, tray,  2);
 531    grid->checkVector (5, tang,  3);
 532    grid->checkVector (6, thaut, 2);
 533
 534    int    nb    = tang.size()-1;
 535    double angle = nb <0 ? 0 : tang[nb] ;
 536    EnumGrid type = angle >= 359.9 ? GR_HEMISPHERIC : GR_PART_SPHERIC;
 537    if (grid->isValid())
 538       {
 539       Real3 orig;
 540       center->getPoint (orig);
 541       grid->makeRind (type, orig, rx, rz, tray, tang, thaut);
 542       }
 543
 544
 545    DumpReturn (grid);
 546    return grid;
 547 }
 548
 549 // ======================================================== makeRindTop
 550 Elements* Document::makeRindTop (int nr, int na, int nh)
 551 {
 552    DumpStart ("makeRindTop", nr << na << nh);
 553    Elements* grid = new Elements (this);
 554
 555    Real3 orig = { 0, 0, 0 };
 556    Real3 rx   = { 1, 0, 0 };
 557    Real3 rz   = { 0, 0, 1 };
 558
 559    double rtrou = 1;
 560    double rint  = 8;
 561    double rext  = 10;
 562    double angle = 360;
 563    double phi0, phi1;
 564
 565    grid->checkPhi (false, orig, rz, rext, rtrou, NULL, phi0, phi1);
 566
 567    if (nr<2)
 568       {
 569       grid->setError (HERR);
 570       Mess << "Argument number 1 :";
 571       Mess << "Nbr slices must be greather or equal than 2";
 572       }
 573
 574    if (na<3)
 575       {
 576       grid->setError (HERR);
 577       Mess << "Argument number 2 :";
 578       Mess << "Nbr sectors must be greather or equal than 3";
 579       }
 580
 581    if (nh<2)
 582       {
 583       grid->setError (HERR);
 584       Mess << "Argument number 3 :";
 585       Mess << "Nbr slices must be greather or equal than 2";
 586       }
 587
 588    if (grid->isValid())
 589       {
 590       double dr = (rext-rint)/nr;
 591       double da = angle / na;
 592       double dh = (phi1-phi0) / nh;
 593       RealVector tray, tang, thaut;
 594       tray .push_back (rtrou);
 595       for (int nro=0 ; nro<=nr; ++nro) tray .push_back (rint + nro*dr);
 596       for (int nro=0 ; nro<=na; ++nro) tang .push_back (nro*da);
 597       for (int nro=0 ; nro<=nh; ++nro) thaut.push_back (phi0  + nro*dh);
 598       grid->makeRind (GR_RIND, orig, rx, rz, tray, tang, thaut);
 599       }
 600
 601    DumpReturn (grid);
 602    return grid;
 603 }
 604 // ======================================================== makeRindUni
 605 Elements* Document::makeRindUni (Vertex* center, Vector* vx, Vector* vz,
 606                           double rtrou, double rint, double rext,
 607                           double ang, Vertex* hplan, int nr, int na, int nh)
 608 {
 609    DumpStart ("makeRindUni", center << vx << vz << rtrou << rint << rext << ang
 610                                     << hplan << nr << na << nh);
 611    Elements* grid = new Elements (this);
 612
 613    Real3 rx, rz;
 614    grid->checkOrig   (1, center);
 615    grid->checkSystem (2, vx, vz, rx, rz);
 616    grid->checkDiams  (4, rtrou, rint);
 617    grid->checkDiams  (5, rint,  rext);
 618
 619    Real3 orig = { 0, 0, 0 };
 620    if (center != NULL)
 621        center->getPoint (orig);
 622
 623    double phi0, phi1;
 624    grid->checkPhi (false, orig, rz, rext, rtrou, hplan, phi0, phi1);
 625
 626    if (rtrou<=1e-5)
 627       {
 628       grid->setError (HERR);
 629       Mess << "Argument number 4 :";
 630       Mess << "Hole radius must be positive";
 631       }
 632
 633    if (rint<rtrou)
 634       {
 635       grid->setError (HERR);
 636       Mess << "Argument number 5 :";
 637       Mess << "Internal radius must be greather than hole radius";
 638       }
 639
 640    if (rext<rtrou)
 641       {
 642       grid->setError (HERR);
 643       Mess << "Argument number 6 :";
 644       Mess << "External radius must be greather than internal radius";
 645       }
 646
 647    if (ang<1.0)
 648       {
 649       grid->setError (HERR);
 650       Mess << "Argument number 7 :";
 651       Mess << "Angle must be positive";
 652       }
 653
 654    if (nr<2)
 655       {
 656       grid->setError (HERR);
 657       Mess << "Argument number 8 :";
 658       Mess << "Nbr slices must be greather or equal than 2";
 659       }
 660
 661    if (na<3)
 662       {
 663       grid->setError (HERR);
 664       Mess << "Argument number 9 :";
 665       Mess << "Nbr sectors must be greather or equal than 3";
 666       }
 667
 668    if (nh<2)
 669       {
 670       grid->setError (HERR);
 671       Mess << "Argument number 10 :";
 672       Mess << "Nbr slices must be greather or equal than 2";
 673       }
 674
 675    EnumGrid type = ang>= 359.9 ? GR_RIND : GR_PART_RIND;
 676    if (grid->isValid())
 677       {
 678       double dr = (rext-rint)/nr;
 679       double dh = (phi1-phi0) /nh;
 680       double da = ang / na;
 681       RealVector tray, tang, thaut;
 682       tray .push_back (rtrou);
 683       for (int nro=0 ; nro<=nr; ++nro) tray .push_back (rint + nro*dr);
 684       for (int nro=0 ; nro<=na; ++nro) tang .push_back (nro*da);
 685       for (int nro=0 ; nro<=nh; ++nro) thaut.push_back (phi0 + nro*dh);
 686       grid->makeRind (type, orig, rx, rz, tray, tang, thaut);
 687       }
 688
 689    DumpReturn (grid);
 690    return grid;
 691 }
 692 // ======================================================== makeRind
 693 Elements* Document::makeRind (Vertex* center, Vector* vx, Vector* vz,
 694                           RealVector tray, RealVector tang, RealVector thaut)
 695 {
 696    DumpStart ("makeRind", center << vx << vz << tray << tang << thaut );
 697    Elements* grid = new Elements (this);
 698
 699    Real3 rx, rz;
 700    grid->checkOrig   (1, center);
 701    grid->checkSystem (2, vx, vz, rx, rz);
 702    grid->checkVector (4, tray,  2);
 703    grid->checkVector (5, tang,  3);
 704    grid->checkVector (6, thaut, 2);
 705
 706    int    nb    = tang.size()-1;
 707    double angle = nb <0 ? 0 : tang[nb] ;
 708    EnumGrid type = angle >= 359.9 ? GR_RIND : GR_PART_RIND;
 709    if (grid->isValid())
 710       {
 711       Real3 orig;
 712       center->getPoint (orig);
 713       grid->makeRind (type, orig, rx, rz, tray, tang, thaut);
 714       }
 715
 716
 717    DumpReturn (grid);
 718    return grid;
 719 }
 720 // ======================================================== makeCylinders
 721 BiCylinder* Document::makeCylinders (Vertex* ori1, Vector* vz1, double ray1,
 722                 double h1, Vertex* ori2, Vector* vz2, double ray2, double h2)
 723 {
 724    DumpStart ("makeCylinders", ori1 << vz1 << ray1 << h1 << ori2 << vz2
 725                                            << ray2<< h2);
 726    BiCylinder* grid = new BiCylinder (this);
 727
 728    Real3 rz1, rz2;
 729    grid->checkOrig   (1, ori1);
 730    grid->checkOrig   (5, ori2);
 731    grid->checkSystem (2, vz1, vz2, rz1, rz2);
 732
 733    if (grid->isValid())
 734       {
 735       if (ray1 < ray2)
 736           grid->makeCylinders (ori1, rz1, ray1, h1, ori2, rz2, ray2, h2);
 737       else
 738           grid->makeCylinders (ori2, rz2, ray2, h2, ori1, rz1, ray1, h1);
 739       }
 740
 741
 742    DumpReturn (grid);
 743    return grid;
 744 }
 745 // ======================================================== makePipes
 746 BiCylinder* Document::makePipes (Vertex* ori1, Vector* vz1, double rint1,
 747                                  double rext1, double h1,
 748       Vertex* ori2, Vector* vz2, double rint2, double rext2, double h2)
 749 {
 750    DumpStart ("makePipes", ori1 << vz1 << rint1 << rext1 << h1
 751                         << ori2 << vz2 << rint2 << rext2 << h2);
 752    BiCylinder* grid = new BiCylinder (this);
 753
 754    Real3 rz1, rz2;
 755    grid->checkOrig   (1, ori1);
 756    grid->checkOrig   (6, ori2);
 757    grid->checkSystem (2, vz1, vz2, rz1, rz2);
 758    grid->checkDiams  (3, rint1,  rext1);
 759    grid->checkDiams  (8, rint2,  rext2);
 760
 761    if (grid->isValid())
 762       {
 763       if (rext1 < rext2)
 764           grid->makePipes (ori1,rz1,rint1,rext1,h1, ori2,rz2,rint2,rext2,h2);
 765       else
 766           grid->makePipes (ori2,rz2,rint2,rext2,h2, ori1,rz1,rint1,rext1,h1);
 767       }
 768
 769    DumpReturn (grid);
 770    return grid;
 771 }
 772 // ------------------------------------------------------------------------
 773 // ------------------------------------------------------------------------
 774 // ------------------------------------------------------------------------
 775 // ========================================================= extrudeQuadTop
 776 Elements* Document::extrudeQuadTop (Quad* start, int nbre)
 777 {
 778    DumpStart ("extrudeQuadTop", start << nbre);
 779
 780    Quads      tstart;
 781    RealVector tlen;
 782    Real3      axe  = { 0, 0, 1 };
 783    double     larg = 1;
 784    Elements*  grid = new Elements (this);
 785
 786    tstart.push_back  (start);
 787    grid->checkQuad   (start);
 788    make_normale (tstart, axe, larg);
 789    for (int nro=0 ; nro<nbre; ++nro) tlen.push_back ((nro+1)*larg);
 790
 791    if (grid->isValid())
 792       {
 793       grid->extrudeQuads  (tstart, axe, tlen);
 794       }
 795
 796    DumpReturn (grid);
 797    return grid;
 798 }
 799 // ========================================================= extrudeQuadUni
 800 Elements* Document::extrudeQuadUni (Quad* start, Vector* axis, double len,
 801                                     int nbre)
 802 {
 803    DumpStart ("extrudeQuadUni", start << axis << len << nbre);
 804
 805    Quads      tstart;
 806    RealVector tlen;
 807    Real3      axe;
 808    Elements*  grid = new Elements (this);
 809
 810    grid->checkQuad  (start);
 811    grid->checkAxis  (axis, axe);
 812
 813    tstart.push_back (start);
 814
 815    for (int nro=0 ; nro<nbre; ++nro) tlen.push_back ((nro+1)*len);
 816
 817    if (grid->isValid())
 818       {
 819       grid->extrudeQuads  (tstart, axe, tlen);
 820       }
 821
 822    DumpReturn (grid);
 823    return grid;
 824 }
 825 // ========================================================= extrudeQuad
 826 Elements* Document::extrudeQuad (Quad* start, Vector* axis, RealVector tlen)
 827 {
 828    DumpStart ("extrudeQuad", start << axis << tlen);
 829
 830    Quads     tstart;
 831    Real3     axe;
 832    Elements* grid = new Elements (this);
 833
 834    grid->checkQuad   (start);
 835    grid->checkAxis   (axis, axe);
 836    grid->checkVector (3, tlen);
 837
 838    tstart.push_back (start);
 839
 840    if (grid->isValid())
 841       {
 842       grid->extrudeQuads  (tstart, axe, tlen);
 843       }
 844
 845    DumpReturn (grid);
 846    return grid;
 847 }
 848
 849 // ========================================================= extrudeQuadsTop
 850 Elements* Document::extrudeQuadsTop (Quads tstart, int nbre)
 851 {
 852    DumpStart ("extrudeQuadsTop", tstart << nbre);
 853
 854    RealVector tlen;
 855    Real3      axe  = { 0, 0, 1 };
 856    double     larg = 1;
 857    Elements* grid = new Elements (this);
 858
 859    grid->checkQuads (tstart);
 860    make_normale     (tstart, axe, larg);
 861    for (int nro=0 ; nro<nbre; ++nro) tlen.push_back ((nro+1)*larg);
 862
 863    if (grid->isValid())
 864       {
 865       grid->extrudeQuads  (tstart, axe, tlen);
 866       }
 867
 868    DumpReturn (grid);
 869    return grid;
 870 }
 871 // ========================================================= extrudeQuadsUni
 872 Elements* Document::extrudeQuadsUni (Quads tstart, Vector* axis, double len,
 873                                      int nbre)
 874 {
 875    DumpStart ("extrudeQuadsUni", tstart << axis << len << nbre);
 876
 877    RealVector tlen;
 878    Real3      axe;
 879    Elements*  grid = new Elements (this);
 880
 881    grid->checkQuads (tstart);
 882    grid->checkAxis  (axis, axe);
 883
 884    for (int nro=0 ; nro<nbre; ++nro) tlen.push_back ((nro+1)*len);
 885
 886    if (grid->isValid())
 887       {
 888       grid->extrudeQuads  (tstart, axe, tlen);
 889       }
 890
 891    DumpReturn (grid);
 892    return grid;
 893 }
 894 // ========================================================= extrudeQuads
 895 Elements* Document::extrudeQuads (Quads tstart, Vector* axis, RealVector tlen)
 896 {
 897    DumpStart ("extrudeQuads", tstart << axis << tlen);
 898
 899    Real3     axe;
 900    Elements* grid = new Elements (this);
 901
 902    grid->checkQuads  (tstart);
 903    grid->checkAxis   (axis, axe);
 904    grid->checkVector (3, tlen);
 905
 906    if (grid->isValid())
 907       {
 908       grid->extrudeQuads  (tstart, axe, tlen);
 909       }
 910
 911    DumpReturn (grid);
 912    return grid;
 913 }
 914 // ========================================================= revolutionQuadUni
 915 Elements* Document::revolutionQuadUni (Quad* start, Vertex* center,
 916                                        Vector* axis, double angle, int nbre)
 917 {
 918    DumpStart ("revolutionQuadUni", start << center << axis << angle << nbre);
 919
 920    Quads      tstart;
 921    RealVector angles;
 922    Real3      axe;
 923    Elements* grid = new Elements (this);
 924    double   alpha = angle / std::max (nbre, 1);
 925
 926    grid->checkQuad (start);
 927    grid->checkOrig (2, center);
 928    grid->checkAxis (axis, axe);
 929
 930    tstart.push_back (start);
 931    for (int nro=0 ; nro<nbre; ++nro) angles.push_back ((nro+1)*alpha);
 932
 933    if (grid->isValid())
 934       {
 935       grid->extrudeQuads  (tstart, axe, angles, true, center);
 936       }
 937
 938    DumpReturn (grid);
 939    return grid;
 940 }
 941 // ========================================================= revolutionQuad
 942 Elements* Document::revolutionQuad (Quad* start, Vertex* center, Vector* axis,
 943                                     RealVector angles)
 944 {
 945    DumpStart ("revolutionQuads", start << center << axis << angles);
 946
 947    Quads      tstart;
 948    Real3      axe;
 949    Elements*  grid = new Elements (this);
 950
 951    grid->checkQuad (start);
 952    grid->checkOrig (2, center);
 953    grid->checkAxis (axis, axe);
 954
 955    tstart.push_back (start);
 956
 957    if (grid->isValid())
 958       {
 959       grid->extrudeQuads  (tstart, axe, angles, true, center);
 960       }
 961
 962    DumpReturn (grid);
 963    return grid;
 964 }
 965 // ========================================================= revolutionQuadsUni
 966 Elements* Document::revolutionQuadsUni (Quads tstart, Vertex* center,
 967                                         Vector* axis, double angle, int nbre)
 968 {
 969    DumpStart ("revolutionQuadsUni", tstart << center << axis << angle << nbre);
 970
 971    RealVector angles;
 972    Real3      axe;
 973    Elements* grid  = new Elements (this);
 974    double    alpha = angle / std::max (nbre, 1);
 975
 976    grid->checkQuads (tstart);
 977    grid->checkOrig  (2, center);
 978    grid->checkAxis  (axis, axe);
 979
 980    for (int nro=0 ; nro<nbre; ++nro) angles.push_back ((nro+1)*alpha);
 981
 982    if (grid->isValid())
 983       {
 984       grid->extrudeQuads  (tstart, axe, angles, true, center);
 985       }
 986
 987    DumpReturn (grid);
 988    return grid;
 989 }
 990 // ========================================================= revolutionQuads
 991 Elements* Document::revolutionQuads (Quads tstart, Vertex* center, Vector* axis,
 992                                      RealVector angles)
 993 {
 994    DumpStart ("revolutionQuads", tstart << center << axis << angles);
 995
 996    Real3     axe;
 997    Elements* grid = new Elements (this);
 998
 999    grid->checkQuads  (tstart);
1000    grid->checkOrig   (2, center);
1001    grid->checkAxis   (axis,   axe);
1002    grid->checkVector (4, angles);
1003
1004    if (grid->isValid())
1005       {
1006       grid->extrudeQuads  (tstart, axe, angles, true, center);
1007       }
1008
1009    DumpReturn (grid);
1010    return grid;
1011 }
1012 // ------------------------------------------------------------------------
1013 // ------------------------------------------------------------------------
1014 // ------------------------------------------------------------------------
1015 // ========================================================= joinQuadUni
1016 Elements* Document::joinQuadUni (Quad* quada, Quad* quadb, Vertex* va1,
1017                               Vertex* vb1, Vertex* va2, Vertex* vb2, int nbre)
1018 {
1019    DumpStart ("joinQuadUni", quada << quadb << va1<<vb1 << va2<<vb2 << nbre);
1020
1021    Quads      tquada;
1022    RealVector tlen;
1023    Elements*  grid = new Elements (this);
1024
1025    update ();
1026    grid->checkQuad  (quada);
1027    grid->checkQuad  (quadb, -1);
1028    grid->checkSense (3, va1, va2, quada);
1029    grid->checkSense (6, vb1, vb2, quadb);
1030
1031    nbre = nbre < 0 ? 0 : nbre;
1032    double dist = 1.0 / (nbre+1);
1033
1034    tquada.push_back (quada);
1035    for (int nro=0 ; nro<nbre; ++nro) tlen.push_back ((nro+1)*dist);
1036
1037    if (grid->isValid())
1038       {
1039       grid->joinQuads  (tquada, quadb, va1, vb1, va2, vb2, tlen);
1040       }
1041
1042    DumpReturn (grid);
1043    return grid;
1044 }
1045 // ========================================================= joinQuad
1046 Elements* Document::joinQuad (Quad* quada, Quad* quadb, Vertex* va1,
1047                     Vertex* vb1, Vertex* va2, Vertex* vb2, RealVector& tlen)
1048 {
1049    DumpStart ("joinQuad", quada << quadb << va1 << vb1 << va2 << vb2 << tlen);
1050
1051    Quads      tquada;
1052    Elements*  grid = new Elements (this);
1053
1054    update ();
1055    grid->checkQuad   (quada);
1056    grid->checkQuad   (quadb, -1);
1057    grid->checkSense  (3, va1, va2, quada);
1058    grid->checkSense  (6, vb1, vb2, quadb);
1059    grid->checkVector (7, tlen, 0, true);
1060
1061    tquada.push_back (quada);
1062
1063    if (grid->isValid())
1064       {
1065       grid->joinQuads  (tquada, quadb, va1, vb1, va2, vb2, tlen);
1066       }
1067
1068    DumpReturn (grid);
1069    return grid;
1070 }
1071 // ========================================================= joinQuadsUni
1072 Elements* Document::joinQuadsUni (Quads tquada, Quad* quadb, Vertex* va1,
1073                            Vertex* vb1, Vertex* va2, Vertex* vb2, int nbre)
1074 {
1075    DumpStart ("joinQuadsUni", tquada << quadb << va1 << vb1 << va2 << vb2
1076                                      << nbre);
1077    RealVector tlen;
1078    Elements* grid = new Elements (this);
1079
1080    update ();
1081    grid->checkQuads (tquada);
1082    grid->checkQuad  (quadb, -1);
1083    if (tquada.size()>0)
1084        grid->checkSense (3, va1, va2, tquada[0]);
1085    grid->checkSense (6, vb1, vb2, quadb);
1086
1087    nbre = nbre < 0 ? 0 : nbre;
1088    double dist = 1.0 / (nbre+1);
1089
1090    for (int nro=0 ; nro<nbre; ++nro) tlen.push_back ((nro+1)*dist);
1091
1092    if (grid->isValid())
1093       {
1094       grid->joinQuads  (tquada, quadb, va1, vb1, va2, vb2, tlen);
1095       }
1096
1097    DumpReturn (grid);
1098    return grid;
1099 }
1100 // ========================================================= joinQuads
1101 Elements* Document::joinQuads (Quads tquada, Quad* quadb, Vertex* va1,
1102                       Vertex* vb1, Vertex* va2, Vertex* vb2, RealVector& tlen)
1103 {
1104    DumpStart ("joinQuads", tquada << quadb << va1 << vb1 << va2 << vb2 << tlen);
1105
1106    Elements* grid = new Elements (this);
1107
1108    update ();
1109    grid->checkQuads (tquada);
1110    grid->checkQuad  (quadb, -1);
1111    if (tquada.size()>0)
1112        grid->checkSense (3, va1, va2, tquada[0]);
1113    grid->checkSense  (6, vb1, vb2, quadb);
1114    grid->checkVector (7, tlen, 1, true);
1115
1116    if (grid->isValid())
1117       {
1118       grid->joinQuads  (tquada, quadb, va1, vb1, va2, vb2, tlen);
1119       }
1120
1121    DumpReturn (grid);
1122    return grid;
1123 }
1124 // ========================================================= cutUni
1125 Elements* Document::cutUni (Edge* edge, int nbre)
1126 {
1127    DumpStart ("cutUni", edge << nbre);
1128
1129    Elements*  grid = new Elements (this);
1130
1131    if (BadElement (edge))
1132       {
1133       grid->setError (HERR);
1134       Mess << "Argument number 1 : edge is not valid";
1135       }
1136
1137    if (nbre <= 0)
1138       {
1139       grid->setError (HERR);
1140       Mess << "Argument number 2 : number od subdivisions must be positive";
1141       nbre = 0;
1142       }
1143
1144    double dist = 1.0 / (nbre+1);
1145
1146    RealVector tlen;
1147    for (int nro=0 ; nro<nbre; ++nro) tlen.push_back ((nro+1)*dist);
1148
1149    if (grid->isValid())
1150       {
1151       grid->cutEdge  (edge, tlen);
1152       }
1153
1154    DumpReturn (grid);
1155    return grid;
1156 }
1157 // ========================================================= cut
1158 Elements* Document::cut (Edge* edge, RealVector& tlen)
1159 {
1160    DumpStart ("cut", edge << tlen);
1161
1162    Elements* grid = new Elements (this);
1163
1164    if (BadElement (edge))
1165       {
1166       grid->setError (HERR);
1167       Mess << "Argument number 1 : edge is not valid";
1168       }
1169
1170    grid->checkVector (2, tlen, 1, true);
1171
1172    if (grid->isValid())
1173       {
1174       grid->cutEdge  (edge, tlen);
1175       }
1176
1177    DumpReturn (grid);
1178    return grid;
1179 }
1180 END_NAMESPACE_HEXA