2 // C++ : Classe Document : Methodes internes 2011
4 // Copyright (C) 2009-2013 CEA/DEN, EDF R&D
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.
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.
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
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 #include "HexDocument.hxx"
28 #include "HexVertex.hxx"
29 #include "HexEdge.hxx"
30 #include "HexQuad.hxx"
31 #include "HexHexa.hxx"
35 #include "HexAnaQuads.hxx"
36 #include "HexElements.hxx"
37 #include "HexCramer.hxx"
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[]);
45 static bool db = false;
47 // ======================================================== copyDocument
48 Document* Document::copyDocument ()
50 string nom = "CopyOf_";
53 Document* clone = new Document (nom.c_str());
55 for (EltBase* elt = doc_first_elt[EL_VERTEX]->next (); elt!=NULL;
58 if (elt !=NULL && elt->isHere())
60 Vertex* node = static_cast <Vertex*> (elt);
61 node->duplicate (clone);
65 for (int type=EL_EDGE ; type <= EL_HEXA ; type++)
67 for (EltBase* elt = doc_first_elt[type]->next (); elt!=NULL;
70 if (elt !=NULL && elt->isHere())
75 for (int nro=0 ; nro<nbr_laws ; nro++)
77 Law* law = new Law (doc_laws [nro]);
78 clone->doc_laws.push_back (law);
83 // ----------------------------------------------------------------------------
84 // ======================================================== countUsedHexa
85 int Document::countUsedHexa ()
90 return nbr_used_hexas;
92 // ======================================================== countUsedQuad
93 int Document::countUsedQuad ()
98 return nbr_used_quads;
100 // ======================================================== countUsedEdge
101 int Document::countUsedEdge ()
106 return nbr_used_edges;
108 // ======================================================== countUsedVertex
109 int Document::countUsedVertex ()
114 return nbr_used_vertex;
117 // ======================================================== getUsedHexa
118 Hexa* Document::getUsedHexa (int nro)
123 if (nro<0 || nro >= nbr_used_hexas)
126 return doc_used_hexas [nro];
128 // ======================================================== getUsedQuad
129 Quad* Document::getUsedQuad (int nro)
134 if (nro<0 || nro >= nbr_used_quads)
137 return doc_used_quads [nro];
139 // ======================================================== getUsedEdge
140 Edge* Document::getUsedEdge (int nro)
145 if (nro<0 || nro >= nbr_used_edges)
148 return doc_used_edges [nro];
150 // ======================================================== getUsedVertex
151 Vertex* Document::getUsedVertex (int nro)
156 if (nro<0 || nro >= nbr_used_vertex)
159 return doc_used_vertex [nro];
161 // ======================================================== renumeroter
162 void Document::renumeroter ()
164 count_modified = false;
165 // -- 1) Raz numerotation precedente
166 markAll (NO_COUNTED);
168 doc_used_hexas .clear ();
169 doc_used_quads .clear ();
170 doc_used_edges .clear ();
171 doc_used_vertex.clear ();
173 for (EltBase* elt = doc_first_elt[EL_HEXA]->next (); elt!=NULL;
176 if (elt!=NULL && elt->isHere())
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);
183 for (int nb=0 ; nb<HE_MAXI ; nb++)
184 cell->getEdge (nb)->setMark (IS_USED);
186 for (int nb=0 ; nb<HV_MAXI ; nb++)
187 cell->getVertex (nb)->setMark (IS_USED);
191 for (EltBase* elt = doc_first_elt[EL_QUAD]->next (); elt!=NULL;
194 if (elt!=NULL && elt->isHere() && elt->getMark()==IS_USED)
196 Quad* cell = static_cast <Quad*> (elt);
197 doc_used_quads.push_back (cell);
201 for (EltBase* elt = doc_first_elt[EL_EDGE]->next (); elt!=NULL;
204 if (elt!=NULL && elt->isHere() && elt->getMark()==IS_USED)
206 Edge* cell = static_cast <Edge*> (elt);
207 doc_used_edges.push_back (cell);
211 for (EltBase* elt = doc_first_elt[EL_VERTEX]->next (); elt!=NULL;
214 if (elt!=NULL && elt->isHere() && elt->getMark()==IS_USED)
216 Vertex* cell = static_cast <Vertex*> (elt);
217 doc_used_vertex.push_back (cell);
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 ();
226 // ---------------------------------------------------------------
227 // ============================================================== addHexa2quads
228 Hexa* Document::addHexa2Quads (Quad* q1, Quad* q2)
230 AnaQuads ana_quads (q1, q2);
233 if (ana_quads.status != HOK)
236 else if (ana_quads.nbr_aretes == 0)
237 hexa = addHexaQuadsAB (ana_quads);
239 else if (ana_quads.nbr_aretes == 1)
240 hexa = addHexaQuadsAC (ana_quads);
244 // ============================================================= addHexa3quads
245 Hexa* Document::addHexa3Quads (Quad* q1, Quad* q2, Quad* q3)
247 AnaQuads ana_quads (q1, q2, q3);
250 if (ana_quads.status != HOK)
253 else if (ana_quads.nbr_aretes == 2)
254 hexa = addHexaQuadsACD (ana_quads);
256 else if (ana_quads.nbr_aretes == 3)
257 hexa = addHexaQuadsACE (ana_quads);
261 // ============================================================= addHexa4quads
262 Hexa* Document::addHexa4Quads (Quad* q1, Quad* q2, Quad* q3, Quad* q4)
264 AnaQuads ana_quads (q1, q2, q3, q4);
267 if (ana_quads.status != HOK)
270 else if (ana_quads.nbr_aretes == 4)
271 hexa = addHexaQuadsABCD (ana_quads);
273 else if (ana_quads.nbr_aretes == 5)
274 hexa = addHexaQuadsACDE (ana_quads);
278 // ============================================================== addHexa5quads
279 Hexa* Document::addHexa5Quads (Quad* q1, Quad* q2, Quad* q3, Quad* q4, Quad* q5)
281 AnaQuads ana_quads (q1, q2, q3, q4, q5);
282 if (ana_quads.status != HOK)
284 else if (ana_quads.nbr_aretes != 8)
288 for (int nquad=0 ; nquad < ana_quads.nbr_quads ; nquad++)
289 if (ana_quads.inter_nbre [nquad] == 4)
292 if (qbase == NOTHING)
297 for (int nedge=0 ; nedge < QUAD4 ; nedge++)
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);
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]);
312 // ---------------------------------------------------------------
313 // ========================================================== addHexaquadsAB
314 Hexa* Document::addHexaQuadsAB (AnaQuads& strquads)
316 Quad* q_a = strquads.tab_quads[0];
317 Quad* q_b = strquads.tab_quads[1];
322 for (int is = 0 ; is<2 ; is++)
325 for (int ndec = 0 ; ndec<QUAD4 ; ndec++)
328 for (int na = 0 ; na<QUAD4 ; na++)
330 int nb = (ndec + QUAD4 + ns*na) MODULO QUAD4;
331 dist += distance (q_a->getVertex (na),
332 q_b->getVertex (nb));
334 if (ndec==0 && is==0)
351 for (int na = 0 ; na<QUAD4 ; na++)
353 int nb = (decal + QUAD4 + sens*na) MODULO QUAD4;
354 tedge [na] = new Edge (q_a->getVertex (na), q_b->getVertex (nb));
357 for (int nal = 0 ; nal<QUAD4 ; nal++)
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]);
369 Hexa* hexa = new Hexa (q_a, q_b, tquad[0], tquad[2], tquad[1], tquad[3]);
372 // ========================================================== addHexaquadsAC
373 Hexa* Document::addHexaQuadsAC (AnaQuads& strquads)
375 Quad* q_a = strquads.tab_quads[0];
376 Quad* q_c = strquads.tab_quads[1];
378 Vertex* tv_bdx [V_TWO]; // x = e ou f
383 int neda0 = strquads.inter_edge [0] [1];
385 Edge* e_ac = q_a->getEdge (neda0);
387 for (int ns=V_AMONT; ns<=V_AVAL ; ns++)
389 Vertex* vx1 = e_ac->getVertex (ns);
390 Vertex* vx2 = e_ac->getVertex (1-ns);
392 int nda2 = q_a->indexVertex (vx2);
393 nda2 = (nda2 +2) MODULO QUAD4;
395 int ndc2 = q_c->indexVertex (vx2);
396 ndc2 = (ndc2 +2) MODULO QUAD4;
398 Vertex* vxa = q_a->getVertex (nda2);
399 Vertex* vxc = q_c->getVertex (ndc2);
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,
406 Edge* edga = q_a->findEdge (vx1, vxa);
407 Edge* edgc = q_c->findEdge (vx1, vxc);
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);
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);
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]);
422 Hexa* hexa = new Hexa (q_a, q_b, q_c, q_d, tq_ef[V_AMONT], tq_ef[V_AVAL]);
425 // ========================================================= addHexaquadsACE
426 // ==== Construction d'un hexaedre a partir d'un triedre
428 6=bed +----bd-----+ bdf=7
432 4=bce +----bc-----+...|...bcf=5
436 2=ade...|...+----ad-|---+ adf=3 | y
440 0=ace +----ac-----+ acf=1 +-----> x
442 On connait les faces A, C, E
443 Il faut determiner le point bdf, intersection des faces B, D, F
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
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
454 Scalaire ((M-bce), norm_b) = 0
455 Scalaire ((M-ade), morm_d) = 0
456 Scalaire ((M-acf), morm_f) = 0
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)
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
469 * ----------------------------------------------------- */
470 Hexa* Document::addHexaQuadsACE (AnaQuads& strquads)
472 Real3 v_bce, v_ade, v_acf, v_bdf; // Sommets
473 Real3 norm_b, norm_d, norm_f; // Normales aux faces
476 Quad* q_a = strquads.tab_quads[0];
477 Quad* q_c = strquads.tab_quads[1];
478 Quad* q_e = strquads.tab_quads[2];
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);
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);
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);
492 e_ce->commonPoint (e_bc, v_bce);
493 e_ae->commonPoint (e_ad, v_ade);
494 e_ac->commonPoint (e_af, v_acf);
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
500 Real3 membre2 = { prod_scalaire (v_bce, norm_b),
501 prod_scalaire (v_ade, norm_d),
502 prod_scalaire (v_acf, norm_f) };
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);
511 printf (" addHexaQuadsACE : systeme impossible\n");
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);
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);
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);
528 Hexa* hexa = new Hexa (q_a, q_b, q_c, q_d, q_e, q_f);
531 // ========================================================= addHexaquadsACD
532 // ==== Construction d'un hexaedre a partir d'un U
533 Hexa* Document::addHexaQuadsACD (AnaQuads& strquads)
536 for (int np=0 ; np<3 && pos_a==NOTHING ; np++)
537 if (strquads.inter_nbre[np]==2)
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];
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
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
556 Edge* e_bc = q_c->getEdge ((nc_ac + 2) MODULO QUAD4);
557 Edge* e_bd = q_d->getEdge ((nd_ad + 2) MODULO QUAD4);
559 Edge* e_ce = q_c->getEdge ((nc_ac + 1) MODULO QUAD4);
560 Edge* e_cf = q_c->getEdge ((nc_ac + 3) MODULO QUAD4);
562 Edge* e_de = q_d->getEdge ((nd_ad + 1) MODULO QUAD4);
563 Edge* e_df = q_d->getEdge ((nd_ad + 3) MODULO QUAD4);
565 Vertex* v_ace = e_ae->commonVertex (e_ce);
566 Vertex* v_ade = e_ae->commonVertex (e_de);
569 permuter_edges (e_ce, e_cf);
570 v_ace = e_ae->commonVertex (e_ce);
575 permuter_edges (e_de, e_df);
576 v_ade = e_ae->commonVertex (e_de);
579 Vertex* v_acf = e_af->commonVertex (e_cf);
580 Vertex* v_adf = e_af->commonVertex (e_df);
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);
587 Edge* e_be = new Edge (v_bce, v_bde);
588 Edge* e_bf = new Edge (v_bcf, v_bdf);
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);
594 Hexa* hexa = new Hexa (q_a, q_b, q_c, q_d, q_e, q_f);
597 // ========================================================= addHexaquadsABCD
598 Hexa* Document::addHexaQuadsABCD (AnaQuads& strquads)
605 for (int np=1 ; np<4 ; np++)
607 if (strquads.inter_edge [pos_a] [np] == NOTHING )
609 else if (pos_c==NOTHING)
615 if (pos_b==NOTHING || pos_c==NOTHING || pos_d==NOTHING)
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];
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
628 Edge* e_ae = q_a->getEdge ((na_ac + 1) MODULO QUAD4);
629 Edge* e_af = q_a->getEdge ((na_ac + 3) MODULO QUAD4);
631 Edge* e_ce = q_c->getEdge ((nc_ac + 1) MODULO QUAD4);
632 Edge* e_cf = q_c->getEdge ((nc_ac + 3) MODULO QUAD4);
634 Edge* e_de = q_d->getEdge ((nd_ad + 1) MODULO QUAD4);
635 Edge* e_df = q_d->getEdge ((nd_ad + 3) MODULO QUAD4);
637 Edge* e_be = q_b->getEdge ((nb_bc + 1) MODULO QUAD4);
638 Edge* e_bf = q_b->getEdge ((nb_bc + 3) MODULO QUAD4);
657 if (e_ae->commonVertex (e_ce) == NULL)
658 PermuterEdges (e_ce, e_cf);
660 if (e_ae->commonVertex (e_de) == NULL)
661 PermuterEdges (e_de, e_df);
663 if (e_ce->commonVertex (e_be) == NULL)
664 PermuterEdges (e_be, e_bf);
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);
669 Hexa* hexa = new Hexa (q_a, q_b, q_c, q_d, q_e, q_f);
672 // ========================================================= addHexaquadsACDE
673 Hexa* Document::addHexaQuadsACDE (AnaQuads& strquads)
680 for (int np=0 ; np<4 ; np++)
681 if (strquads.inter_nbre[np]==3)
683 if (pos_a == NOTHING) pos_a = np;
686 else if (strquads.inter_nbre[np]==2)
688 if (pos_c == NOTHING) pos_c = np;
692 if (pos_a==NOTHING || pos_c==NOTHING || pos_d==NOTHING || pos_e==NOTHING)
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];
700 int na_ac = strquads.inter_edge[pos_a][pos_c]; // Nro dans q_a de e_ac
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
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
708 int ne_ae = strquads.inter_edge[pos_e][pos_a]; // Nro dans q_e de e_ae
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);
717 Vertex* v_bcf = e_cf->opposedVertex (e_cf->commonVertex (e_af));
718 Vertex* v_bdf = e_df->opposedVertex (e_df->commonVertex (e_af));
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);
724 Hexa* hexa = new Hexa (q_a, q_b, q_c, q_d, q_e, q_f);
727 // ======================================================== revolutionQuads
728 Elements* Document::revolutionQuads (Quads& start, Vertex* center, Vector* axis,
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;
736 Elements* prisme = new Elements (this);
737 prisme->revolutionQuads (start, center, axis, angles);
740 // ======================================================== makeCylindricals
741 Elements* Document::makeCylindricals (Vertex* c, Vector* b, Vector* h,
742 RealVector& tdr, RealVector& tda, RealVector& tdl, bool fill)
744 Elements* grille = new Elements (this);
745 grille->makeCylindricalGrid (c, b, h, tdr, tda, tdl, fill);
748 // ========================================================= replace
749 Elements* Document::replace (Quads& pattern, Vertex* p1, Vertex* c1,
750 Vertex* p2, Vertex* c2, Vertex* p3, Vertex* c3)
752 Elements* t_hexas = new Elements (this);
753 int ier = t_hexas->replaceHexas (pattern, p1, c1, p2, c2, p3, c3);
756 cout << " **** Error in Document::replace\n" << endl;
757 t_hexas->setError (ier);
761 // ========================================================= print_replace
762 void print_replace (Edge* zig, Edge* zag)
764 cout << zig->getName() << " = (" << zig->getVertex(0)->getName()
765 << ", " << zig->getVertex(1)->getName() << ") est clone en ";
767 cout << zag->getName() << " = (" << zag->getVertex(0)->getName()
768 << ", " << zag->getVertex(1)->getName() << ")" << endl;
770 // ========================================================= only_in_hexas
771 bool only_in_hexas (Hexas& thexas, Quad* quad)
773 int nbhexas = thexas.size();
774 int nbp = quad->getNbrParents();
775 for (int nh=0 ; nh <nbp ; nh++)
778 Hexa* hexa = quad->getParent (nh);
779 for (int nc=0 ; pasla && nc < nbhexas ; nc++)
780 pasla = hexa != thexas [nc];
786 // ========================================================= only_in_hexas
787 bool only_in_hexas (Hexas& thexas, Edge* edge)
790 int nbp = edge->getNbrParents();
791 for (int nq=0 ; nq <nbp ; nq++)
793 Quad* quad = edge->getParent (nq);
794 if (NOT only_in_hexas (thexas, quad))
796 cout << " ... inMoreHexas " << edge->makeDefinition() << endl;
800 cout << " ... only_in_hexas " << edge->makeDefinition() << endl;
803 // ========================================================= replace_vertex
804 void replace_vertex (Hexas& thexas, Vertex* node, Vertex* par)
806 int nbh = thexas.size();
807 for (int nh=0 ; nh <nbh ; nh++)
808 thexas[nh]->replaceVertex (node, par);
810 // ========================================================= disconnectEdges
811 Elements* Document::disconnectEdges (Hexas& thexas, Edges& tedges)
813 int nbedges = tedges.size();
814 int nbhexas = thexas.size();
817 cout << " +++ Disconnect Edges" << endl;
819 if (nbhexas != nbedges)
821 cout << " **** Error in Document::disconnectEdges\n" << endl;
822 cout << " **** Number of Edges and number of Hexas are different\n"
828 Elements* grille = disconnectEdge (thexas[0], tedges[0]);
832 for (int nro=0 ; nro<nbedges ; nro++)
834 if (BadElement (tedges[nro]))
836 cout << " **** Eddge number " << nro+1 << " is incorrect"
840 if (BadElement (thexas[nro]))
842 cout << " **** Hexa number " << nro+1 << " is incorrect"
847 cout << nro+1 << " hexa = " << thexas[nro]->getName ()
848 << ", edge = " << tedges[nro]->getName ()
849 << " = (" << tedges[nro]->getVertex(0)->getName ()
850 << ", " << tedges[nro]->getVertex(1)->getName ()
854 for (int nro=0 ; nro<nbhexas ; nro++)
856 int ned = thexas[nro]->findEdge (tedges[nro]);
859 cout << " **** Edge number " << nro+1
860 << " doesnt belong to correspondant hexa" << endl;
865 vector <Vertex*> tvertex (nbedges+1);
867 for (int nro=1 ; nro<nbedges ; nro++)
869 tvertex[nro] = tedges[nro]->commonVertex (tedges[nro-1]);
870 if (tvertex[nro]==NULL)
872 cout << " **** Edge number " << nro
873 << " doesnt intesect next edge" << endl;
878 int nv0 = tedges[0] ->inter (tedges[1]);
879 int nvn = tedges[nbedges-1]->inter (tedges[nbedges-2]);
880 tvertex [0] = tedges[0] ->getVertex (1-nv0);
881 tvertex [nbedges] = tedges[nbedges-1]->getVertex (1-nvn);
883 for (int nro=0 ; nro<nbhexas ; nro++)
885 int ned = thexas[nro]->findEdge (tedges[nro]);
888 cout << " **** Edge number " << nro+1
889 << " doesnt belong to correspondant hexa" << endl;
893 // Fin des controles, on peut y aller ...
895 map <Edge*, int> state_edge;
896 map <Quad*, int> state_quad;
897 enum { UNDEFINED, REPLACED, AS_IS };
899 map <Vertex*, Vertex*> new_vertex;
900 map <Edge*, Edge*> new_edge;
901 map <Quad*, Quad*> new_quad;
903 map <Vertex*, Vertex*> :: iterator it_vertex;
904 map <Edge*, Edge*> :: iterator it_edge;
905 map <Quad*, Quad*> :: iterator it_quad;
907 #define VertexIsNew(v) (it_vertex=new_vertex.find(v))!=new_vertex.end()
909 Elements* nouveaux = new Elements (this);
910 Vertex* node1 = NULL;
912 for (int nro=0 ; nro<=nbedges ; nro++)
914 Vertex* node0 = node1;
915 node1 = new Vertex (tvertex[nro]);
916 nouveaux->addVertex (node1);
917 new_vertex [tvertex[nro]] = node1;
920 cout << nro << " : " << tvertex[nro]->getName()
921 << " est clone en " << node1->getName() << endl;
926 Edge* edge = new Edge (node0, node1);
927 nouveaux->addEdge (edge);
928 new_edge [tedges[nro-1]] = edge;
929 state_edge [tedges[nro-1]] = REPLACED;
931 print_replace (tedges[nro-1], edge);
936 cout << "_____________________________ Autres substitutions" << endl;
938 // Un edge non remplace, qui contient un vertex remplace
939 // commun a plus de 2 faces (donc appartenant a un autre hexa)
940 // doit etre duplique
942 for (int nro=0 ; nro<nbhexas ; nro++)
944 Hexa* hexa = thexas [nro];
945 for (int nro=0 ; nro<HE_MAXI ; nro++)
947 Edge* edge = hexa->getEdge (nro);
948 if (state_edge[edge]==UNDEFINED)
950 Vertex* v1 = edge->getVertex (V_AMONT);
951 Vertex* v2 = edge->getVertex (V_AVAL);
953 if (VertexIsNew (v1))
955 if (only_in_hexas (thexas, edge))
957 replace_vertex (thexas, v1, new_vertex[v1]);
961 v1 = new_vertex [v1];
963 else if (VertexIsNew (v2))
965 if (only_in_hexas (thexas, edge))
967 replace_vertex (thexas, v2, new_vertex[v2]);
971 v2 = new_vertex [v2];
978 Edge* arete = new Edge (v1, v2);
979 new_edge [edge] = arete;
980 nouveaux->addEdge (arete);
982 print_replace (edge, arete);
984 state_edge [edge] = etat;
989 // Un quad non remplace, qui contient un edge remplace
990 // commun a plus de 2 Hexas
991 // doit etre duplique
993 for (int nro=0 ; nro<nbhexas ; nro++)
995 Hexa* hexa = thexas [nro];
996 for (int nro=0 ; nro<HQ_MAXI ; nro++)
998 Quad* quad = hexa->getQuad (nro);
999 if (state_quad[quad]==UNDEFINED)
1003 for (int ned=0 ; ned < QUAD4 ; ned++)
1005 Edge* edge = quad->getEdge (ned);
1006 if (state_edge [edge]==AS_IS)
1010 ted[ned] = new_edge[edge];
1016 Quad* face = new Quad (ted[0], ted[1], ted[2], ted[3]);
1017 new_quad [quad] = face;
1018 nouveaux->addQuad (face);
1020 state_quad [quad] = etat;
1025 for (int nro=0 ; nro<nbhexas ; nro++)
1027 Hexa* hexa = thexas [nro];
1028 for (it_quad=new_quad.begin(); it_quad != new_quad.end() ; ++it_quad)
1030 Quad* pile = it_quad->first;
1031 Quad* face = it_quad->second;
1032 hexa->replaceQuad (pile, face);
1035 for (it_edge=new_edge.begin(); it_edge != new_edge.end() ; ++it_edge)
1037 Edge* zig = it_edge->first;
1038 Edge* zag = it_edge->second;
1039 hexa->replaceEdge (zig, zag);
1042 for (it_vertex=new_vertex.begin();
1043 it_vertex != new_vertex.end() ; ++it_vertex)
1045 Vertex* flip = it_vertex->first;
1046 Vertex* flop = it_vertex->second;
1047 hexa->replaceVertex (flip, flop);
1051 Real3 center0, center1;
1055 for (int nro=0 ; nro<=nbhexas ; nro++)
1060 hexa1 = thexas [nro];
1061 hexa1->getCenter (center1);
1063 else if (nro==nbhexas)
1065 hexa1->getCenter (center1);
1069 hexa1 = thexas [nro];
1070 hexa0->getCenter (center0);
1071 hexa1->getCenter (center1);
1072 for (int nc=0 ; nc<DIM3 ; nc++)
1073 center1[nc] = (center0[nc] + center1[nc])/2;
1076 Vertex* node = nouveaux->getVertex (nro);
1077 matrix.defScale (center1, 0.55);
1078 matrix.perform (node);
1083 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1084 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1085 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1086 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1087 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1088 // ========================================================= prod_vectoriel
1089 double* prod_vectoriel (Edge* e1, Edge* e2, double prod[])
1091 prod [dir_x] = prod [dir_y] = prod [dir_z] = 0;
1092 if (e1==NULL || e2==NULL)
1099 prod [dir_x] = v1[dir_y] * v2[dir_z] - v2[dir_y] * v1[dir_z];
1100 prod [dir_y] = v1[dir_z] * v2[dir_x] - v2[dir_z] * v1[dir_x];
1101 prod [dir_z] = v1[dir_x] * v2[dir_y] - v2[dir_x] * v1[dir_y];
1105 // ========================================================= permuter_edges
1106 void permuter_edges (Edge* &e1, Edge* &e2, cpchar nm1, cpchar nm2)
1108 if (db && nm1!=NULL)
1110 printf (" ... permuter_edges %s = %s et %s = %s\n",
1111 nm1, e1->getName(), nm2, e2->getName() );