2 // C++ : Classe Document : Methodes internes 2011
4 // Copyright (C) 2009-2014 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, or (at your option) any later version.
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"
38 #include "HexGlobale.hxx"
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[]);
46 static bool db = false;
48 // ======================================================== copyDocument
49 Document* Document::copyDocument ()
51 string nom = "CopyOf_";
54 Document* clone = new Document (nom.c_str());
56 for (EltBase* elt = doc_first_elt[EL_VERTEX]->next (); elt!=NULL;
59 if (elt !=NULL && elt->isHere())
61 Vertex* node = static_cast <Vertex*> (elt);
62 node->duplicate (clone);
66 for (int type=EL_EDGE ; type <= EL_HEXA ; type++)
68 for (EltBase* elt = doc_first_elt[type]->next (); elt!=NULL;
71 if (elt !=NULL && elt->isHere())
76 for (int nro=0 ; nro<nbr_laws ; nro++)
78 Law* law = new Law (doc_laws [nro]);
79 clone->doc_laws.push_back (law);
84 // ---------------------------------------------------------------
85 // ============================================================== addHexa2quads
86 Hexa* Document::addHexa2Quads (Quad* q1, Quad* q2)
88 DumpStart ("addHexa2Quads", q1 << q2);
89 AnaQuads ana_quads (q1, q2);
92 if (ana_quads.status != HOK)
95 else if (ana_quads.nbr_aretes == 0)
96 hexa = addHexaQuadsAB (ana_quads);
98 else if (ana_quads.nbr_aretes == 1)
99 hexa = addHexaQuadsAC (ana_quads);
104 // ============================================================= addHexa3quads
105 Hexa* Document::addHexa3Quads (Quad* q1, Quad* q2, Quad* q3)
107 DumpStart ("addHexa3Quads", q1 << q2 << q3);
108 AnaQuads ana_quads (q1, q2, q3);
111 if (ana_quads.status != HOK)
114 else if (ana_quads.nbr_aretes == 2)
115 hexa = addHexaQuadsACD (ana_quads);
117 else if (ana_quads.nbr_aretes == 3)
118 hexa = addHexaQuadsACE (ana_quads);
123 // ============================================================= addHexa4quads
124 Hexa* Document::addHexa4Quads (Quad* q1, Quad* q2, Quad* q3, Quad* q4)
126 DumpStart ("addHexa4Quads", q1 << q2 << q3 << q4);
127 AnaQuads ana_quads (q1, q2, q3, q4);
130 if (ana_quads.status != HOK)
133 else if (ana_quads.nbr_aretes == 4)
134 hexa = addHexaQuadsABCD (ana_quads);
136 else if (ana_quads.nbr_aretes == 5)
137 hexa = addHexaQuadsACDE (ana_quads);
142 // ============================================================== addHexa5quads
143 Hexa* Document::addHexa5Quads (Quad* q1, Quad* q2, Quad* q3, Quad* q4, Quad* q5)
145 DumpStart ("addHexa5Quads", q1 << q2 << q3 << q4 << q5);
146 AnaQuads ana_quads (q1, q2, q3, q4, q5);
147 if (ana_quads.status != HOK)
149 else if (ana_quads.nbr_aretes != 8)
153 for (int nquad=0 ; nquad < ana_quads.nbr_quads ; nquad++)
154 if (ana_quads.inter_nbre [nquad] == 4)
157 if (qbase == NOTHING)
162 for (int nedge=0 ; nedge < QUAD4 ; nedge++)
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);
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]);
179 // ---------------------------------------------------------------
180 // ========================================================== addHexaquadsAB
181 Hexa* Document::addHexaQuadsAB (AnaQuads& strquads)
183 Quad* q_a = strquads.tab_quads[0];
184 Quad* q_b = strquads.tab_quads[1];
189 for (int is = 0 ; is<2 ; is++)
192 for (int ndec = 0 ; ndec<QUAD4 ; ndec++)
195 for (int na = 0 ; na<QUAD4 ; na++)
197 int nb = (ndec + QUAD4 + ns*na) MODULO QUAD4;
198 dist += distance (q_a->getVertex (na),
199 q_b->getVertex (nb));
201 if (ndec==0 && is==0)
218 for (int na = 0 ; na<QUAD4 ; na++)
220 int nb = (decal + QUAD4 + sens*na) MODULO QUAD4;
221 tedge [na] = new Edge (q_a->getVertex (na), q_b->getVertex (nb));
224 for (int nal = 0 ; nal<QUAD4 ; nal++)
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]);
236 Hexa* hexa = new Hexa (q_a, q_b, tquad[0], tquad[2], tquad[1], tquad[3]);
239 // ========================================================== addHexaquadsAC
240 Hexa* Document::addHexaQuadsAC (AnaQuads& strquads)
242 Quad* q_a = strquads.tab_quads[0];
243 Quad* q_c = strquads.tab_quads[1];
245 Vertex* tv_bdx [V_TWO]; // x = e ou f
250 int neda0 = strquads.inter_edge [0] [1];
252 Edge* e_ac = q_a->getEdge (neda0);
254 for (int ns=V_AMONT; ns<=V_AVAL ; ns++)
256 Vertex* vx1 = e_ac->getVertex (ns);
257 Vertex* vx2 = e_ac->getVertex (1-ns);
259 int nda2 = q_a->indexVertex (vx2);
260 nda2 = (nda2 +2) MODULO QUAD4;
262 int ndc2 = q_c->indexVertex (vx2);
263 ndc2 = (ndc2 +2) MODULO QUAD4;
265 Vertex* vxa = q_a->getVertex (nda2);
266 Vertex* vxc = q_c->getVertex (ndc2);
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,
273 Edge* edga = q_a->findEdge (vx1, vxa);
274 Edge* edgc = q_c->findEdge (vx1, vxc);
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);
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);
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]);
289 Hexa* hexa = new Hexa (q_a, q_b, q_c, q_d, tq_ef[V_AMONT], tq_ef[V_AVAL]);
292 // ========================================================= addHexaquadsACE
293 // ==== Construction d'un hexaedre a partir d'un triedre
295 6=bed +----bd-----+ bdf=7
299 4=bce +----bc-----+...|...bcf=5
303 2=ade...|...+----ad-|---+ adf=3 | y
307 0=ace +----ac-----+ acf=1 +-----> x
309 On connait les faces A, C, E
310 Il faut determiner le point bdf, intersection des faces B, D, F
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
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
321 Scalaire ((M-bce), norm_b) = 0
322 Scalaire ((M-ade), morm_d) = 0
323 Scalaire ((M-acf), morm_f) = 0
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)
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
336 * ----------------------------------------------------- */
337 Hexa* Document::addHexaQuadsACE (AnaQuads& strquads)
339 Real3 v_bce, v_ade, v_acf, v_bdf; // Sommets
340 Real3 norm_b, norm_d, norm_f; // Normales aux faces
343 Quad* q_a = strquads.tab_quads[0];
344 Quad* q_c = strquads.tab_quads[1];
345 Quad* q_e = strquads.tab_quads[2];
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);
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);
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);
359 e_ce->commonPoint (e_bc, v_bce);
360 e_ae->commonPoint (e_ad, v_ade);
361 e_ac->commonPoint (e_af, v_acf);
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
367 Real3 membre2 = { prod_scalaire (v_bce, norm_b),
368 prod_scalaire (v_ade, norm_d),
369 prod_scalaire (v_acf, norm_f) };
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);
378 printf (" addHexaQuadsACE : systeme impossible\n");
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);
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);
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);
395 Hexa* hexa = new Hexa (q_a, q_b, q_c, q_d, q_e, q_f);
398 // ========================================================= addHexaquadsACD
399 // ==== Construction d'un hexaedre a partir d'un U
400 Hexa* Document::addHexaQuadsACD (AnaQuads& strquads)
403 for (int np=0 ; np<3 && pos_a==NOTHING ; np++)
404 if (strquads.inter_nbre[np]==2)
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];
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
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
423 Edge* e_bc = q_c->getEdge ((nc_ac + 2) MODULO QUAD4);
424 Edge* e_bd = q_d->getEdge ((nd_ad + 2) MODULO QUAD4);
426 Edge* e_ce = q_c->getEdge ((nc_ac + 1) MODULO QUAD4);
427 Edge* e_cf = q_c->getEdge ((nc_ac + 3) MODULO QUAD4);
429 Edge* e_de = q_d->getEdge ((nd_ad + 1) MODULO QUAD4);
430 Edge* e_df = q_d->getEdge ((nd_ad + 3) MODULO QUAD4);
432 Vertex* v_ace = e_ae->commonVertex (e_ce);
433 Vertex* v_ade = e_ae->commonVertex (e_de);
436 permuter_edges (e_ce, e_cf);
437 v_ace = e_ae->commonVertex (e_ce);
442 permuter_edges (e_de, e_df);
443 v_ade = e_ae->commonVertex (e_de);
446 Vertex* v_acf = e_af->commonVertex (e_cf);
447 Vertex* v_adf = e_af->commonVertex (e_df);
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);
454 Edge* e_be = new Edge (v_bce, v_bde);
455 Edge* e_bf = new Edge (v_bcf, v_bdf);
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);
461 Hexa* hexa = new Hexa (q_a, q_b, q_c, q_d, q_e, q_f);
464 // ========================================================= addHexaquadsABCD
465 Hexa* Document::addHexaQuadsABCD (AnaQuads& strquads)
472 for (int np=1 ; np<4 ; np++)
474 if (strquads.inter_edge [pos_a] [np] == NOTHING )
476 else if (pos_c==NOTHING)
482 if (pos_b==NOTHING || pos_c==NOTHING || pos_d==NOTHING)
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];
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
495 Edge* e_ae = q_a->getEdge ((na_ac + 1) MODULO QUAD4);
496 Edge* e_af = q_a->getEdge ((na_ac + 3) MODULO QUAD4);
498 Edge* e_ce = q_c->getEdge ((nc_ac + 1) MODULO QUAD4);
499 Edge* e_cf = q_c->getEdge ((nc_ac + 3) MODULO QUAD4);
501 Edge* e_de = q_d->getEdge ((nd_ad + 1) MODULO QUAD4);
502 Edge* e_df = q_d->getEdge ((nd_ad + 3) MODULO QUAD4);
504 Edge* e_be = q_b->getEdge ((nb_bc + 1) MODULO QUAD4);
505 Edge* e_bf = q_b->getEdge ((nb_bc + 3) MODULO QUAD4);
524 if (e_ae->commonVertex (e_ce) == NULL)
525 PermuterEdges (e_ce, e_cf);
527 if (e_ae->commonVertex (e_de) == NULL)
528 PermuterEdges (e_de, e_df);
530 if (e_ce->commonVertex (e_be) == NULL)
531 PermuterEdges (e_be, e_bf);
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);
536 Hexa* hexa = new Hexa (q_a, q_b, q_c, q_d, q_e, q_f);
539 // ========================================================= addHexaquadsACDE
540 Hexa* Document::addHexaQuadsACDE (AnaQuads& strquads)
547 for (int np=0 ; np<4 ; np++)
548 if (strquads.inter_nbre[np]==3)
550 if (pos_a == NOTHING) pos_a = np;
553 else if (strquads.inter_nbre[np]==2)
555 if (pos_c == NOTHING) pos_c = np;
559 if (pos_a==NOTHING || pos_c==NOTHING || pos_d==NOTHING || pos_e==NOTHING)
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];
567 int na_ac = strquads.inter_edge[pos_a][pos_c]; // Nro dans q_a de e_ac
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
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
575 int ne_ae = strquads.inter_edge[pos_e][pos_a]; // Nro dans q_e de e_ae
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);
584 Vertex* v_bcf = e_cf->opposedVertex (e_cf->commonVertex (e_af));
585 Vertex* v_bdf = e_df->opposedVertex (e_df->commonVertex (e_af));
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);
591 Hexa* hexa = new Hexa (q_a, q_b, q_c, q_d, q_e, q_f);
594 // ========================================================= replaceHexa
595 Elements* Document::replaceHexa (Quads pattern, Vertex* p1, Vertex* c1,
596 Vertex* p2, Vertex* c2, Vertex* p3, Vertex* c3)
598 DumpStart ("replaceHexa", pattern << p1 << c1 << p2 << c2 << p3 << c3);
600 Elements* t_hexas = new Elements (this);
601 Mess << " **** This syntax is deprecated" ;
602 t_hexas->setError (HERR);
604 DumpReturn (t_hexas);
607 // ========================================================= replace
608 Elements* Document::replace (Quads motif, Quads cible, Vertex* p1, Vertex* c1, Vertex* p2, Vertex* c2)
610 DumpStart ("replace", motif << cible << p1 << c1 << p2 << c2);
612 Elements* t_hexas = new Elements (this);
613 int ier = t_hexas->replaceHexas (motif, cible, p1, c1, p2, c2);
617 Mess << " **** Error in Document::replace" ;
618 t_hexas->setError (HERR);
621 DumpReturn (t_hexas);
624 // ========================================================= print_replace
625 void print_replace (Edge* zig, Edge* zag)
627 cout << zig->getName() << " = (" << zig->getVertex(0)->getName()
628 << ", " << zig->getVertex(1)->getName() << ") est clone en ";
630 cout << zag->getName() << " = (" << zag->getVertex(0)->getName()
631 << ", " << zag->getVertex(1)->getName() << ")" << endl;
633 // ========================================================= only_in_hexas
634 bool only_in_hexas (Hexas& thexas, Quad* quad)
636 int nbhexas = thexas.size();
637 int nbp = quad->getNbrParents();
638 for (int nh=0 ; nh <nbp ; nh++)
641 Hexa* hexa = quad->getParent (nh);
642 for (int nc=0 ; pasla && nc < nbhexas ; nc++)
643 pasla = hexa != thexas [nc];
649 // ========================================================= only_in_hexas
650 bool only_in_hexas (Hexas& thexas, Edge* edge)
653 int nbp = edge->getNbrParents();
654 for (int nq=0 ; nq <nbp ; nq++)
656 Quad* quad = edge->getParent (nq);
657 if (NOT only_in_hexas (thexas, quad))
659 cout << " ... inMoreHexas " << edge->makeDefinition() << endl;
663 cout << " ... only_in_hexas " << edge->makeDefinition() << endl;
666 // ========================================================= replace_vertex
667 void replace_vertex (Hexas& thexas, Vertex* node, Vertex* par)
669 int nbh = thexas.size();
670 for (int nh=0 ; nh <nbh ; nh++)
671 thexas[nh]->replaceVertex (node, par);
673 // ========================================================= disconnectEdges
674 Elements* Document::disconnectEdges (Hexas thexas, Edges tedges)
676 DumpStart ("disconnectEdges", thexas << tedges);
679 cout << " +++ Disconnect Edges" << endl;
681 Elements* grid = new Elements (this);
683 grid->checkDisco (thexas, tedges);
685 int nbedges = tedges.size();
686 int nbhexas = thexas.size();
688 if (nbhexas != nbedges)
690 cout << " **** Error in Document::disconnectEdges\n" << endl;
691 cout << " **** Number of Edges and number of Hexas are different\n"
698 thexas[0]->disconnectEdge (tedges[0], grid);
702 for (int nro=0 ; nro<nbedges ; nro++)
704 if (BadElement (tedges[nro]))
706 cout << " **** Eddge number " << nro+1 << " is incorrect"
710 if (BadElement (thexas[nro]))
712 cout << " **** Hexa number " << nro+1 << " is incorrect"
717 cout << nro+1 << " hexa = " << thexas[nro]->getName ()
718 << ", edge = " << tedges[nro]->getName ()
719 << " = (" << tedges[nro]->getVertex(0)->getName ()
720 << ", " << tedges[nro]->getVertex(1)->getName ()
724 for (int nro=0 ; nro<nbhexas ; nro++)
726 int ned = thexas[nro]->findEdge (tedges[nro]);
729 cout << " **** Edge number " << nro+1
730 << " doesnt belong to correspondant hexa" << endl;
735 vector <Vertex*> tvertex (nbedges+1);
737 for (int nro=1 ; nro<nbedges ; nro++)
739 tvertex[nro] = tedges[nro]->commonVertex (tedges[nro-1]);
740 if (tvertex[nro]==NULL)
742 cout << " **** Edge number " << nro
743 << " doesnt intesect next edge" << endl;
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);
753 for (int nro=0 ; nro<nbhexas ; nro++)
755 int ned = thexas[nro]->findEdge (tedges[nro]);
758 cout << " **** Edge number " << nro+1
759 << " doesnt belong to correspondant hexa" << endl;
763 // Fin des controles, on peut y aller ...
765 map <Edge*, int> state_edge;
766 map <Quad*, int> state_quad;
767 enum { UNDEFINED, REPLACED, AS_IS };
769 map <Vertex*, Vertex*> new_vertex;
770 map <Edge*, Edge*> new_edge;
771 map <Quad*, Quad*> new_quad;
773 map <Vertex*, Vertex*> :: iterator it_vertex;
774 map <Edge*, Edge*> :: iterator it_edge;
775 map <Quad*, Quad*> :: iterator it_quad;
777 #define VertexIsNew(v) (it_vertex=new_vertex.find(v))!=new_vertex.end()
779 Vertex* node1 = NULL;
781 for (int nro=0 ; nro<=nbedges ; nro++)
783 Vertex* node0 = node1;
784 node1 = new Vertex (tvertex[nro]);
785 grid->addVertex (node1);
786 new_vertex [tvertex[nro]] = node1;
789 cout << nro << " : " << tvertex[nro]->getName()
790 << " est clone en " << node1->getName() << endl;
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;
800 print_replace (tedges[nro-1], edge);
805 cout << "_____________________________ Autres substitutions" << endl;
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
811 for (int nro=0 ; nro<nbhexas ; nro++)
813 Hexa* hexa = thexas [nro];
814 for (int nro=0 ; nro<HE_MAXI ; nro++)
816 Edge* edge = hexa->getEdge (nro);
817 if (state_edge[edge]==UNDEFINED)
819 Vertex* v1 = edge->getVertex (V_AMONT);
820 Vertex* v2 = edge->getVertex (V_AVAL);
822 if (VertexIsNew (v1))
824 if (only_in_hexas (thexas, edge))
826 replace_vertex (thexas, v1, new_vertex[v1]);
830 v1 = new_vertex [v1];
832 else if (VertexIsNew (v2))
834 if (only_in_hexas (thexas, edge))
836 replace_vertex (thexas, v2, new_vertex[v2]);
840 v2 = new_vertex [v2];
847 Edge* arete = new Edge (v1, v2);
848 new_edge [edge] = arete;
849 grid->addEdge (arete);
851 print_replace (edge, arete);
853 state_edge [edge] = etat;
858 // Un quad non remplace, qui contient un edge remplace
859 // commun a plus de 2 Hexas
860 // doit etre duplique
862 for (int nro=0 ; nro<nbhexas ; nro++)
864 Hexa* hexa = thexas [nro];
865 for (int nro=0 ; nro<HQ_MAXI ; nro++)
867 Quad* quad = hexa->getQuad (nro);
868 if (state_quad[quad]==UNDEFINED)
872 for (int ned=0 ; ned < QUAD4 ; ned++)
874 Edge* edge = quad->getEdge (ned);
875 if (state_edge [edge]==AS_IS)
879 ted[ned] = new_edge[edge];
885 Quad* face = new Quad (ted[0], ted[1], ted[2], ted[3]);
886 new_quad [quad] = face;
887 grid->addQuad (face);
889 state_quad [quad] = etat;
894 for (int nro=0 ; nro<nbhexas ; nro++)
896 Hexa* hexa = thexas [nro];
897 for (it_quad=new_quad.begin(); it_quad != new_quad.end() ; ++it_quad)
899 Quad* pile = it_quad->first;
900 Quad* face = it_quad->second;
901 hexa->replaceQuad (pile, face);
904 for (it_edge=new_edge.begin(); it_edge != new_edge.end() ; ++it_edge)
906 Edge* zig = it_edge->first;
907 Edge* zag = it_edge->second;
908 hexa->replaceEdge (zig, zag);
911 for (it_vertex=new_vertex.begin();
912 it_vertex != new_vertex.end() ; ++it_vertex)
914 Vertex* flip = it_vertex->first;
915 Vertex* flop = it_vertex->second;
916 hexa->replaceVertex (flip, flop);
920 Real3 center0, center1;
924 for (int nro=0 ; nro<=nbhexas ; nro++)
929 hexa1 = thexas [nro];
930 hexa1->getCenter (center1);
932 else if (nro==nbhexas)
934 hexa1->getCenter (center1);
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;
945 Vertex* node = grid->getVertex (nro);
946 matrix.defScale (center1, 0.55);
947 matrix.perform (node);
953 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
954 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
955 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
956 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
957 // . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
958 // ========================================================= prod_vectoriel
959 double* prod_vectoriel (Edge* e1, Edge* e2, double prod[])
961 prod [dir_x] = prod [dir_y] = prod [dir_z] = 0;
962 if (e1==NULL || e2==NULL)
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];
975 // ========================================================= permuter_edges
976 void permuter_edges (Edge* &e1, Edge* &e2, cpchar nm1, cpchar nm2)
980 printf (" ... permuter_edges %s = %s et %s = %s\n",
981 nm1, e1->getName(), nm2, e2->getName() );