src/HEXABLOCK/HexDocument_quads.cxx

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