src/HEXABLOCK/HexDocument_quads.cxx

   1
   2 // C++ : Classe Document : Methodes internes 2011
   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/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include "HexDocument.hxx"
  24
  25 #include <cmath>
  26 #include <map>
  27
  28 #include "HexVertex.hxx"
  29 #include "HexEdge.hxx"
  30 #include "HexQuad.hxx"
  31 #include "HexHexa.hxx"
  32
  33 #include "HexLaw.hxx"
  34
  35 #include "HexAnaQuads.hxx"
  36 #include "HexElements.hxx"
  37 #include "HexCramer.hxx"
  38 #include "HexGlobale.hxx"
  39
  40 BEGIN_NAMESPACE_HEXA
  41
  42 #define PermuterEdges(e1,e2) permuter_edges (e1, e2, #e1, #e2)
  43 void    permuter_edges  (Edge* &e1, Edge* &e2, cpchar n1=NULL, cpchar n2=NULL);
  44 double* prod_vectoriel (Edge* e1, Edge* e2, double result[]);
  45
  46 static bool db = false;
  47
  48 // ======================================================== copyDocument
  49 Document* Document::copyDocument ()
  50 {
  51    string nom = "CopyOf_";
  52    nom += el_name;
  53
  54    Document* clone = new Document (nom.c_str());
  55
  56    for (EltBase* elt = doc_first_elt[EL_VERTEX]->next (); elt!=NULL;
  57                  elt = elt->next())
  58        {
  59        if (elt !=NULL && elt->isHere())
  60           {
  61           Vertex* node = static_cast <Vertex*> (elt);
  62           node->duplicate (clone);
  63           }
  64        }
  65
  66    for (int type=EL_EDGE ; type <= EL_HEXA ; type++)
  67        {
  68        for (EltBase* elt = doc_first_elt[type]->next (); elt!=NULL;
  69                      elt = elt->next())
  70            {
  71            if (elt !=NULL && elt->isHere())
  72               elt->duplicate ();
  73            }
  74        }
  75
  76    for (int nro=0 ; nro<nbr_laws ; nro++)
  77        {
  78        Law* law = new Law (doc_laws [nro]);
  79        clone->doc_laws.push_back (law);
  80        }
  81
  82    return clone;
  83 }
  84 // ---------------------------------------------------------------
  85 // ============================================================== addHexa2quads
  86 Hexa* Document::addHexa2Quads (Quad* q1, Quad* q2)
  87 {
  88    DumpStart ("addHexa2Quads", q1 << q2);
  89    AnaQuads ana_quads (q1, q2);
  90
  91    Hexa* hexa = NULL;
  92    if (ana_quads.status != HOK)
  93       hexa = NULL;
  94
  95    else if (ana_quads.nbr_aretes == 0)
  96       hexa = addHexaQuadsAB (ana_quads);
  97
  98    else if (ana_quads.nbr_aretes == 1)
  99       hexa = addHexaQuadsAC (ana_quads);
 100
 101    DumpReturn (hexa);
 102    return hexa;
 103 }
 104 // ============================================================= addHexa3quads
 105 Hexa* Document::addHexa3Quads (Quad* q1, Quad* q2, Quad* q3)
 106 {
 107    DumpStart ("addHexa3Quads", q1 << q2 << q3);
 108    AnaQuads ana_quads (q1, q2, q3);
 109
 110    Hexa* hexa = NULL;
 111    if (ana_quads.status != HOK)
 112       hexa = NULL;
 113
 114    else if (ana_quads.nbr_aretes == 2)
 115       hexa = addHexaQuadsACD (ana_quads);
 116
 117    else if (ana_quads.nbr_aretes == 3)
 118       hexa = addHexaQuadsACE (ana_quads);
 119
 120    DumpReturn (hexa);
 121    return hexa;
 122 }
 123 // ============================================================= addHexa4quads
 124 Hexa* Document::addHexa4Quads (Quad* q1, Quad* q2, Quad* q3, Quad* q4)
 125 {
 126    DumpStart ("addHexa4Quads", q1 << q2 << q3 << q4);
 127    AnaQuads ana_quads (q1, q2, q3, q4);
 128
 129    Hexa* hexa = NULL;
 130    if (ana_quads.status != HOK)
 131       hexa = NULL;
 132
 133    else if (ana_quads.nbr_aretes == 4)
 134       hexa = addHexaQuadsABCD (ana_quads);
 135
 136    else if (ana_quads.nbr_aretes == 5)
 137       hexa = addHexaQuadsACDE (ana_quads);
 138
 139    DumpReturn (hexa);
 140    return hexa;
 141 }
 142 // ============================================================== addHexa5quads
 143 Hexa* Document::addHexa5Quads (Quad* q1, Quad* q2, Quad* q3, Quad* q4, Quad* q5)
 144 {
 145    DumpStart ("addHexa5Quads", q1 << q2 << q3 << q4 << q5);
 146    AnaQuads ana_quads (q1, q2, q3, q4, q5);
 147    if (ana_quads.status != HOK)
 148       return NULL;
 149    else if (ana_quads.nbr_aretes != 8)
 150       return NULL;
 151
 152    int qbase = NOTHING;
 153    for (int nquad=0 ; nquad < ana_quads.nbr_quads ; nquad++)
 154        if (ana_quads.inter_nbre [nquad] == 4)
 155           qbase = nquad;
 156
 157    if (qbase == NOTHING)
 158       return NULL;
 159
 160    Edge* tedge [QUAD4];
 161    Quad* tquad [QUAD4];
 162    for (int nedge=0 ; nedge < QUAD4 ; nedge++)
 163        {
 164        int nq    = ana_quads.inter_quad [qbase] [nedge];
 165        int ned1  = ana_quads.inter_edge [nq]    [qbase];
 166        int ned2  = (ned1 + 2) MODULO QUAD4;
 167        Quad* mur = ana_quads.tab_quads[nq];
 168        tedge [nedge] = mur->getEdge (ned2);
 169        tquad [nedge] = mur;
 170        }
 171
 172    Quad*  q_a  = ana_quads.tab_quads[qbase];
 173    Quad*  q_b  = new Quad (tedge[0], tedge[1], tedge[2], tedge[3]);
 174    Hexa*  hexa = new Hexa (q_a, q_b, tquad[0], tquad[2], tquad[1], tquad[3]);
 175
 176    DumpReturn (hexa);
 177    return hexa;
 178 }
 179 // ---------------------------------------------------------------
 180 // ========================================================== addHexaquadsAB
 181 Hexa* Document::addHexaQuadsAB (AnaQuads& strquads)
 182 {
 183    Quad* q_a = strquads.tab_quads[0];
 184    Quad* q_b = strquads.tab_quads[1];
 185
 186    double dmin  = 0;
 187    int    sens  = 1;
 188    int    decal = 0;
 189    for (int is = 0 ; is<2 ; is++)
 190        {
 191        int ns = 1-2*is;
 192        for (int ndec = 0 ; ndec<QUAD4 ; ndec++)
 193            {
 194            double dist = 0;
 195            for (int na = 0 ; na<QUAD4 ; na++)
 196                {
 197                int nb = (ndec + QUAD4 + ns*na) MODULO QUAD4;
 198                dist += distance (q_a->getVertex (na),
 199                                  q_b->getVertex (nb));
 200                }
 201            if (ndec==0 && is==0)
 202               {
 203               decal = ndec;
 204               dmin  = dist;
 205               sens  = ns;
 206               }
 207            else if (dist<dmin)
 208               {
 209               dmin  = dist;
 210               decal = ndec;
 211               sens  = ns;
 212               }
 213            }
 214        }
 215
 216    Edge* tedge [QUAD4];
 217    Quad* tquad [QUAD4];
 218    for (int na = 0 ; na<QUAD4 ; na++)
 219        {
 220        int nb = (decal + QUAD4 + sens*na) MODULO QUAD4;
 221        tedge [na] = new Edge (q_a->getVertex (na), q_b->getVertex (nb));
 222        }
 223
 224    for (int nal = 0 ; nal<QUAD4 ; nal++)
 225        {
 226        int nar = (nal+1) MODULO QUAD4;
 227        Edge* e_left  = tedge [nal];
 228        Edge* e_right = tedge [nar];
 229        Edge* e_ax  = q_a->findEdge (e_left ->getVertex (V_AMONT),
 230                                     e_right->getVertex (V_AMONT));
 231        Edge* e_bx  = q_b->findEdge (e_left ->getVertex (V_AVAL),
 232                                     e_right->getVertex (V_AVAL));
 233        tquad [nal] = new Quad (e_ax, tedge [nal], e_bx, tedge [nar]);
 234        }
 235
 236    Hexa* hexa = new Hexa (q_a, q_b, tquad[0], tquad[2], tquad[1], tquad[3]);
 237    return hexa;
 238 }
 239 // ========================================================== addHexaquadsAC
 240 Hexa* Document::addHexaQuadsAC (AnaQuads& strquads)
 241 {
 242    Quad* q_a = strquads.tab_quads[0];
 243    Quad* q_c = strquads.tab_quads[1];
 244
 245    Vertex* tv_bdx [V_TWO];   // x = e ou f
 246    Edge*   te_dX  [V_TWO];
 247    Edge*   te_bX  [V_TWO];
 248    Quad*   tq_ef  [V_TWO];
 249
 250    int neda0  = strquads.inter_edge [0] [1];
 251
 252    Edge* e_ac = q_a->getEdge (neda0);
 253
 254    for (int ns=V_AMONT; ns<=V_AVAL ; ns++)
 255        {
 256        Vertex* vx1 = e_ac->getVertex (ns);
 257        Vertex* vx2 = e_ac->getVertex (1-ns);
 258
 259        int nda2 = q_a->indexVertex (vx2);
 260        nda2 = (nda2 +2) MODULO QUAD4;
 261
 262        int ndc2 = q_c->indexVertex (vx2);
 263        ndc2 = (ndc2 +2) MODULO QUAD4;
 264
 265        Vertex* vxa = q_a->getVertex (nda2);
 266        Vertex* vxc = q_c->getVertex (ndc2);
 267
 268        double dx  = (vxa->getX() - vx1->getX()) + (vxc->getX() - vx1->getX());
 269        double dy  = (vxa->getY() - vx1->getY()) + (vxc->getY() - vx1->getY());
 270        double dz  = (vxa->getZ() - vx1->getZ()) + (vxc->getZ() - vx1->getZ());
 271        tv_bdx [ns] = new Vertex (this, vx1->getX()+dx, vx1->getY()+dy,
 272                                                       vx1->getZ()+dz);
 273        Edge* edga = q_a->findEdge (vx1, vxa);
 274        Edge* edgc = q_c->findEdge (vx1, vxc);
 275
 276        te_dX [ns] = new Edge (vxa, tv_bdx[ns]);
 277        te_bX [ns] = new Edge (vxc, tv_bdx[ns]);
 278        tq_ef [ns] = new Quad (edga, te_dX[ns], te_bX[ns], edgc);
 279        }
 280
 281    int ff = 0;
 282    Edge* e_bd = new Edge (tv_bdx[V_AMONT], tv_bdx[V_AVAL]);
 283    Edge* e_ad = q_a->getOpposEdge (e_ac, ff);
 284    Edge* e_bc = q_c->getOpposEdge (e_ac, ff);
 285
 286    Quad* q_d = new Quad (e_bd, te_dX[V_AMONT], e_ad, te_dX[V_AVAL]);
 287    Quad* q_b = new Quad (e_bd, te_bX[V_AMONT], e_bc, te_bX[V_AVAL]);
 288
 289    Hexa*  hexa  = new Hexa (q_a, q_b, q_c, q_d, tq_ef[V_AMONT], tq_ef[V_AVAL]);
 290    return hexa;
 291 }
 292 // ========================================================= addHexaquadsACE
 293 // ==== Construction d'un hexaedre a partir d'un triedre
 294 /*
 295        6=bed  +----bd-----+ bdf=7
 296              /|          /|
 297            be |   B    bf |
 298            /  |        /  |
 299     4=bce +----bc-----+...|...bcf=5
 300           |  de     D |   df
 301           | E |       | F |             z
 302          ce   | C     cf  |             ^
 303   2=ade...|...+----ad-|---+ adf=3       |   y
 304           |  /        |  /              |  /
 305           | ae    A   | af              | /
 306           |/          |/                |/
 307     0=ace +----ac-----+ acf=1           +----->  x
 308
 309 On connait les faces A, C, E
 310 Il faut determiner le point bdf, intersection des faces B, D, F
 311
 312 bdf in B : Scalaire ((bdf-bce), Vectoriel (be, bc)) = 0
 313 bdf in D : Scalaire ((bdf-ade), Vectoriel (ad, de)) = 0
 314 bdf in F : Scalaire ((bdf-acf), Vectoriel (af, cf)) = 0
 315
 316 Soit 3 inconnues (Xbdf, Ybdf, Zbdf) et 3 equations
 317 Un merci a Francois Mougery qui m'a rappele quelques notions de geometrie 3D
 318
 319 Le systeme s'ecrit :
 320
 321   Scalaire ((M-bce), norm_b) = 0
 322   Scalaire ((M-ade), morm_d) = 0
 323   Scalaire ((M-acf), morm_f) = 0
 324
 325 <=>
 326
 327   Scalaire (M, norm_b) = Scalaire (bce, norm_b)
 328   scalaire (M, norm_d) = Scalaire (ade, norm_b)
 329   scalaire (M, norm_f) = Scalaire (acf, norm_b)
 330
 331 <=>
 332    norme_b.x*X + norme_b.y*Y + norme_b.z*Z = K1
 333    norme_d.x*X + norme_d.y*Y + norme_d.z*Z = K2
 334    norme_f.x*X + norme_f.y*Y + norme_f.z*Z = K3
 335
 336  * ----------------------------------------------------- */
 337 Hexa* Document::addHexaQuadsACE (AnaQuads& strquads)
 338 {
 339    Real3 v_bce,  v_ade,  v_acf, v_bdf;    // Sommets
 340    Real3 norm_b, norm_d, norm_f;   // Normales aux faces
 341    int   ff = 0;
 342
 343    Quad* q_a = strquads.tab_quads[0];
 344    Quad* q_c = strquads.tab_quads[1];
 345    Quad* q_e = strquads.tab_quads[2];
 346
 347    Edge* e_ac = q_a->commonEdge (q_c);
 348    Edge* e_ae = q_a->commonEdge (q_e);
 349    Edge* e_ce = q_c->commonEdge (q_e);
 350
 351    Edge* e_ad = q_a->getOpposEdge (e_ac, ff);
 352    Edge* e_af = q_a->getOpposEdge (e_ae, ff);
 353    Edge* e_cf = q_c->getOpposEdge (e_ce, ff);
 354
 355    Edge* e_bc = q_c->getOpposEdge (e_ac, ff);
 356    Edge* e_be = q_e->getOpposEdge (e_ae, ff);
 357    Edge* e_de = q_e->getOpposEdge (e_ce, ff);
 358
 359    e_ce->commonPoint (e_bc, v_bce);
 360    e_ae->commonPoint (e_ad, v_ade);
 361    e_ac->commonPoint (e_af, v_acf);
 362
 363    prod_vectoriel (e_be, e_bc, norm_b);  // Calcul normale a la face B
 364    prod_vectoriel (e_ad, e_de, norm_d);  // Calcul normale a la face D
 365    prod_vectoriel (e_af, e_cf, norm_f);  // Calcul normale a la face F
 366
 367    Real3 membre2 = { prod_scalaire (v_bce, norm_b),
 368                      prod_scalaire (v_ade, norm_d),
 369                      prod_scalaire (v_acf, norm_f) };
 370
 371    double matrix [] = { norm_b[dir_x],  norm_b[dir_y],  norm_b[dir_z],
 372                         norm_d[dir_x],  norm_d[dir_y],  norm_d[dir_z],
 373                         norm_f[dir_x],  norm_f[dir_y],  norm_f[dir_z] };
 374    Cramer systeme (DIM3);
 375    int ier = systeme.resoudre (matrix, membre2, v_bdf);
 376    if (ier != HOK)
 377       {
 378       printf (" addHexaQuadsACE : systeme impossible\n");
 379       return NULL;
 380       }
 381
 382    Vertex* s_bdf = new Vertex (this, v_bdf[dir_x], v_bdf[dir_y], v_bdf[dir_z]);
 383    Vertex* s_bde = e_be -> commonVertex (e_de);
 384    Vertex* s_bcf = e_bc -> commonVertex (e_cf);
 385    Vertex* s_adf = e_af -> commonVertex (e_ad);
 386
 387    Edge* e_bd = new Edge (s_bdf, s_bde);
 388    Edge* e_bf = new Edge (s_bdf, s_bcf);
 389    Edge* e_df = new Edge (s_bdf, s_adf);
 390
 391    Quad* q_b = new Quad (e_bc, e_be, e_bd, e_bf);
 392    Quad* q_d = new Quad (e_de, e_ad, e_df, e_bd);
 393    Quad* q_f = new Quad (e_af, e_cf, e_bf, e_df);
 394
 395    Hexa*  hexa = new Hexa (q_a, q_b, q_c, q_d, q_e, q_f);
 396    return hexa;
 397 }
 398 // ========================================================= addHexaquadsACD
 399 // ==== Construction d'un hexaedre a partir d'un U
 400 Hexa* Document::addHexaQuadsACD (AnaQuads& strquads)
 401 {
 402    int pos_a = NOTHING;
 403    for (int np=0 ; np<3 && pos_a==NOTHING ; np++)
 404        if (strquads.inter_nbre[np]==2)
 405           pos_a = np;
 406
 407    if (pos_a==NOTHING)
 408       return NULL;
 409
 410    int pos_c = (pos_a+1) MODULO 3;
 411    int pos_d = (pos_a+2) MODULO 3;
 412    Quad* q_a = strquads.tab_quads[pos_a];
 413    Quad* q_c = strquads.tab_quads[pos_c];
 414    Quad* q_d = strquads.tab_quads[pos_d];
 415
 416    int   na_ac = strquads.inter_edge[pos_a][pos_c]; // Nro dans  q_a de e_ac
 417    int   nc_ac = strquads.inter_edge[pos_c][pos_a]; // Nro dans  q_c de e_ac
 418    int   nd_ad = strquads.inter_edge[pos_d][pos_a]; // Nro dans  q_d de e_ad
 419
 420    Edge* e_ae  = q_a->getEdge ((na_ac + 1) MODULO QUAD4); // Arbitraire
 421    Edge* e_af  = q_a->getEdge ((na_ac + 3) MODULO QUAD4); // Arbitraire
 422
 423    Edge* e_bc  = q_c->getEdge ((nc_ac + 2) MODULO QUAD4);
 424    Edge* e_bd  = q_d->getEdge ((nd_ad + 2) MODULO QUAD4);
 425
 426    Edge* e_ce  = q_c->getEdge ((nc_ac + 1) MODULO QUAD4);
 427    Edge* e_cf  = q_c->getEdge ((nc_ac + 3) MODULO QUAD4);
 428
 429    Edge* e_de  = q_d->getEdge ((nd_ad + 1) MODULO QUAD4);
 430    Edge* e_df  = q_d->getEdge ((nd_ad + 3) MODULO QUAD4);
 431
 432    Vertex* v_ace = e_ae->commonVertex (e_ce);
 433    Vertex* v_ade = e_ae->commonVertex (e_de);
 434    if (v_ace==NULL)
 435       {
 436       permuter_edges (e_ce, e_cf);
 437       v_ace = e_ae->commonVertex (e_ce);
 438       }
 439
 440    if (v_ade==NULL)
 441       {
 442       permuter_edges (e_de, e_df);
 443       v_ade = e_ae->commonVertex (e_de);
 444       }
 445
 446    Vertex* v_acf = e_af->commonVertex (e_cf);
 447    Vertex* v_adf = e_af->commonVertex (e_df);
 448
 449    Vertex* v_bce = e_ce->opposedVertex (v_ace);
 450    Vertex* v_bde = e_de->opposedVertex (v_ade);
 451    Vertex* v_bcf = e_cf->opposedVertex (v_acf);
 452    Vertex* v_bdf = e_df->opposedVertex (v_adf);
 453
 454    Edge*  e_be = new Edge (v_bce, v_bde);
 455    Edge*  e_bf = new Edge (v_bcf, v_bdf);
 456
 457    Quad* q_b = new Quad  (e_be, e_bc, e_bf, e_bd);
 458    Quad* q_e = new Quad  (e_ae, e_ce, e_be, e_de);
 459    Quad* q_f = new Quad  (e_af, e_cf, e_bf, e_df);
 460
 461    Hexa*  hexa = new Hexa (q_a, q_b, q_c, q_d, q_e, q_f);
 462    return hexa;
 463 }
 464 // ========================================================= addHexaquadsABCD
 465 Hexa* Document::addHexaQuadsABCD (AnaQuads& strquads)
 466 {
 467    int pos_a = 0;
 468    int pos_b = NOTHING;
 469    int pos_c = NOTHING;
 470    int pos_d = NOTHING;
 471
 472    for (int np=1 ; np<4 ; np++)
 473        {
 474        if (strquads.inter_edge [pos_a] [np] == NOTHING )
 475           pos_b = np;
 476        else if (pos_c==NOTHING)
 477           pos_c = np;
 478        else
 479           pos_d = np;
 480        }
 481
 482    if (pos_b==NOTHING || pos_c==NOTHING || pos_d==NOTHING)
 483       return NULL;
 484
 485    Quad* q_a = strquads.tab_quads [pos_a];
 486    Quad* q_b = strquads.tab_quads [pos_b];
 487    Quad* q_c = strquads.tab_quads [pos_c];
 488    Quad* q_d = strquads.tab_quads [pos_d];
 489
 490    int   na_ac = strquads.inter_edge[pos_a][pos_c]; // Nro dans  q_a de e_ac
 491    int   nc_ac = strquads.inter_edge[pos_c][pos_a]; // Nro dans  q_c de e_ac
 492    int   nd_ad = strquads.inter_edge[pos_d][pos_a]; // Nro dans  q_d de e_ad
 493    int   nb_bc = strquads.inter_edge[pos_b][pos_c]; // Nro dans  q_b de e_bc
 494
 495    Edge* e_ae  = q_a->getEdge ((na_ac + 1) MODULO QUAD4);
 496    Edge* e_af  = q_a->getEdge ((na_ac + 3) MODULO QUAD4);
 497
 498    Edge* e_ce  = q_c->getEdge ((nc_ac + 1) MODULO QUAD4);
 499    Edge* e_cf  = q_c->getEdge ((nc_ac + 3) MODULO QUAD4);
 500
 501    Edge* e_de  = q_d->getEdge ((nd_ad + 1) MODULO QUAD4);
 502    Edge* e_df  = q_d->getEdge ((nd_ad + 3) MODULO QUAD4);
 503
 504    Edge*  e_be = q_b->getEdge ((nb_bc + 1) MODULO QUAD4);
 505    Edge*  e_bf = q_b->getEdge ((nb_bc + 3) MODULO QUAD4);
 506
 507    if (db)
 508       {
 509       HexDump (q_a);
 510       HexDump (q_b);
 511       HexDump (q_c);
 512       HexDump (q_d);
 513
 514       HexDump (e_ae);
 515       HexDump (e_af);
 516       HexDump (e_ce);
 517       HexDump (e_cf);
 518       HexDump (e_de);
 519       HexDump (e_df);
 520       HexDump (e_be);
 521       HexDump (e_bf);
 522       }
 523
 524    if (e_ae->commonVertex (e_ce) == NULL)
 525       PermuterEdges (e_ce, e_cf);
 526
 527    if (e_ae->commonVertex (e_de) == NULL)
 528       PermuterEdges (e_de, e_df);
 529
 530    if (e_ce->commonVertex (e_be) == NULL)
 531       PermuterEdges (e_be, e_bf);
 532
 533    Quad* q_e = new Quad  (e_ae, e_ce, e_be, e_de);
 534    Quad* q_f = new Quad  (e_af, e_cf, e_bf, e_df);
 535
 536    Hexa*  hexa = new Hexa (q_a, q_b, q_c, q_d, q_e, q_f);
 537    return hexa;
 538 }
 539 // ========================================================= addHexaquadsACDE
 540 Hexa* Document::addHexaQuadsACDE (AnaQuads& strquads)
 541 {
 542    int pos_a = NOTHING;
 543    int pos_c = NOTHING;
 544    int pos_d = NOTHING;
 545    int pos_e = NOTHING;
 546
 547    for (int np=0 ; np<4  ; np++)
 548        if (strquads.inter_nbre[np]==3)
 549           {
 550           if (pos_a == NOTHING) pos_a = np;
 551              else               pos_e = np;
 552           }
 553        else if (strquads.inter_nbre[np]==2)
 554           {
 555           if (pos_c == NOTHING) pos_c = np;
 556              else               pos_d = np;
 557           }
 558
 559    if (pos_a==NOTHING || pos_c==NOTHING  || pos_d==NOTHING || pos_e==NOTHING)
 560       return NULL;
 561
 562    Quad* q_a = strquads.tab_quads[pos_a];
 563    Quad* q_c = strquads.tab_quads[pos_c];
 564    Quad* q_d = strquads.tab_quads[pos_d];
 565    Quad* q_e = strquads.tab_quads[pos_e];
 566
 567    int   na_ac = strquads.inter_edge[pos_a][pos_c]; // Nro dans  q_a de e_ac
 568
 569    int   nc_ac = strquads.inter_edge[pos_c][pos_a]; // Nro dans  q_c de e_ac
 570    int   nc_ce = strquads.inter_edge[pos_c][pos_e]; // Nro dans  q_c de e_ce
 571
 572    int   nd_ad = strquads.inter_edge[pos_d][pos_a]; // Nro dans  q_d de e_ad
 573    int   nd_de = strquads.inter_edge[pos_d][pos_e]; // Nro dans  q_d de e_de
 574
 575    int   ne_ae = strquads.inter_edge[pos_e][pos_a]; // Nro dans  q_e de e_ae
 576
 577    Edge* e_af  = q_a->getEdge ((na_ac + 3) MODULO QUAD4);
 578    Edge* e_bc  = q_c->getEdge ((nc_ac + 2) MODULO QUAD4);
 579    Edge* e_cf  = q_c->getEdge ((nc_ce + 2) MODULO QUAD4);
 580    Edge* e_bd  = q_d->getEdge ((nd_ad + 2) MODULO QUAD4);
 581    Edge* e_df  = q_d->getEdge ((nd_de + 2) MODULO QUAD4);
 582    Edge* e_be  = q_e->getEdge ((ne_ae + 2) MODULO QUAD4);
 583
 584    Vertex* v_bcf = e_cf->opposedVertex (e_cf->commonVertex (e_af));
 585    Vertex* v_bdf = e_df->opposedVertex (e_df->commonVertex (e_af));
 586
 587    Edge*   e_bf = new Edge (v_bcf, v_bdf);
 588    Quad*   q_b  = new Quad (e_be, e_bc, e_bf, e_bd);
 589    Quad*   q_f  = new Quad (e_af, e_cf, e_bf, e_df);
 590
 591    Hexa*  hexa = new Hexa (q_a, q_b, q_c, q_d, q_e, q_f);
 592    return hexa;
 593 }
 594 // ========================================================= replace
 595 Elements* Document::replace (Quads pattern, Vertex* p1, Vertex* c1,
 596                              Vertex* p2, Vertex* c2, Vertex* p3, Vertex* c3)
 597 {
 598    DumpStart ("replace", pattern << p1 << c1 << p2 << c2 << p3 << c3);
 599
 600    Elements* t_hexas = new Elements (this);
 601    int ier = t_hexas->replaceHexas (pattern, p1, c1, p2, c2, p3, c3);
 602    if (ier!=HOK)
 603       {
 604       Mess << " **** Error in Document::replace" ;
 605       t_hexas->setError (ier);
 606       }
 607
 608    DumpReturn (t_hexas);
 609    return      t_hexas;
 610 }
 611 // ========================================================= print_replace
 612 void print_replace (Edge* zig, Edge*  zag)
 613 {
 614    cout << zig->getName() << " = (" << zig->getVertex(0)->getName()
 615         << ", " << zig->getVertex(1)->getName() << ") est clone en ";
 616
 617    cout << zag->getName() << " = (" << zag->getVertex(0)->getName()
 618         << ", " << zag->getVertex(1)->getName() << ")" << endl;
 619 }
 620 // ========================================================= only_in_hexas
 621 bool only_in_hexas (Hexas& thexas, Quad* quad)
 622 {
 623    int nbhexas = thexas.size();
 624    int nbp     = quad->getNbrParents();
 625    for (int nh=0 ; nh <nbp ; nh++)
 626        {
 627        bool pasla = true;
 628        Hexa* hexa = quad->getParent (nh);
 629        for (int nc=0 ; pasla && nc < nbhexas ; nc++)
 630            pasla = hexa != thexas [nc];
 631        if (pasla)
 632            return false;
 633        }
 634    return true;
 635 }
 636 // ========================================================= only_in_hexas
 637 bool only_in_hexas (Hexas& thexas, Edge*  edge)
 638 {
 639    return false;
 640    int nbp = edge->getNbrParents();
 641    for (int nq=0 ; nq <nbp ; nq++)
 642        {
 643        Quad* quad = edge->getParent   (nq);
 644        if (NOT only_in_hexas (thexas, quad))
 645           {
 646           cout << " ... inMoreHexas " << edge->makeDefinition() << endl;
 647           return false;
 648           }
 649        }
 650    cout << " ... only_in_hexas " << edge->makeDefinition() << endl;
 651    return true;
 652 }
 653 // ========================================================= replace_vertex
 654 void replace_vertex (Hexas& thexas, Vertex* node,  Vertex* par)
 655 {
 656    int nbh = thexas.size();
 657    for (int nh=0 ; nh <nbh ; nh++)
 658        thexas[nh]->replaceVertex (node, par);
 659 }
 660 // ========================================================= disconnectEdges
 661 Elements* Document::disconnectEdges (Hexas thexas, Edges  tedges)
 662 {
 663    DumpStart ("disconnectEdges",  thexas << tedges);
 664
 665    if (db)
 666       cout << " +++ Disconnect Edges" << endl;
 667
 668    Elements* grid  = new Elements (this);
 669
 670    grid->checkDisco (thexas, tedges);
 671
 672    int nbedges = tedges.size();
 673    int nbhexas = thexas.size();
 674
 675    if (nbhexas != nbedges)
 676       {
 677       cout << " **** Error in Document::disconnectEdges\n" << endl;
 678       cout << " **** Number of Edges and number of Hexas are different\n"
 679            << endl;
 680       return NULL;
 681       }
 682    else if (nbhexas==1)
 683       {
 684       update ();
 685       thexas[0]->disconnectEdge (tedges[0], grid);
 686       return    grid;
 687       }
 688
 689    for (int nro=0 ; nro<nbedges ; nro++)
 690        {
 691        if (BadElement (tedges[nro]))
 692           {
 693           cout << " **** Eddge number " << nro+1 << " is incorrect"
 694                << endl;
 695           return NULL;
 696           }
 697        if (BadElement (thexas[nro]))
 698           {
 699           cout << " **** Hexa number " << nro+1 << " is incorrect"
 700                << endl;
 701           return NULL;
 702           }
 703        if (db)
 704           cout << nro+1 << " hexa = " << thexas[nro]->getName ()
 705                         << ", edge = " << tedges[nro]->getName ()
 706                         << " = (" << tedges[nro]->getVertex(0)->getName ()
 707                         << ", "   << tedges[nro]->getVertex(1)->getName ()
 708                         << ")" << endl;
 709        }
 710
 711    for (int nro=0 ; nro<nbhexas ; nro++)
 712        {
 713        int ned = thexas[nro]->findEdge (tedges[nro]);
 714        if (ned==NOTHING)
 715           {
 716           cout << " **** Edge number " << nro+1
 717                << " doesnt belong to correspondant hexa" << endl;
 718           return NULL;
 719           }
 720        }
 721
 722    vector <Vertex*> tvertex (nbedges+1);
 723
 724    for (int nro=1 ; nro<nbedges ; nro++)
 725        {
 726        tvertex[nro] = tedges[nro]->commonVertex (tedges[nro-1]);
 727        if (tvertex[nro]==NULL)
 728           {
 729           cout << " **** Edge number " << nro
 730                << " doesnt intesect next edge" << endl;
 731           return NULL;
 732           }
 733        }
 734
 735    int nv0 = tedges[0]        ->inter (tedges[1]);
 736    int nvn = tedges[nbedges-1]->inter (tedges[nbedges-2]);
 737    tvertex [0]       = tedges[0]        ->getVertex (1-nv0);
 738    tvertex [nbedges] = tedges[nbedges-1]->getVertex (1-nvn);
 739
 740    for (int nro=0 ; nro<nbhexas ; nro++)
 741        {
 742        int ned = thexas[nro]->findEdge (tedges[nro]);
 743        if (ned==NOTHING)
 744           {
 745           cout << " **** Edge number " << nro+1
 746                << " doesnt belong to correspondant hexa" << endl;
 747           return NULL;
 748           }
 749        }
 750                       // Fin des controles, on peut y aller ...
 751
 752    map <Edge*, int> state_edge;
 753    map <Quad*, int> state_quad;
 754    enum { UNDEFINED, REPLACED, AS_IS };
 755
 756    map <Vertex*, Vertex*> new_vertex;
 757    map <Edge*,   Edge*>   new_edge;
 758    map <Quad*,   Quad*>   new_quad;
 759
 760    map <Vertex*, Vertex*> :: iterator it_vertex;
 761    map <Edge*,   Edge*>   :: iterator it_edge;
 762    map <Quad*,   Quad*>   :: iterator it_quad;
 763
 764 #define VertexIsNew(v) (it_vertex=new_vertex.find(v))!=new_vertex.end()
 765
 766    Vertex*   node1     = NULL;
 767
 768    for (int nro=0 ; nro<=nbedges ; nro++)
 769        {
 770        Vertex* node0 = node1;
 771        node1 = new Vertex (tvertex[nro]);
 772        grid->addVertex  (node1);
 773        new_vertex [tvertex[nro]] = node1;
 774        if (db)
 775           {
 776           cout << nro << " : "         << tvertex[nro]->getName()
 777                << " est clone en " << node1->getName() << endl;
 778           }
 779
 780        if (nro>0)
 781           {
 782           Edge* edge = new Edge (node0, node1);
 783           grid->addEdge  (edge);
 784           new_edge   [tedges[nro-1]] = edge;
 785           state_edge [tedges[nro-1]] = REPLACED;
 786           if (db)
 787              print_replace (tedges[nro-1], edge);
 788           }
 789        }
 790
 791    if (db)
 792       cout << "_____________________________ Autres substitutions" << endl;
 793
 794    // Un edge non remplace, qui contient un vertex remplace
 795    //         commun a plus de 2 faces (donc appartenant a un autre hexa)
 796    //         doit etre duplique
 797
 798    for (int nro=0 ; nro<nbhexas ; nro++)
 799        {
 800        Hexa* hexa = thexas [nro];
 801        for (int nro=0 ; nro<HE_MAXI ; nro++)
 802            {
 803            Edge* edge = hexa->getEdge (nro);
 804            if (state_edge[edge]==UNDEFINED)
 805               {
 806               Vertex* v1 = edge->getVertex (V_AMONT);
 807               Vertex* v2 = edge->getVertex (V_AVAL);
 808               int etat = REPLACED;
 809               if (VertexIsNew (v1))
 810                  {
 811                  if (only_in_hexas (thexas, edge))
 812                     {
 813                     replace_vertex (thexas, v1, new_vertex[v1]);
 814                     etat = AS_IS;
 815                     }
 816                  else
 817                     v1 = new_vertex [v1];
 818                  }
 819               else if (VertexIsNew (v2))
 820                  {
 821                  if (only_in_hexas (thexas, edge))
 822                     {
 823                     replace_vertex (thexas, v2, new_vertex[v2]);
 824                     etat = AS_IS;
 825                     }
 826                  else
 827                     v2 = new_vertex [v2];
 828                  }
 829               else
 830                  etat = AS_IS;
 831
 832               if (etat==REPLACED)
 833                  {
 834                  Edge* arete = new Edge (v1, v2);
 835                  new_edge   [edge] = arete;
 836                  grid->addEdge  (arete);
 837                  if (db)
 838                     print_replace (edge, arete);
 839                  }
 840               state_edge [edge] = etat;
 841               }
 842            }
 843        }
 844
 845    // Un quad non remplace, qui contient un edge remplace
 846    //         commun a plus de 2 Hexas
 847    //         doit etre duplique
 848
 849    for (int nro=0 ; nro<nbhexas ; nro++)
 850        {
 851        Hexa* hexa = thexas [nro];
 852        for (int nro=0 ; nro<HQ_MAXI ; nro++)
 853            {
 854            Quad* quad = hexa->getQuad (nro);
 855            if (state_quad[quad]==UNDEFINED)
 856               {
 857               Edge* ted [QUAD4];
 858               int etat = AS_IS;
 859               for (int ned=0 ; ned < QUAD4 ; ned++)
 860                   {
 861                   Edge* edge = quad->getEdge (ned);
 862                   if (state_edge [edge]==AS_IS)
 863                       ted[ned] = edge;
 864                   else
 865                       {
 866                       ted[ned] = new_edge[edge];
 867                       etat = REPLACED;
 868                       }
 869                   }
 870               if (etat==REPLACED)
 871                  {
 872                  Quad* face = new Quad (ted[0], ted[1], ted[2], ted[3]);
 873                  new_quad   [quad] = face;
 874                  grid->addQuad  (face);
 875                  }
 876               state_quad [quad] = etat;
 877               }
 878            }
 879        }
 880
 881    for (int nro=0 ; nro<nbhexas ; nro++)
 882        {
 883        Hexa* hexa = thexas [nro];
 884        for (it_quad=new_quad.begin(); it_quad != new_quad.end() ; ++it_quad)
 885            {
 886            Quad* pile = it_quad->first;
 887            Quad* face = it_quad->second;
 888            hexa->replaceQuad (pile, face);
 889            }
 890
 891        for (it_edge=new_edge.begin(); it_edge != new_edge.end() ; ++it_edge)
 892            {
 893            Edge* zig = it_edge->first;
 894            Edge* zag = it_edge->second;
 895            hexa->replaceEdge (zig, zag);
 896            }
 897
 898        for (it_vertex=new_vertex.begin();
 899             it_vertex != new_vertex.end() ; ++it_vertex)
 900            {
 901            Vertex* flip = it_vertex->first;
 902            Vertex* flop = it_vertex->second;
 903            hexa->replaceVertex (flip, flop);
 904            }
 905        }
 906
 907    Real3  center0, center1;
 908    Matrix matrix;
 909    Hexa*  hexa0 = NULL;
 910    Hexa*  hexa1 = NULL;
 911    for (int nro=0 ; nro<=nbhexas ; nro++)
 912        {
 913        hexa0 = hexa1;
 914        if (nro==0)
 915           {
 916           hexa1 = thexas [nro];
 917           hexa1->getCenter (center1);
 918           }
 919        else if (nro==nbhexas)
 920           {
 921           hexa1->getCenter (center1);
 922           }
 923        else
 924           {
 925           hexa1 = thexas [nro];
 926           hexa0->getCenter (center0);
 927           hexa1->getCenter (center1);
 928           for (int nc=0 ; nc<DIM3 ; nc++)
 929               center1[nc] = (center0[nc] + center1[nc])/2;
 930           }
 931
 932        Vertex* node  = grid->getVertex (nro);
 933        matrix.defScale (center1, 0.55);
 934        matrix.perform  (node);
 935        }
 936
 937    DumpReturn (grid);
 938    return      grid;
 939 }
 940 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
 941 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
 942 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
 943 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
 944 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
 945 // ========================================================= prod_vectoriel
 946 double*  prod_vectoriel (Edge* e1, Edge* e2, double prod[])
 947 {
 948    prod [dir_x] = prod [dir_y] = prod [dir_z] = 0;
 949    if (e1==NULL || e2==NULL)
 950       return prod;
 951
 952    Real3 v1, v2;
 953    e1->getVector (v1);
 954    e2->getVector (v2);
 955
 956    prod [dir_x] = v1[dir_y] * v2[dir_z] - v2[dir_y] * v1[dir_z];
 957    prod [dir_y] = v1[dir_z] * v2[dir_x] - v2[dir_z] * v1[dir_x];
 958    prod [dir_z] = v1[dir_x] * v2[dir_y] - v2[dir_x] * v1[dir_y];
 959
 960    return prod;
 961 }
 962 // ========================================================= permuter_edges
 963 void permuter_edges (Edge* &e1, Edge* &e2, cpchar nm1, cpchar nm2)
 964 {
 965    if (db && nm1!=NULL)
 966       {
 967       printf (" ... permuter_edges %s = %s et %s = %s\n",
 968                     nm1, e1->getName(), nm2, e2->getName() );
 969       }
 970
 971    Edge* foo = e1;
 972    e1  = e2;
 973    e2  = foo;
 974 }
 975 END_NAMESPACE_HEXA