2 // C++ : Table d'hexaedres (Evol Versions 3)
4 // Copyright (C) 2009-2012 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"
36 static bool db = false;
38 // --------------------------------------------------------
44 // --------------------------------------------------------
50 // --------------------------------------------------------
58 // --------------------------------------------------------
61 friend class Elements;
63 int initialize (Vertex* v1, Vertex* v2, Vertex* v3);
64 int addQuad (Quad* quad);
65 int verify (int &nbed, int &nbver);
68 int addEdge (Edge* edge);
69 int addVertex (Vertex* vertex);
72 enum EnumProj { ProjXY, ProjYZ, ProjZX };
74 vector <PatVertex> pat_vertex;
75 vector <PatEdge > pat_edge;
76 vector <PatQuad > pat_quad;
78 int nbr_vertex, nbr_edges, nbr_quads;
81 Real3 base_i, base_j, origine;
84 int pos_edge3, pos_edge4;
87 // ====================================================== initialize
88 int Pattern::initialize (Vertex* v1, Vertex* v2, Vertex* v3)
90 nbr_vertex = nbr_edges = nbr_quads = 0;
93 base_i[0] = base_i[1] = base_i[2] = 0;
94 base_j[0] = base_j[1] = base_j[2] = 0;
96 if (v1==NULL || v2==NULL || v3==NULL)
105 Document* doc = v1->dad ();
106 Edge* edc1 = doc->findEdge (v1, v2);
107 Edge* edc2 = doc->findEdge (v2, v3);
108 if (edc1==NULL || edc2==NULL)
116 pat_edge[0].nbr_refer = pat_edge[1].nbr_refer = 0;
118 pat_vertex[0].v_x = pat_vertex[1].v_x = 0;
119 pat_vertex[2].v_x = 1;
120 pat_vertex[0].v_y = 1;
121 pat_vertex[1].v_y = pat_vertex[2].v_y = 0;
125 v2->getPoint (origine);
128 calc_vecteur (origine, pb, base_i);
129 calc_vecteur (origine, pc, base_j);
136 PutData (determinant);
139 /* ******************************
142 * vx = pu.vI[0] + pv.vJ[0] (1) *+vJ[1] ) };
144 * (1 & 2 ) : pu = (vx.vJ[1] - vy.vJ[0]) / detxy
145 * pv = -(vx.vI[1] - vy.vI[0]) / detxy
147 * (2 & 3 ) : pu = (vy.vJ[2] - vz.vJ[1]) / detyz
148 * pv = -(vy.vI[2] - vz.vI[1]) / detyz
150 * (3 & 1 ) : pu = (vz.vJ[0] - vx.vJ[2]) / detzx
151 * pv = -(vz.vI[0] - vx.vI[2]) / detzx
153 * Les 3 systemes d'equations sont valides.
154 * On va choisir celui dont la valeur absolue du determinant est maximale
155 ****************************** */
157 double detxy = base_i[dir_x]*base_j[dir_y] - base_i[dir_y]*base_j[dir_x];
158 double detyz = base_i[dir_y]*base_j[dir_z] - base_i[dir_z]*base_j[dir_y];
159 double detzx = base_i[dir_z]*base_j[dir_x] - base_i[dir_x]*base_j[dir_z];
164 if (fabs (detyz) > fabs (determinant))
169 if (fabs (detzx) > fabs (determinant))
177 // ====================================================== verify
178 int Pattern::verify (int &nbed, int &nbver)
182 pos_edge3 = pos_edge4 = pos_vertex4 = NOTHING;
184 if (pat_edge[0].nbr_refer!=1 || pat_edge[1].nbr_refer!=1)
187 for (int nro=2 ; nro<nbr_edges; nro++)
189 if (pat_edge[nro].nbr_refer==1)
191 int pv1 = pat_edge[nro].v_amont;
192 int pv2 = pat_edge[nro].v_aval;
197 if (pos_vertex4 == NOTHING)
199 else if (pos_vertex4 != pv2)
207 if (pos_vertex4 == NOTHING)
209 else if (pos_vertex4 != pv1)
217 if (pos_vertex4 == NOTHING)
219 else if (pos_vertex4 != pv2)
227 if (pos_vertex4 == NOTHING)
229 else if (pos_vertex4 != pv1)
237 if (pos_edge3==NOTHING || pos_edge4==NOTHING)
242 // ====================================================== addQuad
243 int Pattern::addQuad (Quad* elt)
250 for (int nro=0; nro<QUAD4 ; nro++)
252 Edge* edge = elt->getEdge (nro);
253 quad.q_edge [nro] = addEdge (edge);
256 pat_quad.push_back (quad);
260 // ====================================================== addEdge
261 int Pattern::addEdge (Edge* elt)
263 for (int nro=0; nro<nbr_edges ; nro++)
265 if (elt==pat_edge [nro].refer)
267 pat_edge [nro].nbr_refer++;
275 edge.v_amont = addVertex (elt->getVertex (V_AMONT));
276 edge.v_aval = addVertex (elt->getVertex (V_AVAL));
278 pat_edge.push_back (edge);
282 // ====================================================== addVertex
283 int Pattern::addVertex (Vertex* elt)
285 for (int nro=0; nro<nbr_vertex ; nro++)
287 if (elt==pat_vertex [nro].refer)
293 double vx = elt->getX() - origine [dir_x];
294 double vy = elt->getY() - origine [dir_y];
295 double vz = elt->getZ() - origine [dir_z];
298 case ProjXY : default :
299 vertex.v_x = (vx*base_j[dir_y] - vy*base_j[dir_x]) / determinant;
300 vertex.v_y = -(vx*base_i[dir_y] - vy*base_i[dir_x]) / determinant;
304 vertex.v_y = (vy*base_j[dir_z] - vz*base_j[dir_y]) / determinant;
305 vertex.v_x = -(vy*base_i[dir_z] - vz*base_i[dir_y]) / determinant;
309 vertex.v_x = (vz*base_j[dir_x] - vx*base_j[dir_z]) / determinant;
310 vertex.v_y = -(vz*base_i[dir_x] - vx*base_i[dir_z]) / determinant;
315 printf (" Vertex nro %d : (%g, %g, %g)\t -> (%g, %g)\n",
316 nbr_vertex, vx, vy, vz, vertex.v_x, vertex.v_y);
317 pat_vertex.push_back (vertex);
321 // -------------------------------------------------------------------
322 // -------------------------------------------------------------------
323 // -------------------------------------------------------------------
324 // -------------------------------------------------------------------
325 // ====================================================== replaceHexas
326 int Elements::replaceHexas (Quads& liste, Vertex* p1, Vertex* c1,
327 Vertex* p2, Vertex* c2, Vertex* p3, Vertex* c3)
329 Edge* edp1 = el_root->findEdge (p1, p2);
330 Edge* edp2 = el_root->findEdge (p2, p3);
331 Edge* edc1 = el_root->findEdge (c1, c2);
332 Edge* edc2 = el_root->findEdge (c2, c3);
334 Quad* quadc = el_root->findQuad (c1, c3);
336 if (edp1==NULL || edp2==NULL || edc1==NULL || edc2==NULL || quadc==NULL)
338 printf ("... Error in HexaBlock function \n");
339 printf ("... doc.replace (lquads, p1,c1, p2,c2, p3,c3)\n");
341 printf ("Vertices p1 and p2 don't define an edge\n");
343 printf ("Vertices p2 and p3 don't define an edge\n");
345 printf ("Vertices c1 and c2 don't define an edge\n");
347 printf ("Vertices c2 and c3 don't define an edge\n");
348 else if (quadc==NULL)
349 printf ("Vertices c1 and c3 don't define a quad\n");
353 int np = quadc->getNbrParents ();
354 Hexa* hexac = quadc->getParent (0);
356 if (np!=1 || hexac==NULL)
358 printf ("... Error in HexaBlock function \n");
359 printf ("... doc.replace (lquads, p1,c1, p2,c2, p3,c3)\n");
360 printf ("Quad (c1,c2,c3) is not an external quad\n");
363 // Analyse du pattern
364 int nbquads = liste.size();
366 int ier = pattern.initialize (p1, p2, p3);
369 printf ("... Error in HexaBlock function \n");
370 printf ("... doc.replace (lquads, p1,c1, p2,c2, p3,c3)\n");
371 printf ("Vertices (p1,p2,p3) don't define a virtual quad\n");
375 for (int nq=0 ; nq<nbquads ; nq++)
377 ier = pattern.addQuad (liste[nq]);
380 printf ("... Error in HexaBlock function \n");
381 printf ("... doc.replace (lquads, p1,c1, p2,c2, p3,c3)\n");
382 printf ("Quad nr %d of the list is NULL\n", nq+1);
389 ier = pattern.verify (nbed, nbv);
392 printf ("... Error in HexaBlock function \n");
393 printf ("... doc.replace (lquads, p1,c1, p2,c2, p3,c3)\n");
394 printf ("Vertices (p1,p2,p3) don't define a virtual quad\n");
397 // Analyse de la cible
401 pil_hexa.push_back (hexac);
402 pil_quad.push_back (quadc);
404 // Constitution de la pile des hexaedres
408 Quad* oppos = hexac->getOpposedQuad (quadc);
409 pil_quad.push_back (oppos);
410 more = oppos->getNbrParents() == 2;
414 if (oppos->getParent (0)==hexac)
415 hexac = oppos->getParent(1);
417 hexac = oppos->getParent(0);
418 pil_hexa.push_back (hexac);
423 resize (GR_REPLACE, nbquads, nbh, nbv, nbed);
425 // 1) Constitution des 4 coins de la cible
426 // 2) Quadriller la face du haut
427 // 3) Trouver les 4 coins de la face du bas
428 // 4) Quadriller la face du bas
429 // 6) Decouper l'hexaedre defini par ces 2 faces
430 // 7) S'il y a un hexaedre dessous goto 3)
432 Vertex* tvert[QUAD4] = { c1, c2, c3, pil_quad[0]->getOpposVertex(c2) };
433 replaceQuad (1, &pattern, pil_quad[0], tvert);
434 for (int nh=1 ; nh<gr_hauteur ; nh++)
436 for (int nro=0 ; nro<QUAD4 ; nro++)
438 Edge* edv = pil_hexa[nh-1]->getPerpendicularEdge (pil_quad[nh-1],
440 tvert [nro] = edv ->opposedVertex(tvert[nro]);
442 replaceQuad (nh+1, &pattern, pil_quad[nh], tvert);
443 replaceHexa (nh, &pattern, pil_hexa[nh-1]);
446 for (int nh=0 ; nh<=nbh ; nh++)
447 pil_quad[nh]->remove ();
449 extrudeQuad (&pattern);
450 replaceHexa (0, &pattern, NULL);
453 // ====================================================== repVertex
454 void Elements::repVertex (int nh, int nro, Vertex* elt)
456 int addr = nh * pat_nbvertex + nro;
457 if (tab_vertex[addr] == NULL)
458 tab_vertex[addr] = elt;
460 // ====================================================== repVertex
461 Vertex* Elements::repVertex (int nh, int nro, double px, double py, double pz)
463 int addr = nh * pat_nbvertex + nro;
465 if (tab_vertex[addr] == NULL)
466 tab_vertex[addr] = el_root->addVertex (px, py, pz);
468 return tab_vertex[addr];
470 // ====================================================== repVertex
471 Vertex* Elements::repVertex (int nh, int nro)
473 int addr = nh * pat_nbvertex + nro;
474 return tab_vertex[addr];
476 // ====================================================== repEdgeH
477 void Elements::repEdgeH (int nh, int nro, Edge* elt)
479 int addr = nh * (pat_nbedges + pat_nbvertex) + nro;
480 if (tab_edge[addr] == NULL)
481 tab_edge[addr] = elt;
483 // ====================================================== repEdgeH
484 Edge* Elements::repEdgeH (int nh, int nro, Vertex* v1, Vertex* v2)
486 int addr = nh * (pat_nbedges + pat_nbvertex) + nro;
487 if (tab_edge[addr] == NULL)
488 tab_edge[addr] = el_root->addEdge (v1, v2);
489 return tab_edge[addr];
491 // ====================================================== repEdgeH
492 Edge* Elements::repEdgeH (int nh, int nro)
494 int addr = nh * (pat_nbedges + pat_nbvertex) + nro;
495 return tab_edge[addr];
497 // ====================================================== repEdgeV
498 void Elements::repEdgeV (int nh, int nro, Edge* elt)
500 int addr = nh * (pat_nbedges + pat_nbvertex) + pat_nbedges + nro;
501 if (tab_edge[addr] == NULL)
502 tab_edge[addr] = elt;
504 // ====================================================== repEdgeV
505 Edge* Elements::repEdgeV (int nh, int nro, Vertex* v1, Vertex* v2)
507 int addr = nh * (pat_nbedges + pat_nbvertex) + pat_nbedges + nro;
508 if (tab_edge[addr] == NULL)
509 tab_edge[addr] = el_root->addEdge (v1, v2);
510 return tab_edge[addr];
512 // ====================================================== repEdgeV
513 Edge* Elements::repEdgeV (int nh, int nro)
515 int addr = nh * (pat_nbedges + pat_nbvertex) + pat_nbedges + nro;
516 return tab_edge[addr];
518 // ====================================================== repQuadH
519 void Elements::repQuadH (int nh, int nro, Quad* elt)
521 int addr = nh * (nbr_orig + pat_nbedges) + nro;
522 if (tab_quad[addr] == NULL)
523 tab_quad[addr] = elt;
525 // ====================================================== repQuadH
526 Quad* Elements::repQuadH (int nh, int nro)
528 int addr = nh * (nbr_orig + pat_nbedges) + nro;
529 return tab_quad [addr];
531 // ====================================================== repQuadH
532 Quad* Elements::repQuadH (int nh, int nro, Edge* e1, Edge* e2, Edge* e3,
535 int addr = nh * (nbr_orig + pat_nbedges) + nro;
536 if (tab_quad[addr] == NULL)
537 tab_quad[addr] = el_root->addQuad (e1, e2, e3, e4);
538 return tab_quad [addr];
540 // ====================================================== repQuadV
541 void Elements::repQuadV (int nh, int nro, Quad* elt)
543 int addr = nh * (nbr_orig + pat_nbedges) + nbr_orig + nro;
544 if (tab_quad[addr] == NULL)
545 tab_quad[addr] = elt;
547 // ====================================================== repQuadV
548 Quad* Elements::repQuadV (int nh, int nro)
550 int addr = nh * (nbr_orig + pat_nbedges) + nbr_orig + nro;
551 return tab_quad [addr];
553 // ====================================================== repQuadV
554 Quad* Elements::repQuadV (int nh, int nro, Edge* e1, Edge* e2, Edge* e3,
557 int addr = nh * (nbr_orig + pat_nbedges) + nbr_orig + nro;
558 if (tab_quad[addr] == NULL)
559 tab_quad[addr] = el_root->addQuad (e1, e2, e3, e4);
560 return tab_quad [addr];
562 // ====================================================== repHexa
563 Hexa* Elements::repHexa (int nh, int nro, Quad* qa, Quad* qb, Quad* qc,
564 Quad* qd, Quad* qe, Quad* qf)
566 int addr = nh * nbr_orig + nro;
567 if (tab_hexa[addr] == NULL)
568 tab_hexa[addr] = el_root->addHexa (qa, qb, qc, qd, qe, qf);
569 return tab_hexa [addr];
571 // ====================================================== replaceQuad
572 // ==== Creation des quads horizontaux
573 int Elements::replaceQuad (int nh, Pattern* pat, Quad* quad, Vertex* tvert[])
575 Vertex* pvert[QUAD4] = { tvert[0], tvert[1], tvert[2], tvert[3] };
577 int vnro [4] = { 0, 1, 2, pat->pos_vertex4};
578 int ednro[4] = { 0, 1, pat->pos_edge3, pat->pos_edge4 };
580 // Enregistrement des vertex & edges existant
581 for (int nro=0 ; nro<QUAD4 ; nro++)
583 Vertex* vh = pvert [nro];
584 Vertex* vh1 = pvert [(nro+1) MODULO QUAD4];
585 Edge* edh = quad->findEdge (vh, vh1);
587 repVertex (nh, vnro [nro], vh);
588 repEdgeH (nh, ednro [nro], edh);
591 Real3 orig, ib, jb, pb, pc;
592 // Creation des vertex
593 tvert[0]->getPoint (pc);
594 tvert[1]->getPoint (orig);
595 tvert[2]->getPoint (pb);
596 calc_vecteur (orig, pb, ib);
597 calc_vecteur (orig, pc, jb);
599 for (int nro=0 ; nro<pat->nbr_vertex ; nro++)
601 double lambda = pat->pat_vertex [nro].v_x;
602 double mu = pat->pat_vertex [nro].v_y;
603 double px = orig[dir_x] + lambda*ib[dir_x] + mu*jb[dir_x];
604 double py = orig[dir_y] + lambda*ib[dir_y] + mu*jb[dir_y];
605 double pz = orig[dir_z] + lambda*ib[dir_z] + mu*jb[dir_z];
606 repVertex (nh, nro, px, py, pz);
608 // Creation des edges horizontaux
609 for (int nro=0 ; nro<pat->nbr_edges ; nro++)
611 int nv1 = pat->pat_edge [nro].v_amont;
612 int nv2 = pat->pat_edge [nro].v_aval;
613 Vertex* v1 = repVertex (nh, nv1);
614 Vertex* v2 = repVertex (nh, nv2);
615 repEdgeH (nh, nro, v1, v2);
617 // Creation des quads horizontaux
618 for (int nro=0 ; nro<pat->nbr_quads ; nro++)
620 Edge* eda = repEdgeH (nh, pat->pat_quad [nro].q_edge [0]);
621 Edge* edb = repEdgeH (nh, pat->pat_quad [nro].q_edge [1]);
622 Edge* edc = repEdgeH (nh, pat->pat_quad [nro].q_edge [2]);
623 Edge* edd = repEdgeH (nh, pat->pat_quad [nro].q_edge [3]);
625 repQuadH (nh, nro, eda, edb, edc, edd);
629 // ====================================================== extrudeQuad
630 // ==== Creation des quads horizontaux
631 int Elements::extrudeQuad (Pattern* pat)
633 // Creation des vertex de niveau 0
634 for (int nro=0 ; nro<pat->nbr_vertex ; nro++)
636 Vertex* v1 = repVertex (1, nro);
637 Vertex* v2 = repVertex (2, nro);
638 double px = 2*v1->getX() - v2->getX();
639 double py = 2*v1->getY() - v2->getY();
640 double pz = 2*v1->getZ() - v2->getZ();
641 repVertex (0, nro, px, py, pz);
643 // Creation des edges horizontaux
644 for (int nro=0 ; nro<pat->nbr_edges ; nro++)
646 int nv1 = pat->pat_edge [nro].v_amont;
647 int nv2 = pat->pat_edge [nro].v_aval;
648 Vertex* v1 = repVertex (0, nv1);
649 Vertex* v2 = repVertex (0, nv2);
650 repEdgeH (0, nro, v1, v2);
652 // Creation des quads horizontaux
653 for (int nro=0 ; nro<pat->nbr_quads ; nro++)
655 Edge* eda = repEdgeH (0, pat->pat_quad [nro].q_edge [0]);
656 Edge* edb = repEdgeH (0, pat->pat_quad [nro].q_edge [1]);
657 Edge* edc = repEdgeH (0, pat->pat_quad [nro].q_edge [2]);
658 Edge* edd = repEdgeH (0, pat->pat_quad [nro].q_edge [3]);
660 repQuadH (0, nro, eda, edb, edc, edd);
664 // ====================================================== replaceHexa
665 int Elements::replaceHexa (int nh, Pattern* pat, Hexa* hexa)
667 int vnro [4] = { 0, 1, 2, pat->pos_vertex4};
668 int ednro[4] = { 0, 1, pat->pos_edge3, pat->pos_edge4 };
670 // Enregistrement des edges & quads existants
673 for (int nro=0 ; nro<QUAD4 ; nro++)
675 Vertex* vh = repVertex (nh, vnro [nro]);
676 Vertex* vh1 = repVertex (nh+1, vnro [nro]);
677 Edge* edv = hexa->findEdge (vh, vh1);
678 repEdgeV (nh, vnro [nro], edv);
680 Edge* edh = repEdgeH (nh, ednro [nro]);
681 Edge* edh1 = repEdgeH (nh+1, ednro [nro]);
682 Quad* quadv = hexa->findQuad (edh, edh1);
683 repQuadV (nh, ednro [nro], quadv);
686 // Creation des edges verticaux
687 for (int nro=0 ; nro<pat->nbr_vertex ; nro++)
689 Vertex* v1 = repVertex (nh, nro);
690 Vertex* v2 = repVertex (nh+1, nro);
691 repEdgeV (nh, nro, v1, v2);
693 // Creation des quads verticaux
694 for (int nro=0 ; nro<pat->nbr_edges ; nro++)
696 int nv1 = pat->pat_edge [nro].v_amont;
697 int nv2 = pat->pat_edge [nro].v_aval;
699 Edge* eda = repEdgeH (nh, nro);
700 Edge* edb = repEdgeV (nh, nv1);
701 Edge* edc = repEdgeH (nh+1, nro);
702 Edge* edd = repEdgeV (nh, nv2);
703 repQuadV (nh, nro, eda, edb, edc, edd);
705 // Creation des hexaedres
706 for (int nro=0 ; nro<pat->nbr_quads ; nro++)
708 Quad* qa = repQuadH (nh, nro);
709 Quad* qb = repQuadH (nh+1, nro);
710 Quad* qc = repQuadV (nh, pat->pat_quad [nro].q_edge [0]);
711 Quad* qd = repQuadV (nh, pat->pat_quad [nro].q_edge [2]);
712 Quad* qe = repQuadV (nh, pat->pat_quad [nro].q_edge [1]);
713 Quad* qf = repQuadV (nh, pat->pat_quad [nro].q_edge [3]);
714 repHexa (nh, nro, qa, qb, qc, qd, qe, qf);