2 // C++ : Table d'hexaedres (Evol Versions 3)
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 "HexElements.hxx"
25 #include "HexVector.hxx"
26 #include "HexVertex.hxx"
27 #include "HexEdge.hxx"
28 #include "HexDiagnostics.hxx"
29 #include "HexDocument.hxx"
30 #include "HexHexa.hxx"
37 static bool db = false;
39 // --------------------------------------------------------
45 // --------------------------------------------------------
51 // --------------------------------------------------------
59 // --------------------------------------------------------
62 friend class Elements;
64 int initialize (Vertex* v1, Vertex* v2, Vertex* v3);
65 int addQuad (Quad* quad);
66 int verify (int &nbed, int &nbver);
69 int addEdge (Edge* edge);
70 int addVertex (Vertex* vertex);
73 enum EnumProj { ProjXY, ProjYZ, ProjZX };
75 vector <PatVertex> pat_vertex;
76 vector <PatEdge > pat_edge;
77 vector <PatQuad > pat_quad;
79 int nbr_vertex, nbr_edges, nbr_quads;
82 Real3 base_i, base_j, origine;
85 int pos_edge3, pos_edge4;
88 // ====================================================== initialize
89 int Pattern::initialize (Vertex* v1, Vertex* v2, Vertex* v3)
91 nbr_vertex = nbr_edges = nbr_quads = 0;
94 base_i[0] = base_i[1] = base_i[2] = 0;
95 base_j[0] = base_j[1] = base_j[2] = 0;
97 if (v1==NULL || v2==NULL || v3==NULL)
106 Document* doc = v1->dad ();
107 Edge* edc1 = doc->findEdge (v1, v2);
108 Edge* edc2 = doc->findEdge (v2, v3);
109 if (edc1==NULL || edc2==NULL)
117 pat_edge[0].nbr_refer = pat_edge[1].nbr_refer = 0;
119 pat_vertex[0].v_x = pat_vertex[1].v_x = 0;
120 pat_vertex[2].v_x = 1;
121 pat_vertex[0].v_y = 1;
122 pat_vertex[1].v_y = pat_vertex[2].v_y = 0;
126 v2->getPoint (origine);
129 calc_vecteur (origine, pb, base_i);
130 calc_vecteur (origine, pc, base_j);
137 PutData (determinant);
140 /* ******************************
143 * vx = pu.vI[0] + pv.vJ[0] (1) *+vJ[1] ) };
145 * (1 & 2 ) : pu = (vx.vJ[1] - vy.vJ[0]) / detxy
146 * pv = -(vx.vI[1] - vy.vI[0]) / detxy
148 * (2 & 3 ) : pu = (vy.vJ[2] - vz.vJ[1]) / detyz
149 * pv = -(vy.vI[2] - vz.vI[1]) / detyz
151 * (3 & 1 ) : pu = (vz.vJ[0] - vx.vJ[2]) / detzx
152 * pv = -(vz.vI[0] - vx.vI[2]) / detzx
154 * Les 3 systemes d'equations sont valides.
155 * On va choisir celui dont la valeur absolue du determinant est maximale
156 ****************************** */
158 double detxy = base_i[dir_x]*base_j[dir_y] - base_i[dir_y]*base_j[dir_x];
159 double detyz = base_i[dir_y]*base_j[dir_z] - base_i[dir_z]*base_j[dir_y];
160 double detzx = base_i[dir_z]*base_j[dir_x] - base_i[dir_x]*base_j[dir_z];
165 if (fabs (detyz) > fabs (determinant))
170 if (fabs (detzx) > fabs (determinant))
178 // ====================================================== verify
179 int Pattern::verify (int &nbed, int &nbver)
183 pos_edge3 = pos_edge4 = pos_vertex4 = NOTHING;
185 if (pat_edge[0].nbr_refer!=1 || pat_edge[1].nbr_refer!=1)
188 for (int nro=2 ; nro<nbr_edges; nro++)
190 if (pat_edge[nro].nbr_refer==1)
192 int pv1 = pat_edge[nro].v_amont;
193 int pv2 = pat_edge[nro].v_aval;
198 if (pos_vertex4 == NOTHING)
200 else if (pos_vertex4 != pv2)
208 if (pos_vertex4 == NOTHING)
210 else if (pos_vertex4 != pv1)
218 if (pos_vertex4 == NOTHING)
220 else if (pos_vertex4 != pv2)
228 if (pos_vertex4 == NOTHING)
230 else if (pos_vertex4 != pv1)
238 if (pos_edge3==NOTHING || pos_edge4==NOTHING)
243 // ====================================================== addQuad
244 int Pattern::addQuad (Quad* elt)
251 for (int nro=0; nro<QUAD4 ; nro++)
253 Edge* edge = elt->getEdge (nro);
254 quad.q_edge [nro] = addEdge (edge);
257 pat_quad.push_back (quad);
261 // ====================================================== addEdge
262 int Pattern::addEdge (Edge* elt)
264 for (int nro=0; nro<nbr_edges ; nro++)
266 if (elt==pat_edge [nro].refer)
268 pat_edge [nro].nbr_refer++;
276 edge.v_amont = addVertex (elt->getVertex (V_AMONT));
277 edge.v_aval = addVertex (elt->getVertex (V_AVAL));
279 pat_edge.push_back (edge);
283 // ====================================================== addVertex
284 int Pattern::addVertex (Vertex* elt)
286 for (int nro=0; nro<nbr_vertex ; nro++)
288 if (elt==pat_vertex [nro].refer)
294 double vx = elt->getX() - origine [dir_x];
295 double vy = elt->getY() - origine [dir_y];
296 double vz = elt->getZ() - origine [dir_z];
299 case ProjXY : default :
300 vertex.v_x = (vx*base_j[dir_y] - vy*base_j[dir_x]) / determinant;
301 vertex.v_y = -(vx*base_i[dir_y] - vy*base_i[dir_x]) / determinant;
305 vertex.v_y = (vy*base_j[dir_z] - vz*base_j[dir_y]) / determinant;
306 vertex.v_x = -(vy*base_i[dir_z] - vz*base_i[dir_y]) / determinant;
310 vertex.v_x = (vz*base_j[dir_x] - vx*base_j[dir_z]) / determinant;
311 vertex.v_y = -(vz*base_i[dir_x] - vx*base_i[dir_z]) / determinant;
316 printf (" Vertex nro %d : (%g, %g, %g)\t -> (%g, %g)\n",
317 nbr_vertex, vx, vy, vz, vertex.v_x, vertex.v_y);
318 pat_vertex.push_back (vertex);
322 // -------------------------------------------------------------------
323 // -------------------------------------------------------------------
324 // -------------------------------------------------------------------
325 // -------------------------------------------------------------------
326 // ====================================================== replaceHexas
327 int Elements::replaceHexas (Quads& liste, Vertex* p1, Vertex* c1,
328 Vertex* p2, Vertex* c2, Vertex* p3, Vertex* c3)
330 Edge* edp1 = el_root->findEdge (p1, p2);
331 Edge* edp2 = el_root->findEdge (p2, p3);
332 Edge* edc1 = el_root->findEdge (c1, c2);
333 Edge* edc2 = el_root->findEdge (c2, c3);
335 Quad* quadc = el_root->findQuad (c1, c3);
337 if (edp1==NULL || edp2==NULL || edc1==NULL || edc2==NULL || quadc==NULL)
339 printf ("... Error in HexaBlock function \n");
340 printf ("... doc.replace (lquads, p1,c1, p2,c2, p3,c3)\n");
342 printf ("Vertices p1 and p2 don't define an edge\n");
344 printf ("Vertices p2 and p3 don't define an edge\n");
346 printf ("Vertices c1 and c2 don't define an edge\n");
348 printf ("Vertices c2 and c3 don't define an edge\n");
349 else if (quadc==NULL)
350 printf ("Vertices c1 and c3 don't define a quad\n");
354 int np = quadc->getNbrParents ();
355 Hexa* hexac = quadc->getParent (0);
357 if (np!=1 || hexac==NULL)
359 printf ("... Error in HexaBlock function \n");
360 printf ("... doc.replace (lquads, p1,c1, p2,c2, p3,c3)\n");
361 printf ("Quad (c1,c2,c3) is not an external quad\n");
364 // Analyse du pattern
365 int nbquads = liste.size();
367 int ier = pattern.initialize (p1, p2, p3);
370 printf ("... Error in HexaBlock function \n");
371 printf ("... doc.replace (lquads, p1,c1, p2,c2, p3,c3)\n");
372 printf ("Vertices (p1,p2,p3) don't define a virtual quad\n");
376 for (int nq=0 ; nq<nbquads ; nq++)
378 ier = pattern.addQuad (liste[nq]);
381 printf ("... Error in HexaBlock function \n");
382 printf ("... doc.replace (lquads, p1,c1, p2,c2, p3,c3)\n");
383 printf ("Quad nr %d of the list is NULL\n", nq+1);
390 ier = pattern.verify (nbed, nbv);
393 printf ("... Error in HexaBlock function \n");
394 printf ("... doc.replace (lquads, p1,c1, p2,c2, p3,c3)\n");
395 printf ("Vertices (p1,p2,p3) don't define a virtual quad\n");
398 // Analyse de la cible
402 pil_hexa.push_back (hexac);
403 pil_quad.push_back (quadc);
405 // Constitution de la pile des hexaedres
409 Quad* oppos = hexac->getOpposedQuad (quadc);
410 pil_quad.push_back (oppos);
411 more = oppos->getNbrParents() == 2;
415 if (oppos->getParent (0)==hexac)
416 hexac = oppos->getParent(1);
418 hexac = oppos->getParent(0);
419 pil_hexa.push_back (hexac);
424 resize (GR_REPLACE, nbquads, nbh, nbv, nbed);
426 // 1) Constitution des 4 coins de la cible
427 // 2) Quadriller la face du haut
428 // 3) Trouver les 4 coins de la face du bas
429 // 4) Quadriller la face du bas
430 // 6) Decouper l'hexaedre defini par ces 2 faces
431 // 7) S'il y a un hexaedre dessous goto 3)
433 Vertex* tvert[QUAD4] = { c1, c2, c3, pil_quad[0]->getOpposVertex(c2) };
434 replaceQuad (1, &pattern, pil_quad[0], tvert);
435 for (int nh=1 ; nh<gr_hauteur ; nh++)
437 for (int nro=0 ; nro<QUAD4 ; nro++)
439 Edge* edv = pil_hexa[nh-1]->getPerpendicularEdge (pil_quad[nh-1],
441 tvert [nro] = edv ->opposedVertex(tvert[nro]);
443 replaceQuad (nh+1, &pattern, pil_quad[nh], tvert);
444 replaceHexa (nh, &pattern, pil_hexa[nh-1]);
447 for (int nh=0 ; nh<=nbh ; nh++)
448 pil_quad[nh]->remove ();
450 extrudeQuad (&pattern);
451 replaceHexa (0, &pattern, NULL);
454 // ====================================================== repVertex
455 void Elements::repVertex (int nh, int nro, Vertex* elt)
457 int addr = nh * pat_nbvertex + nro;
458 if (tab_vertex[addr] == NULL)
459 tab_vertex[addr] = elt;
461 // ====================================================== repVertex
462 Vertex* Elements::repVertex (int nh, int nro, double px, double py, double pz)
464 int addr = nh * pat_nbvertex + nro;
466 if (tab_vertex[addr] == NULL)
467 tab_vertex[addr] = el_root->addVertex (px, py, pz);
469 return tab_vertex[addr];
471 // ====================================================== repVertex
472 Vertex* Elements::repVertex (int nh, int nro)
474 int addr = nh * pat_nbvertex + nro;
475 return tab_vertex[addr];
477 // ====================================================== repEdgeH
478 void Elements::repEdgeH (int nh, int nro, Edge* elt)
480 int addr = nh * (pat_nbedges + pat_nbvertex) + nro;
481 if (tab_edge[addr] == NULL)
482 tab_edge[addr] = elt;
484 // ====================================================== repEdgeH
485 Edge* Elements::repEdgeH (int nh, int nro, Vertex* v1, Vertex* v2)
487 int addr = nh * (pat_nbedges + pat_nbvertex) + nro;
488 if (tab_edge[addr] == NULL)
489 tab_edge[addr] = el_root->addEdge (v1, v2);
490 return tab_edge[addr];
492 // ====================================================== repEdgeH
493 Edge* Elements::repEdgeH (int nh, int nro)
495 int addr = nh * (pat_nbedges + pat_nbvertex) + nro;
496 return tab_edge[addr];
498 // ====================================================== repEdgeV
499 void Elements::repEdgeV (int nh, int nro, Edge* elt)
501 int addr = nh * (pat_nbedges + pat_nbvertex) + pat_nbedges + nro;
502 if (tab_edge[addr] == NULL)
503 tab_edge[addr] = elt;
505 // ====================================================== repEdgeV
506 Edge* Elements::repEdgeV (int nh, int nro, Vertex* v1, Vertex* v2)
508 int addr = nh * (pat_nbedges + pat_nbvertex) + pat_nbedges + nro;
509 if (tab_edge[addr] == NULL)
510 tab_edge[addr] = el_root->addEdge (v1, v2);
511 return tab_edge[addr];
513 // ====================================================== repEdgeV
514 Edge* Elements::repEdgeV (int nh, int nro)
516 int addr = nh * (pat_nbedges + pat_nbvertex) + pat_nbedges + nro;
517 return tab_edge[addr];
519 // ====================================================== repQuadH
520 void Elements::repQuadH (int nh, int nro, Quad* elt)
522 int addr = nh * (nbr_orig + pat_nbedges) + nro;
523 if (tab_quad[addr] == NULL)
524 tab_quad[addr] = elt;
526 // ====================================================== repQuadH
527 Quad* Elements::repQuadH (int nh, int nro)
529 int addr = nh * (nbr_orig + pat_nbedges) + nro;
530 return tab_quad [addr];
532 // ====================================================== repQuadH
533 Quad* Elements::repQuadH (int nh, int nro, Edge* e1, Edge* e2, Edge* e3,
536 int addr = nh * (nbr_orig + pat_nbedges) + nro;
537 if (tab_quad[addr] == NULL)
538 tab_quad[addr] = el_root->addQuad (e1, e2, e3, e4);
539 return tab_quad [addr];
541 // ====================================================== repQuadV
542 void Elements::repQuadV (int nh, int nro, Quad* elt)
544 int addr = nh * (nbr_orig + pat_nbedges) + nbr_orig + nro;
545 if (tab_quad[addr] == NULL)
546 tab_quad[addr] = elt;
548 // ====================================================== repQuadV
549 Quad* Elements::repQuadV (int nh, int nro)
551 int addr = nh * (nbr_orig + pat_nbedges) + nbr_orig + nro;
552 return tab_quad [addr];
554 // ====================================================== repQuadV
555 Quad* Elements::repQuadV (int nh, int nro, Edge* e1, Edge* e2, Edge* e3,
558 int addr = nh * (nbr_orig + pat_nbedges) + nbr_orig + nro;
559 if (tab_quad[addr] == NULL)
560 tab_quad[addr] = el_root->addQuad (e1, e2, e3, e4);
561 return tab_quad [addr];
563 // ====================================================== repHexa
564 Hexa* Elements::repHexa (int nh, int nro, Quad* qa, Quad* qb, Quad* qc,
565 Quad* qd, Quad* qe, Quad* qf)
567 int addr = nh * nbr_orig + nro;
568 if (tab_hexa[addr] == NULL)
569 tab_hexa[addr] = el_root->addHexa (qa, qb, qc, qd, qe, qf);
570 return tab_hexa [addr];
572 // ====================================================== replaceQuad
573 // ==== Creation des quads horizontaux
574 int Elements::replaceQuad (int nh, Pattern* pat, Quad* quad, Vertex* tvert[])
576 Vertex* pvert[QUAD4] = { tvert[0], tvert[1], tvert[2], tvert[3] };
578 int vnro [4] = { 0, 1, 2, pat->pos_vertex4};
579 int ednro[4] = { 0, 1, pat->pos_edge3, pat->pos_edge4 };
581 // Enregistrement des vertex & edges existant
582 for (int nro=0 ; nro<QUAD4 ; nro++)
584 Vertex* vh = pvert [nro];
585 Vertex* vh1 = pvert [(nro+1) MODULO QUAD4];
586 Edge* edh = quad->findEdge (vh, vh1);
588 repVertex (nh, vnro [nro], vh);
589 repEdgeH (nh, ednro [nro], edh);
592 Real3 orig, ib, jb, pb, pc;
593 // Creation des vertex
594 tvert[0]->getPoint (pc);
595 tvert[1]->getPoint (orig);
596 tvert[2]->getPoint (pb);
597 calc_vecteur (orig, pb, ib);
598 calc_vecteur (orig, pc, jb);
600 for (int nro=0 ; nro<pat->nbr_vertex ; nro++)
602 double lambda = pat->pat_vertex [nro].v_x;
603 double mu = pat->pat_vertex [nro].v_y;
604 double px = orig[dir_x] + lambda*ib[dir_x] + mu*jb[dir_x];
605 double py = orig[dir_y] + lambda*ib[dir_y] + mu*jb[dir_y];
606 double pz = orig[dir_z] + lambda*ib[dir_z] + mu*jb[dir_z];
607 repVertex (nh, nro, px, py, pz);
609 // Creation des edges horizontaux
610 for (int nro=0 ; nro<pat->nbr_edges ; nro++)
612 int nv1 = pat->pat_edge [nro].v_amont;
613 int nv2 = pat->pat_edge [nro].v_aval;
614 Vertex* v1 = repVertex (nh, nv1);
615 Vertex* v2 = repVertex (nh, nv2);
616 repEdgeH (nh, nro, v1, v2);
618 // Creation des quads horizontaux
619 for (int nro=0 ; nro<pat->nbr_quads ; nro++)
621 Edge* eda = repEdgeH (nh, pat->pat_quad [nro].q_edge [0]);
622 Edge* edb = repEdgeH (nh, pat->pat_quad [nro].q_edge [1]);
623 Edge* edc = repEdgeH (nh, pat->pat_quad [nro].q_edge [2]);
624 Edge* edd = repEdgeH (nh, pat->pat_quad [nro].q_edge [3]);
626 repQuadH (nh, nro, eda, edb, edc, edd);
630 // ====================================================== extrudeQuad
631 // ==== Creation des quads horizontaux
632 int Elements::extrudeQuad (Pattern* pat)
634 // Creation des vertex de niveau 0
635 for (int nro=0 ; nro<pat->nbr_vertex ; nro++)
637 Vertex* v1 = repVertex (1, nro);
638 Vertex* v2 = repVertex (2, nro);
639 double px = 2*v1->getX() - v2->getX();
640 double py = 2*v1->getY() - v2->getY();
641 double pz = 2*v1->getZ() - v2->getZ();
642 repVertex (0, nro, px, py, pz);
644 // Creation des edges horizontaux
645 for (int nro=0 ; nro<pat->nbr_edges ; nro++)
647 int nv1 = pat->pat_edge [nro].v_amont;
648 int nv2 = pat->pat_edge [nro].v_aval;
649 Vertex* v1 = repVertex (0, nv1);
650 Vertex* v2 = repVertex (0, nv2);
651 repEdgeH (0, nro, v1, v2);
653 // Creation des quads horizontaux
654 for (int nro=0 ; nro<pat->nbr_quads ; nro++)
656 Edge* eda = repEdgeH (0, pat->pat_quad [nro].q_edge [0]);
657 Edge* edb = repEdgeH (0, pat->pat_quad [nro].q_edge [1]);
658 Edge* edc = repEdgeH (0, pat->pat_quad [nro].q_edge [2]);
659 Edge* edd = repEdgeH (0, pat->pat_quad [nro].q_edge [3]);
661 repQuadH (0, nro, eda, edb, edc, edd);
665 // ====================================================== replaceHexa
666 int Elements::replaceHexa (int nh, Pattern* pat, Hexa* hexa)
668 int vnro [4] = { 0, 1, 2, pat->pos_vertex4};
669 int ednro[4] = { 0, 1, pat->pos_edge3, pat->pos_edge4 };
671 // Enregistrement des edges & quads existants
674 for (int nro=0 ; nro<QUAD4 ; nro++)
676 Vertex* vh = repVertex (nh, vnro [nro]);
677 Vertex* vh1 = repVertex (nh+1, vnro [nro]);
678 Edge* edv = hexa->findEdge (vh, vh1);
679 repEdgeV (nh, vnro [nro], edv);
681 Edge* edh = repEdgeH (nh, ednro [nro]);
682 Edge* edh1 = repEdgeH (nh+1, ednro [nro]);
683 Quad* quadv = hexa->findQuad (edh, edh1);
684 repQuadV (nh, ednro [nro], quadv);
687 // Creation des edges verticaux
688 for (int nro=0 ; nro<pat->nbr_vertex ; nro++)
690 Vertex* v1 = repVertex (nh, nro);
691 Vertex* v2 = repVertex (nh+1, nro);
692 repEdgeV (nh, nro, v1, v2);
694 // Creation des quads verticaux
695 for (int nro=0 ; nro<pat->nbr_edges ; nro++)
697 int nv1 = pat->pat_edge [nro].v_amont;
698 int nv2 = pat->pat_edge [nro].v_aval;
700 Edge* eda = repEdgeH (nh, nro);
701 Edge* edb = repEdgeV (nh, nv1);
702 Edge* edc = repEdgeH (nh+1, nro);
703 Edge* edd = repEdgeV (nh, nv2);
704 repQuadV (nh, nro, eda, edb, edc, edd);
706 // Creation des hexaedres
707 for (int nro=0 ; nro<pat->nbr_quads ; nro++)
709 Quad* qa = repQuadH (nh, nro);
710 Quad* qb = repQuadH (nh+1, nro);
711 Quad* qc = repQuadV (nh, pat->pat_quad [nro].q_edge [0]);
712 Quad* qd = repQuadV (nh, pat->pat_quad [nro].q_edge [2]);
713 Quad* qe = repQuadV (nh, pat->pat_quad [nro].q_edge [1]);
714 Quad* qf = repQuadV (nh, pat->pat_quad [nro].q_edge [3]);
715 repHexa (nh, nro, qa, qb, qc, qd, qe, qf);