src/HEXABLOCK/HexDocument_quads.cxx

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