src/HEXABLOCK/HexDocument_v6.cxx

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