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 // ========================================================= replaceHexa
 595 Elements* Document::replaceHexa (Quads pattern, Vertex* p1, Vertex* c1,
 596                              Vertex* p2, Vertex* c2, Vertex* p3, Vertex* c3)
 597 {
 598    DumpStart ("replaceHexa", 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 // ========================================================= repla
 612 Elements* Document::replace (Quads motif, Quads cible, Vertex* p1, Vertex* c1,                               Vertex* p2, Vertex* c2)
 613 {
 614    DumpStart ("replace", motif << cible << p1 << c1 << p2 << c2);
 615
 616    Elements* t_hexas = new Elements (this);
 617
 618    //if (BadElement (edge))
 619       {
 620       t_hexas->setError (HERR);
 621       Mess << "This function is not yet implemented";
 622       }
 623
 624    DumpReturn (t_hexas);
 625    return t_hexas;
 626 }
 627 // ========================================================= print_replace
 628 void print_replace (Edge* zig, Edge*  zag)
 629 {
 630    cout << zig->getName() << " = (" << zig->getVertex(0)->getName()
 631         << ", " << zig->getVertex(1)->getName() << ") est clone en ";
 632
 633    cout << zag->getName() << " = (" << zag->getVertex(0)->getName()
 634         << ", " << zag->getVertex(1)->getName() << ")" << endl;
 635 }
 636 // ========================================================= only_in_hexas
 637 bool only_in_hexas (Hexas& thexas, Quad* quad)
 638 {
 639    int nbhexas = thexas.size();
 640    int nbp     = quad->getNbrParents();
 641    for (int nh=0 ; nh <nbp ; nh++)
 642        {
 643        bool pasla = true;
 644        Hexa* hexa = quad->getParent (nh);
 645        for (int nc=0 ; pasla && nc < nbhexas ; nc++)
 646            pasla = hexa != thexas [nc];
 647        if (pasla)
 648            return false;
 649        }
 650    return true;
 651 }
 652 // ========================================================= only_in_hexas
 653 bool only_in_hexas (Hexas& thexas, Edge*  edge)
 654 {
 655    return false;
 656    int nbp = edge->getNbrParents();
 657    for (int nq=0 ; nq <nbp ; nq++)
 658        {
 659        Quad* quad = edge->getParent   (nq);
 660        if (NOT only_in_hexas (thexas, quad))
 661           {
 662           cout << " ... inMoreHexas " << edge->makeDefinition() << endl;
 663           return false;
 664           }
 665        }
 666    cout << " ... only_in_hexas " << edge->makeDefinition() << endl;
 667    return true;
 668 }
 669 // ========================================================= replace_vertex
 670 void replace_vertex (Hexas& thexas, Vertex* node,  Vertex* par)
 671 {
 672    int nbh = thexas.size();
 673    for (int nh=0 ; nh <nbh ; nh++)
 674        thexas[nh]->replaceVertex (node, par);
 675 }
 676 // ========================================================= disconnectEdges
 677 Elements* Document::disconnectEdges (Hexas thexas, Edges  tedges)
 678 {
 679    DumpStart ("disconnectEdges",  thexas << tedges);
 680
 681    if (db)
 682       cout << " +++ Disconnect Edges" << endl;
 683
 684    Elements* grid  = new Elements (this);
 685
 686    grid->checkDisco (thexas, tedges);
 687
 688    int nbedges = tedges.size();
 689    int nbhexas = thexas.size();
 690
 691    if (nbhexas != nbedges)
 692       {
 693       cout << " **** Error in Document::disconnectEdges\n" << endl;
 694       cout << " **** Number of Edges and number of Hexas are different\n"
 695            << endl;
 696       return NULL;
 697       }
 698    else if (nbhexas==1)
 699       {
 700       update ();
 701       thexas[0]->disconnectEdge (tedges[0], grid);
 702       return    grid;
 703       }
 704
 705    for (int nro=0 ; nro<nbedges ; nro++)
 706        {
 707        if (BadElement (tedges[nro]))
 708           {
 709           cout << " **** Eddge number " << nro+1 << " is incorrect"
 710                << endl;
 711           return NULL;
 712           }
 713        if (BadElement (thexas[nro]))
 714           {
 715           cout << " **** Hexa number " << nro+1 << " is incorrect"
 716                << endl;
 717           return NULL;
 718           }
 719        if (db)
 720           cout << nro+1 << " hexa = " << thexas[nro]->getName ()
 721                         << ", edge = " << tedges[nro]->getName ()
 722                         << " = (" << tedges[nro]->getVertex(0)->getName ()
 723                         << ", "   << tedges[nro]->getVertex(1)->getName ()
 724                         << ")" << endl;
 725        }
 726
 727    for (int nro=0 ; nro<nbhexas ; nro++)
 728        {
 729        int ned = thexas[nro]->findEdge (tedges[nro]);
 730        if (ned==NOTHING)
 731           {
 732           cout << " **** Edge number " << nro+1
 733                << " doesnt belong to correspondant hexa" << endl;
 734           return NULL;
 735           }
 736        }
 737
 738    vector <Vertex*> tvertex (nbedges+1);
 739
 740    for (int nro=1 ; nro<nbedges ; nro++)
 741        {
 742        tvertex[nro] = tedges[nro]->commonVertex (tedges[nro-1]);
 743        if (tvertex[nro]==NULL)
 744           {
 745           cout << " **** Edge number " << nro
 746                << " doesnt intesect next edge" << endl;
 747           return NULL;
 748           }
 749        }
 750
 751    int nv0 = tedges[0]        ->inter (tedges[1]);
 752    int nvn = tedges[nbedges-1]->inter (tedges[nbedges-2]);
 753    tvertex [0]       = tedges[0]        ->getVertex (1-nv0);
 754    tvertex [nbedges] = tedges[nbedges-1]->getVertex (1-nvn);
 755
 756    for (int nro=0 ; nro<nbhexas ; nro++)
 757        {
 758        int ned = thexas[nro]->findEdge (tedges[nro]);
 759        if (ned==NOTHING)
 760           {
 761           cout << " **** Edge number " << nro+1
 762                << " doesnt belong to correspondant hexa" << endl;
 763           return NULL;
 764           }
 765        }
 766                       // Fin des controles, on peut y aller ...
 767
 768    map <Edge*, int> state_edge;
 769    map <Quad*, int> state_quad;
 770    enum { UNDEFINED, REPLACED, AS_IS };
 771
 772    map <Vertex*, Vertex*> new_vertex;
 773    map <Edge*,   Edge*>   new_edge;
 774    map <Quad*,   Quad*>   new_quad;
 775
 776    map <Vertex*, Vertex*> :: iterator it_vertex;
 777    map <Edge*,   Edge*>   :: iterator it_edge;
 778    map <Quad*,   Quad*>   :: iterator it_quad;
 779
 780 #define VertexIsNew(v) (it_vertex=new_vertex.find(v))!=new_vertex.end()
 781
 782    Vertex*   node1     = NULL;
 783
 784    for (int nro=0 ; nro<=nbedges ; nro++)
 785        {
 786        Vertex* node0 = node1;
 787        node1 = new Vertex (tvertex[nro]);
 788        grid->addVertex  (node1);
 789        new_vertex [tvertex[nro]] = node1;
 790        if (db)
 791           {
 792           cout << nro << " : "         << tvertex[nro]->getName()
 793                << " est clone en " << node1->getName() << endl;
 794           }
 795
 796        if (nro>0)
 797           {
 798           Edge* edge = new Edge (node0, node1);
 799           grid->addEdge  (edge);
 800           new_edge   [tedges[nro-1]] = edge;
 801           state_edge [tedges[nro-1]] = REPLACED;
 802           if (db)
 803              print_replace (tedges[nro-1], edge);
 804           }
 805        }
 806
 807    if (db)
 808       cout << "_____________________________ Autres substitutions" << endl;
 809
 810    // Un edge non remplace, qui contient un vertex remplace
 811    //         commun a plus de 2 faces (donc appartenant a un autre hexa)
 812    //         doit etre duplique
 813
 814    for (int nro=0 ; nro<nbhexas ; nro++)
 815        {
 816        Hexa* hexa = thexas [nro];
 817        for (int nro=0 ; nro<HE_MAXI ; nro++)
 818            {
 819            Edge* edge = hexa->getEdge (nro);
 820            if (state_edge[edge]==UNDEFINED)
 821               {
 822               Vertex* v1 = edge->getVertex (V_AMONT);
 823               Vertex* v2 = edge->getVertex (V_AVAL);
 824               int etat = REPLACED;
 825               if (VertexIsNew (v1))
 826                  {
 827                  if (only_in_hexas (thexas, edge))
 828                     {
 829                     replace_vertex (thexas, v1, new_vertex[v1]);
 830                     etat = AS_IS;
 831                     }
 832                  else
 833                     v1 = new_vertex [v1];
 834                  }
 835               else if (VertexIsNew (v2))
 836                  {
 837                  if (only_in_hexas (thexas, edge))
 838                     {
 839                     replace_vertex (thexas, v2, new_vertex[v2]);
 840                     etat = AS_IS;
 841                     }
 842                  else
 843                     v2 = new_vertex [v2];
 844                  }
 845               else
 846                  etat = AS_IS;
 847
 848               if (etat==REPLACED)
 849                  {
 850                  Edge* arete = new Edge (v1, v2);
 851                  new_edge   [edge] = arete;
 852                  grid->addEdge  (arete);
 853                  if (db)
 854                     print_replace (edge, arete);
 855                  }
 856               state_edge [edge] = etat;
 857               }
 858            }
 859        }
 860
 861    // Un quad non remplace, qui contient un edge remplace
 862    //         commun a plus de 2 Hexas
 863    //         doit etre duplique
 864
 865    for (int nro=0 ; nro<nbhexas ; nro++)
 866        {
 867        Hexa* hexa = thexas [nro];
 868        for (int nro=0 ; nro<HQ_MAXI ; nro++)
 869            {
 870            Quad* quad = hexa->getQuad (nro);
 871            if (state_quad[quad]==UNDEFINED)
 872               {
 873               Edge* ted [QUAD4];
 874               int etat = AS_IS;
 875               for (int ned=0 ; ned < QUAD4 ; ned++)
 876                   {
 877                   Edge* edge = quad->getEdge (ned);
 878                   if (state_edge [edge]==AS_IS)
 879                       ted[ned] = edge;
 880                   else
 881                       {
 882                       ted[ned] = new_edge[edge];
 883                       etat = REPLACED;
 884                       }
 885                   }
 886               if (etat==REPLACED)
 887                  {
 888                  Quad* face = new Quad (ted[0], ted[1], ted[2], ted[3]);
 889                  new_quad   [quad] = face;
 890                  grid->addQuad  (face);
 891                  }
 892               state_quad [quad] = etat;
 893               }
 894            }
 895        }
 896
 897    for (int nro=0 ; nro<nbhexas ; nro++)
 898        {
 899        Hexa* hexa = thexas [nro];
 900        for (it_quad=new_quad.begin(); it_quad != new_quad.end() ; ++it_quad)
 901            {
 902            Quad* pile = it_quad->first;
 903            Quad* face = it_quad->second;
 904            hexa->replaceQuad (pile, face);
 905            }
 906
 907        for (it_edge=new_edge.begin(); it_edge != new_edge.end() ; ++it_edge)
 908            {
 909            Edge* zig = it_edge->first;
 910            Edge* zag = it_edge->second;
 911            hexa->replaceEdge (zig, zag);
 912            }
 913
 914        for (it_vertex=new_vertex.begin();
 915             it_vertex != new_vertex.end() ; ++it_vertex)
 916            {
 917            Vertex* flip = it_vertex->first;
 918            Vertex* flop = it_vertex->second;
 919            hexa->replaceVertex (flip, flop);
 920            }
 921        }
 922
 923    Real3  center0, center1;
 924    Matrix matrix;
 925    Hexa*  hexa0 = NULL;
 926    Hexa*  hexa1 = NULL;
 927    for (int nro=0 ; nro<=nbhexas ; nro++)
 928        {
 929        hexa0 = hexa1;
 930        if (nro==0)
 931           {
 932           hexa1 = thexas [nro];
 933           hexa1->getCenter (center1);
 934           }
 935        else if (nro==nbhexas)
 936           {
 937           hexa1->getCenter (center1);
 938           }
 939        else
 940           {
 941           hexa1 = thexas [nro];
 942           hexa0->getCenter (center0);
 943           hexa1->getCenter (center1);
 944           for (int nc=0 ; nc<DIM3 ; nc++)
 945               center1[nc] = (center0[nc] + center1[nc])/2;
 946           }
 947
 948        Vertex* node  = grid->getVertex (nro);
 949        matrix.defScale (center1, 0.55);
 950        matrix.perform  (node);
 951        }
 952
 953    DumpReturn (grid);
 954    return      grid;
 955 }
 956 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
 957 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
 958 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
 959 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
 960 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
 961 // ========================================================= prod_vectoriel
 962 double*  prod_vectoriel (Edge* e1, Edge* e2, double prod[])
 963 {
 964    prod [dir_x] = prod [dir_y] = prod [dir_z] = 0;
 965    if (e1==NULL || e2==NULL)
 966       return prod;
 967
 968    Real3 v1, v2;
 969    e1->getVector (v1);
 970    e2->getVector (v2);
 971
 972    prod [dir_x] = v1[dir_y] * v2[dir_z] - v2[dir_y] * v1[dir_z];
 973    prod [dir_y] = v1[dir_z] * v2[dir_x] - v2[dir_z] * v1[dir_x];
 974    prod [dir_z] = v1[dir_x] * v2[dir_y] - v2[dir_x] * v1[dir_y];
 975
 976    return prod;
 977 }
 978 // ========================================================= permuter_edges
 979 void permuter_edges (Edge* &e1, Edge* &e2, cpchar nm1, cpchar nm2)
 980 {
 981    if (db && nm1!=NULL)
 982       {
 983       printf (" ... permuter_edges %s = %s et %s = %s\n",
 984                     nm1, e1->getName(), nm2, e2->getName() );
 985       }
 986
 987    Edge* foo = e1;
 988    e1  = e2;
 989    e2  = foo;
 990 }
 991 END_NAMESPACE_HEXA