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"
24 #include "HexDocument.hxx"
25 #include "HexVector.hxx"
26 #include "HexVertex.hxx"
27 #include "HexHexa.hxx"
28 #include "HexEdge.hxx"
29 #include "HexGlobale.hxx"
38 // ====================================================== Constructeur
39 Elements::Elements (Document* doc) : EltBase (doc, EL_GRID)
41 glob = Globale::getInstance ();
44 size_qx = size_ex = size_vx = size_hx = 0;
45 size_qy = size_ey = size_vy = size_hy = 0;
46 size_qz = size_ez = size_vz = size_hz = 0;
47 size_qvplus = size_qhplus = size_ehplus = size_evplus = size_hplus = 0;
49 nbr_hexas = nbr_quads = nbr_edges = nbr_vertex = 0;
54 cyl_dispo = CYL_NOFILL;
58 // ====================================================== Constructeur
59 Elements::Elements (Document* doc, int nx, int ny, int nz)
60 : EltBase (doc, EL_GRID)
62 glob = Globale::getInstance ();
65 size_qx = size_ex = size_vx = size_hx = 0;
66 size_qy = size_ey = size_vy = size_hy = 0;
67 size_qz = size_ez = size_vz = size_hz = 0;
68 size_qvplus = size_qhplus = size_ehplus = size_evplus = size_hplus = 0;
70 nbr_hexas = nbr_quads = nbr_edges = nbr_vertex = 0;
73 cyl_dispo = CYL_NOFILL;
75 resize (GR_CYLINDRIC, nx, ny, nz);
78 // ====================================================== Constructeur (clonage)
79 Elements::Elements (Elements* orig) : EltBase (orig->el_root)
81 glob = Globale::getInstance ();
83 grid_type = orig->grid_type;
84 cyl_closed = orig->cyl_closed;
85 cyl_fill = orig->cyl_fill;
86 cyl_dispo = orig->cyl_dispo;
88 size_qx = size_ex = size_vx = size_hx = 0;
89 size_qy = size_ey = size_vy = size_hy = 0;
90 size_qz = size_ez = size_vz = size_hz = 0;
91 size_qvplus = size_qhplus = size_ehplus = size_evplus = size_hplus = 0;
92 nbr_hexas = nbr_quads = nbr_edges = nbr_vertex = 0;
94 resize (orig->grid_type, orig->size_hx, orig->size_hy, orig->size_hz);
95 cyl_closed = orig->cyl_closed;
97 // ====================================================== resize
98 void Elements::resize (EnumGrid type, int nx, int ny, int nz, int nplus)
109 size_hx = std::max (nx, 1);
110 size_hy = std::max (ny, 1);
111 size_hz = std::max (nz, 1);
113 size_qx = size_ex = size_vx = size_hx + 1;
114 size_qy = size_ey = size_vy = size_hy + 1;
115 size_qz = size_ez = size_vz = size_hz + 1;
117 nbr_hexas = size_hx * size_hy * size_hz;
118 nbr_quads = size_qx * size_qy * size_qz * DIM3;
119 nbr_edges = size_ex * size_ey * size_ez * DIM3;
120 nbr_vertex = size_vx * size_vy * size_vz;
124 size_qx = size_ex = size_vx = size_hx = 0;
125 size_qy = size_ey = size_vy = size_hy = 0;
126 size_qz = size_ez = size_vz = size_hz = 0;
127 nbr_quads = nbr_edges = nbr_vertex = 0;
128 gr_rayon = std::max (nx, 1);
130 nbr_hexas = 1 + gr_rayon*HQ_MAXI;
132 ker_vertex = nbr_vertex;
136 nbr_orig = std::max (nx, 1);
137 gr_hauteur = std::max (ny, 1);
140 size_hz = gr_hauteur;
142 nbr_hexas = nbr_orig * gr_hauteur;
143 nbr_vertex = nbr_hexas * QUAD4;
144 nbr_vertex = nbr_orig * (gr_hauteur+1)*QUAD4;
145 nbr_quads = nbr_vertex;
146 nbr_edges = 2*nbr_vertex;
150 nbr_orig = std::max (nx, 1); // nb quads du pattern
151 gr_hauteur = ny + 1; // Hauteur des hexas
152 pat_nbvertex = std::max (nz, 1); // nb vertex du pattern
153 pat_nbedges = std::max (nplus, 1); // nb edges du pattern
156 size_hz = gr_hauteur;
158 nbr_hexas = nbr_orig * gr_hauteur;
159 nbr_vertex = pat_nbvertex * (gr_hauteur+1);
160 nbr_edges = pat_nbedges * (gr_hauteur+1) + pat_nbvertex*gr_hauteur;
161 nbr_quads = nbr_orig * (gr_hauteur+1) + pat_nbedges *gr_hauteur;
170 size_qx = size_ex = size_vx = size_hx + 1;
171 size_qy = size_ey = size_vy = size_hy + 1;
172 size_qz = size_ez = size_vz = size_hz + 1;
174 nbr_hexas = size_hx * size_hy * size_hz;
175 nbr_quads = size_qx * size_qy * size_qz * DIM3;
176 nbr_edges = size_ex * size_ey * size_ez * DIM3;
177 nbr_vertex = size_vx * size_vy * size_vz;
181 tab_hexa .resize (nbr_hexas);
182 tab_quad .resize (nbr_quads);
183 tab_edge .resize (nbr_edges);
184 tab_vertex.resize (nbr_vertex);
186 ker_vertex = nbr_vertex;
188 for (int nc=0 ; nc< nbr_hexas ; nc++) tab_hexa [nc] = NULL;
189 for (int nc=0 ; nc< nbr_quads ; nc++) tab_quad [nc] = NULL;
190 for (int nc=0 ; nc< nbr_edges ; nc++) tab_edge [nc] = NULL;
191 for (int nc=0 ; nc< nbr_vertex ; nc++) tab_vertex [nc] = NULL;
194 // ====================================================== makeCartesianGrid
195 int Elements::makeCartesianGrid (Vertex* orig, Vector* v1, Vector* v2,
196 Vector* v3, int px, int py, int pz, int mx, int my, int mz)
198 if (BadElement (orig) || BadElement(v1) || BadElement(v2) || BadElement(v3)
199 || px<=0 || py<=0 || pz <= 0
200 || v1->getNorm () <= Epsil || v2->getNorm () <= Epsil
201 || v3->getNorm () <= Epsil)
207 resize (GR_CARTESIAN, px+mx, py+my, pz+mz);
209 makeCartesianNodes (orig, v1, v2, v3, px, py, pz, mx, my, mz);
214 // ====================================================== makeCylindricalGrid
215 int Elements::makeCylindricalGrid (Vertex* c, Vector* b, Vector* h,
216 double dr, double da, double dl, int nr, int na, int nl, bool fill)
218 if (BadElement (c) || BadElement(b) || BadElement(h)
219 || nr<=0 || na<=0 || nl <= 0 || dr < Epsil || da < Epsil || dl < Epsil
220 || b->getNorm () <= Epsil || h->getNorm () <= Epsil)
225 resize (GR_CYLINDRIC, nr, na, nl);
226 cyl_closed = da >= 360.0;
227 makeCylindricalNodes (c, b, h, dr, da, dl, nr, na, nl, fill);
229 assoCylinder (c, h, da);
232 // ====================================================== makeSphericalGrid
233 int Elements::makeSphericalGrid (Vertex* c, Vector* dv, int nb, double k)
235 if (BadElement (c) || BadElement(dv)
236 || nb<=0 || k < Epsil
237 || dv->getDx()<=Epsil || dv->getDy()<=Epsil || dv->getDz()<=Epsil)
243 resize (GR_SPHERIC, nb);
245 Vertex* i_node [HV_MAXI]; // Les noeuds de l'hexa englobant
246 Edge* i_edge [HE_MAXI]; // Les noeuds de l'hexa englobant
247 Quad* i_quad [HQ_MAXI]; // Les noeuds de l'hexa englobant
249 for (int nro=0 ; nro<HV_MAXI; nro++)
251 double dx = glob->CoordVertex (nro, dir_x) * dv->getDx();
252 double dy = glob->CoordVertex (nro, dir_y) * dv->getDy();
253 double dz = glob->CoordVertex (nro, dir_z) * dv->getDz();
255 i_node [nro] = el_root->addVertex (c->getX ()+dx, c->getY ()+dy,
259 for (int nro=0 ; nro<HE_MAXI; nro++)
261 int v1 = glob->EdgeVertex (nro, V_AMONT);
262 int v2 = glob->EdgeVertex (nro, V_AVAL);
263 i_edge[nro] = newEdge (i_node[v1], i_node[v2]);
267 char nm0[8], nm1 [8], nm2 [8];
268 printf (" %2d : %s = %s = [%s, %s] = [%d,%d] = [%s,%s]\n", nro,
269 glob->namofHexaEdge(nro), i_edge[nro]->getName(nm0),
270 glob->namofHexaVertex(v1), glob->namofHexaVertex(v2), v1, v2,
271 i_node[v1]->getName(nm1), i_node[v2]->getName(nm2));
275 for (int nro=0 ; nro<HQ_MAXI; nro++)
276 i_quad[nro] = newQuad (i_edge[glob->QuadEdge (nro, E_A)],
277 i_edge[glob->QuadEdge (nro, E_B)],
278 i_edge[glob->QuadEdge (nro, E_C)],
279 i_edge[glob->QuadEdge (nro, E_D)]);
281 tab_hexa.push_back (newHexa (i_quad[Q_A], i_quad[Q_B], i_quad[Q_C],
282 i_quad[Q_D], i_quad[Q_E], i_quad[Q_F]));
285 for (int niv=0; niv<gr_rayon ; niv++)
287 double lambda0 = lambda;
290 addStrate (i_quad, i_edge, i_node, c, lambda/lambda0);
295 // ====================================================== addStrate
296 int Elements::addStrate (Quad* i_quad[], Edge* i_edge[], Vertex* i_node[],
297 Vertex* center, double lambda)
299 Vertex* e_node [HV_MAXI]; // Les noeuds de l'hexa englobant
300 Edge* e_edge [HE_MAXI]; // Les noeuds de l'hexa englobant
301 Quad* e_quad [HQ_MAXI]; // Les noeuds de l'hexa englobant
303 Edge* d_edge [HV_MAXI]; // Les aretes diagonales (1 par sommet)
304 Quad* d_quad [HE_MAXI]; // Les faces diagonales (1 par arete)
307 // + les aretes diagonales
308 for (int nv=0 ; nv<HV_MAXI ; nv++)
310 double px0 = center->getX ();
311 double py0 = center->getY ();
312 double pz0 = center->getZ ();
313 e_node[nv] = el_root->addVertex (px0+lambda*(i_node[nv]->getX()-px0),
314 py0+lambda*(i_node[nv]->getY()-py0),
315 pz0+lambda*(i_node[nv]->getZ()-pz0));
317 d_edge[nv] = newEdge (i_node[nv], e_node[nv]);
319 // Les aretes exterieures
320 // + les faces diagonales
321 for (int nro=0 ; nro<HE_MAXI ; nro++)
323 int nv0 = glob->EdgeVertex (nro, V_AMONT);
324 int nv1 = glob->EdgeVertex (nro, V_AVAL );
325 e_edge[nro] = newEdge (e_node [nv0], e_node [nv1]);
326 d_quad[nro] = newQuad (i_edge [nro], d_edge [nv0],
327 e_edge [nro], d_edge [nv1]);
329 // Les faces exterieures
332 for (int nro=0 ; nro<HQ_MAXI ; nro++)
334 int ne0 = glob->QuadEdge (nro, E_A);
335 int ne1 = glob->QuadEdge (nro, E_B);
336 int ne2 = glob->QuadEdge (nro, E_C);
337 int ne3 = glob->QuadEdge (nro, E_D);
339 e_quad[nro] = newQuad (e_edge[ne0], e_edge[ne1],
340 e_edge[ne2], e_edge[ne3]);
341 strate = newHexa (i_quad[nro], e_quad[nro], d_quad[ne0],
342 d_quad[ne2], d_quad[ne1], d_quad[ne3]);
343 tab_hexa.push_back (strate);
346 for (int nv=0 ; nv<HV_MAXI ; nv++) i_node [nv] = e_node [nv];
347 for (int ns=0 ; ns<HE_MAXI ; ns++) i_edge [ns] = e_edge [ns];
348 for (int nq=0 ; nq<HQ_MAXI ; nq++) i_quad [nq] = e_quad [nq];
352 // ====================================================== saveVtk
353 int Elements::saveVtk (cpchar nomfic)
355 // -- 1) Raz numerotation precedente
356 for (int nro=0 ; nro<nbr_hexas ; nro++)
358 Hexa* cell = tab_hexa[nro];
359 if (cell!=NULL && cell->isHere())
366 for (int nro=0 ; nro<nbr_hexas ; nro++)
368 Hexa* cell = tab_hexa[nro];
369 if (cell!=NULL && cell->isHere())
372 nbnodes += cell->countNodes ();
376 pfile vtk = fopen (nomfic, "w");
377 fprintf (vtk, "# vtk DataFile Version 3.1\n");
378 fprintf (vtk, "%s \n", nomfic);
379 fprintf (vtk, "ASCII\n");
380 fprintf (vtk, "DATASET UNSTRUCTURED_GRID\n");
381 fprintf (vtk, "POINTS %d float\n", nbnodes);
385 for (int nro=0 ; nro<nbr_hexas ; nro++)
387 Hexa* cell = tab_hexa[nro];
388 if (cell!=NULL && cell->isHere())
389 cell->printNodes (vtk, nronode);
393 fprintf (vtk, "CELLS %d %d\n", nbcells, nbcells*(HV_MAXI+1));
395 for (int nro=0 ; nro<nbr_hexas ; nro++)
397 Hexa* cell = tab_hexa[nro];
398 if (cell!=NULL && cell->isHere())
399 cell->printHexa (vtk);
402 fprintf (vtk, "CELL_TYPES %d\n", nbcells);
403 for (int nro=0 ; nro<nbcells ; nro++)
404 fprintf (vtk, "%d\n", HE_MAXI);
409 // ============ = = = = = = = = = Vertices ..........
410 // ====================================================== newEdge
411 Edge* Elements::newEdge (Vertex* v1, Vertex* v2)
413 if (v1==NULL || v2==NULL)
416 Edge* elt = new Edge (v1, v2);
419 // ====================================================== newQuad
420 Quad* Elements::newQuad (Edge* e1, Edge* e2, Edge* e3, Edge* e4)
422 if (e1==NULL || e2==NULL || e3==NULL|| e4==NULL)
425 Quad* elt = new Quad (e1, e2, e3, e4);
428 // ====================================================== newHexa
429 Hexa* Elements::newHexa (Quad* qa, Quad* qb, Quad* qc, Quad* qd,
432 if (qa==NULL || qb==NULL || qc==NULL|| qd==NULL || qe==NULL|| qf==NULL)
435 Hexa* elt = new Hexa (qa, qb, qc, qd, qe, qf);
438 // ====================================================== joinQuads
439 int Elements::joinQuads (Quads& orig, int nb, Vertex* v1, Vertex* v2,
440 Vertex* v3, Vertex* v4, Quad* cible)
442 resize (GR_JOINT, orig.size(), nb);
444 el_root->markAll (IS_NONE);
446 // db = el_root->debug ();
449 nbr_orig = orig.size();
451 for (int nro=0 ; nro<nbr_orig ; nro++)
453 Quad* face = orig[nro];
457 printf (" *** joinQuads : donnees incorrectes\n");
458 printf (" *** le %deme quadrangle de depart est NULL\n", nro);
462 else if (face->getNbrParents()>1)
465 printf (" *** joinQuads : donnees incorrectes\n");
466 printf (" *** le %deme quadrangle de depart n'est pas une "
467 "face externe\n", nro);
472 orig [nro]->setMark (nro);
473 tab_orig.push_back (orig[nro]);
476 Edge* e_orig = tab_orig[0] -> findEdge (v1, v3);
477 Edge* e_dest = cible -> findEdge (v2, v4);
482 printf (" *** joinQuads : donnees incorrectes\n");
483 printf (" *** Les vertex v1 et v3 passes en argument ne definissent\n");
484 printf (" *** pas une arete du 1er quadrangle de depart\n");
494 printf (" *** joinQuads : donnees incorrectes\n");
495 printf (" *** Les vertex v2 et v4 passes en argument ne definissent\n");
496 printf (" *** pas une arete du quadrangle cible\n");
503 if (e_orig==NULL || e_dest==NULL)
509 StrOrient orient (v1, v3, v2, v4);
510 int ier = this ->coupler (0, cible, &orient);
516 ier = orig[0]->coupler (cible, &orient, this);
519 // ======================================================== coupler
520 int Elements::coupler (int nquad, Quad* dest, StrOrient* orient)
522 Quad* orig = tab_orig [nquad];
524 setQuad (orig, dir_z, 0, 0, 0);
525 setQuad (dest, dir_z, nquad, 0, 0);
527 int n11 = orig->indexVertex (orient->v11);
528 int n12 = orig->indexVertex (orient->v12);
529 int n21 = dest->indexVertex (orient->v21);
530 int n22 = dest->indexVertex (orient->v22);
532 // ---------------- Les 4 sommets initiaux
533 Vertex* vorig[QUAD4] = { orient->v11, orient->v12,
534 orig->getVertex((n11+2) MODULO QUAD4),
535 orig->getVertex((n12+2) MODULO QUAD4) };
537 Vertex* vdest[QUAD4] = { orient->v21, orient->v22,
538 dest->getVertex((n21+2) MODULO QUAD4),
539 dest->getVertex((n22+2) MODULO QUAD4) };
542 printf ("Quad nro %d : ", nquad);
543 orig->printName (" est couple avec ");
544 dest->printName ("\n");
545 printf ("Orientation : (");
546 for (int ii=0 ; ii<QUAD4 ; ii++) printf("%s ", vorig[ii]->getName());
549 for (int ii=0 ; ii<QUAD4 ; ii++) printf("%s ", vdest[ii]->getName());
553 // ---------------- Les sommets + les aretes verticales
554 for (int ns=0 ; ns< QUAD4 ; ns++)
556 Vertex* nd1 = vorig[ns];
557 Vertex* nd2 = vdest[ns];
558 int nref0 = nd1->getMark();
559 tab_vertex [nroVertex (nquad,ns,0)] = nd1;
563 double px0 = nd1->getX();
564 double py0 = nd1->getY();
565 double pz0 = nd1->getZ();
567 double dx = (nd2->getX() - nd1->getX()) / gr_hauteur;
568 double dy = (nd2->getY() - nd1->getY()) / gr_hauteur;
569 double dz = (nd2->getZ() - nd1->getZ()) / gr_hauteur;
571 nd1->setMark (indVertex (nquad, ns, 0));
574 for (int nh=0 ; nh<gr_hauteur ; nh++)
578 if (nh == gr_hauteur-1)
581 nd = el_root->addVertex (px0 + nh1*dx, py0 + nh1*dy,
583 int nv = indVertex (nquad, ns, nh);
584 tab_vertex [nv] = nd;
585 tab_edge [nv] = el_root->addEdge (ndp, nd);
587 printf (" Edge vertical nro %d = %s = (%s, %s)\n", nv,
588 tab_edge[nv]->getName(),
589 ndp->getName(), nd->getName());
594 for (int nh=0 ; nh<gr_hauteur ; nh++)
596 int nv = indVertex (nquad, ns, nh);
597 int nv0 = indVertex (nref0, nh);
598 tab_vertex [nv] = tab_vertex [nv0];
599 tab_edge [nv] = tab_edge [nv0];
603 // ---------------- Les quads verticaux + aretes horiz
604 for (int ns=0 ; ns< QUAD4 ; ns++)
606 int next = (ns+1) MODULO QUAD4;
607 Edge* arete = orig->findEdge (vorig[ns], vorig[next]);
608 int nref0 = arete->getMark();
609 int nref = indVertex (nquad, ns, 0);
610 // Construction des faces & aretes H
613 arete->setMark (nref);
618 for (int nh=0 ; nh< gr_hauteur ; nh++)
620 nva = indVertex (nquad, ns, nh);
621 nvb = indVertex (nquad, next, nh);
622 nha = nroEdgeH (nva);
627 if (nh==gr_hauteur-1)
628 ea = dest->findEdge (vdest[ns], vdest[next]);
630 ea = el_root->addEdge (tab_vertex [nva], tab_vertex [nvb]);
632 propagateAssociation (ec, ea, eb);
633 tab_quad [nva] = newQuad (ea, eb, ec, ed);
634 if (BadElement (tab_quad [nva]))
639 // L'arete a deja ete traitee
642 for (int nh=0 ; nh<gr_hauteur ; nh++)
644 int nva = indVertex (nquad, ns, nh);
645 int nha = nroEdgeH (nva);
647 int nvb = indVertex (nref0, nh);
648 int nhb = nroEdgeH (nvb);
650 tab_quad [nva] = tab_quad [nvb];
651 tab_edge [nha] = tab_edge [nhb];
655 // -------------------------------- Les Hexas
657 Quad *fa, *fc, *fd, *fe, *ff;
658 for (int nh=0 ; nh< gr_hauteur ; nh++)
661 if (nh == gr_hauteur-1)
664 fb = newQuad (tab_edge [nroEdgeH (nquad, E_A, nh)],
665 tab_edge [nroEdgeH (nquad, E_B, nh)],
666 tab_edge [nroEdgeH (nquad, E_C, nh)],
667 tab_edge [nroEdgeH (nquad, E_D, nh)]);
669 fc = tab_quad [indVertex (nquad, E_A, nh)];
670 fd = tab_quad [indVertex (nquad, E_C, nh)];
671 fe = tab_quad [indVertex (nquad, E_B, nh)];
672 ff = tab_quad [indVertex (nquad, E_D, nh)];
674 tab_hexa [nroHexa(nquad,nh)] = newHexa (fa,fb,fc,fd,fe,ff);
675 if (BadElement (tab_hexa [nroHexa(nquad,nh)]))
680 // ====================================================== makeCartesianNodes
681 int Elements::makeCartesianNodes (Vertex* orig, Vector* v1, Vector* v2,
682 Vector* v3, int px, int py, int pz, int mx, int my, int mz)
684 double dx = v1->getDx() + v2->getDx() + v3->getDx();
685 double dy = v1->getDy() + v2->getDy() + v3->getDy();
686 double dz = v1->getDz() + v2->getDz() + v3->getDz();
688 double px0 = orig->getX () - mx * dx;
689 double py0 = orig->getY () - my * dy;
690 double pz0 = orig->getZ () - mz * dz;
694 for (int nz=0 ; nz<size_vz ; nz++)
695 for (int ny=0 ; ny<size_vy ; ny++)
696 for (int nx=0 ; nx<size_vx ; nx++)
699 if (nx!=mx || ny!=my || nz!=mz)
700 node = el_root->addVertex (px0 + nx*dx, py0 + ny*dy,
702 setVertex (node, nx, ny, nz);
707 // ====================================================== makeCylindricalNodes
708 int Elements::makeCylindricalNodes (Vertex* orig, Vector* base, Vector* haut,
709 double dr, double da, double dl, int nr, int na, int nl, bool fill)
711 int ier = makeBasicCylinder (dr, da, dl, nr, na, nl, fill);
715 transfoVertices (orig, base, haut);
718 // ====================================================== transfoVertices
719 void Elements::transfoVertices (Vertex* orig, Vector* base, Vector* haut)
721 Vector* iprim = new Vector (base);
722 Vector* jprim = new Vector (base);
723 Vector* kprim = new Vector (haut);
725 int ier = kprim->renormer ();
729 jprim->vectoriel (kprim, base);
730 ier = jprim->renormer ();
734 iprim->vectoriel (jprim, kprim);
735 transfoVertices (orig, iprim, jprim, kprim);
737 // ====================================================== transfoVertices
738 void Elements::transfoVertices (Vertex* orig, Vector* iprim, Vector* jprim,
741 double matrice[DIM3][DIM3]={{iprim->getDx(),jprim->getDx(),kprim->getDx()},
742 {iprim->getDy(),jprim->getDy(),kprim->getDy()},
743 {iprim->getDz(),jprim->getDz(),kprim->getDz()}};
745 double matkx = orig->getX();
746 double matky = orig->getY();
747 double matkz = orig->getZ();
749 int nbre = tab_vertex.size ();
750 for (int nro=0 ; nro<nbre ; nro++)
752 if (tab_vertex[nro] != NULL)
753 tab_vertex[nro]->setMark (NO_USED);
756 for (int nro=0 ; nro<nbre ; nro++)
758 Vertex* node = tab_vertex[nro];
759 if (node != NULL && node->getMark() == NO_USED)
761 double point [DIM3] = {node->getX(), node->getY(), node->getZ()};
762 double result[DIM3] = {matkx, matky, matkz};
764 for (int ni=0 ; ni<DIM3; ni++)
765 for (int nj=0 ; nj<DIM3; nj++)
766 result [ni] += matrice[ni][nj] * point[nj];
768 node->setCoord (result[dir_x], result[dir_y], result[dir_z]);
769 node->setMark (IS_USED);
773 // ====================================================== transform
774 int Elements::transform (Matrix* matrice)
777 for (int nro=0 ; nro<nbr_hexas ; nro++)
779 Hexa* cell = tab_hexa[nro];
780 if (cell!=NULL && cell->isHere())
784 for (int nro=0 ; nro<nbr_hexas ; nro++)
786 Hexa* cell = tab_hexa[nro];
787 if (cell!=NULL && cell->isHere())
788 cell->moveNodes (matrice);
793 // -- Cas pathologique : il n'y avait pas d'hexas
794 // -- On se rabat sur les sommets
796 for (int nro=0 ; nro<nbr_vertex ; nro++)
798 Vertex* node = tab_vertex[nro];
799 if (node!=NULL && node->isHere())
800 matrice->perform (node);
805 // ====================================================== cutHexas
806 int Elements::cutHexas (const Edges& t_edges, int nbcuts)
808 // 1) marquage des hexas
809 el_root->markAll (NO_USED);
811 vector <Quad*> q_amont;
812 vector <Quad*> q_aval;
813 map <Vertex*, Vertex*> vis_a_vis;
815 int nbnodes = t_edges.size();
816 vector <Vertex*> v_amont (nbnodes);
817 vector <Vertex*> v_aval (nbnodes);
820 for (int nro=0; nro<nbnodes ; nro++)
822 Edge* arete = t_edges [nro];
823 v_amont [nro] = arete->getAmont ();
824 v_aval [nro] = arete->getAval ();
827 printf (" %3d : Edge = (", nro);
828 v_amont[nro]->printName (", ");
829 v_aval [nro]->printName (")\n");
832 vis_a_vis [v_amont[nro]] = v_aval[nro];
833 vis_a_vis [v_aval[nro]] = v_amont[nro];
834 int nbcells = arete->getNbrParents ();
836 for (int nq=0 ; nq<nbcells ; nq++)
838 Quad* quad = arete->getParent (nq);
839 if (quad->getMark () != IS_USED)
841 quad->setMark (IS_USED);
842 int nbcubes = quad->getNbrParents ();
843 for (int nh=0 ; nh<nbcubes ; nh++)
845 Hexa* hexa = quad->getParent (nh);
846 if (hexa->getMark () != IS_USED)
848 hexa->setMark (IS_USED);
849 int namont = hexa->getBase (v_amont[nro], arete);
850 int naval = glob->getOpposedQuad (namont);
851 q_amont.push_back (hexa->getQuad (namont));
852 q_aval .push_back (hexa->getQuad (naval ));
856 printf (" %3d : Quad = ", nbfaces);
857 hexa->printName (", ");
858 printf (" Faces = (");
859 hexa->getQuad (namont)->printName (", ");
860 hexa->getQuad (naval )->printName (")\n");
868 // ------------------- Dimensionnement
869 int nbcells = q_amont.size ();
870 nbr_vertex = nbnodes*(nbcuts+2);
871 int nbpiliers = nbnodes*(nbcuts+1); // aretes verticales
872 int nbpoutres = nbcells*(nbcuts+2)*QUAD4; // aretes horizontales
873 nbr_edges = nbpoutres;
874 nbr_quads = nbcells*(nbcuts+1)*QUAD4; // faces Verticales
875 nbr_hexas = nbcells*(nbcuts+1);
877 // ------------------- Les noeuds et les aretes verticales
878 tab_quad.resize (nbr_quads);
879 tab_edge.resize (nbr_edges);
880 tab_hexa.resize (nbr_hexas);
881 tab_vertex.resize (nbr_vertex);
882 vector <Edge*> tab_pilier (nbpiliers);
884 int nbinter = nbcuts + 1;
885 for (int ned=0; ned<nbnodes ; ned++)
887 Shapes ass_shapes = t_edges[ned] -> getAssociations ();
889 t_edges [ned]->remove ();
890 Vertex* ndamont = v_amont [ned];
891 Vertex* ndaval = v_aval [ned];
893 double dx = (ndaval->getX() - ndamont->getX()) / nbinter;
894 double dy = (ndaval->getY() - ndamont->getY()) / nbinter;
895 double dz = (ndaval->getZ() - ndamont->getZ()) / nbinter;
897 Vertex* nd0 = tab_vertex [ned] = ndamont;
898 for (int nc=0; nc<nbcuts ; nc++)
901 Vertex* nd1 = el_root->addVertex (ndamont->getX() + nc1*dx,
902 ndamont->getY() + nc1*dy,
903 ndamont->getZ() + nc1*dz);
904 tab_vertex [nc1*nbnodes + ned] = nd1;
905 tab_pilier [nc *nbnodes + ned] = newEdge (nd0, nd1);
906 ass_edges.push_back (tab_pilier [nc *nbnodes + ned]);
909 tab_vertex [nbinter*nbnodes + ned] = ndaval;
910 tab_pilier [nbcuts *nbnodes + ned] = newEdge (nd0, ndaval);
911 ass_edges.push_back (tab_pilier[nbcuts *nbnodes + ned]);
912 ndamont->setMark (ned);
913 cutAssociation (ass_shapes, ass_edges);
915 // ------------------- Les aretes horizontales
916 // ------------------- Les faces verticales
917 HexDisplay (nbcells);
918 int sizelig = nbcells*QUAD4;
919 for (int nro=0; nro<nbcells ; nro++)
921 Quad* sol = q_amont [nro];
922 Quad* toit = q_aval [nro];
923 for (int ns=0; ns<QUAD4 ; ns++)
925 Edge* plinthe = sol->getEdge (ns);
926 int nmur = nro*QUAD4 + ns;
927 int nmur0 = plinthe->getMark();
930 for (int nc=0 ; nc<nbinter ; nc++)
932 tab_edge [nc*sizelig + nmur] = tab_edge [nc*sizelig + nmur0];
933 tab_quad [nc*sizelig + nmur] = tab_quad [nc *sizelig + nmur0];
936 printf (" %2d : %d quad_vertical [%02d] =", nro, ns,
938 printf (" quad_vertical [%02d]\n", nc*sizelig + nmur0);
941 tab_edge [nbinter*sizelig+nmur] = tab_edge[nbinter*sizelig+nmur0];
945 plinthe->setMark (nmur);
946 Vertex* vs1 = sol->getVertex (ns);
947 Vertex* vs2 = sol->getVertex ((ns+1) MODULO QUAD4);
948 int nd1 = vs1->getMark ();
949 int nd2 = vs2->getMark ();
950 Edge* ed0 = tab_edge [nmur] = plinthe;
954 for (int nc=0 ; nc<nbinter ; nc++)
959 v1 = tab_vertex [nc1*nbnodes + nd1];
960 v2 = tab_vertex [nc1*nbnodes + nd2];
961 ed2 = newEdge (v1, v2);
965 v1 = vis_a_vis [vs1];
966 v2 = vis_a_vis [vs2];
967 ed2 = toit->findEdge (v1, v2);
970 tab_edge [nc1*sizelig + nmur] = ed2;
971 tab_quad [nc *sizelig + nmur] = newQuad (ed0,
972 tab_pilier [nc*nbnodes + nd1], ed2,
973 tab_pilier [nc*nbnodes + nd2]);
977 printf (" %2d : %d quad_vertical [%02d] = ", nro, ns,
979 PrintName (tab_quad [nc *sizelig + nmur]);
986 // ------------------- Les faces horizontales
987 // ------------------- Les hexas
988 // Rappel : sizelig = nbcells*QUAD4
989 for (int nro=0; nro<nbcells ; nro++)
991 Quad* qa = q_amont [nro];
992 for (int nc=0 ; nc<nbinter ; nc++)
994 int quadv = nc*nbcells*QUAD4 + nro*QUAD4;
995 Quad* qb = q_aval[nro];
998 int edh = (nc+1)*nbcells*QUAD4 + nro*QUAD4;
999 qb = newQuad (tab_edge[edh], tab_edge[edh+1],
1000 tab_edge[edh+2], tab_edge[edh+3]);
1004 Quad* qc = tab_quad [quadv];
1005 Quad* qd = tab_quad [quadv + 2];
1006 Quad* qe = tab_quad [quadv + 1];
1007 Quad* qf = tab_quad [quadv + 3];
1008 Hexa* hexa = newHexa (qa, qb, qc, qd, qe, qf);
1009 if (BadElement(hexa))
1011 tab_hexa [nc*nbcells + nro] = hexa;
1017 // ====================================================== clear_associations
1018 void clear_associations (std::vector<Quad*>& table)
1020 int nbelts = table.size();
1021 for (int nro=0 ; nro<nbelts ; nro++)
1023 Quad* elt = table [nro];
1024 if (elt != NULL && elt->isValid())
1025 elt->clearAssociation ();
1028 // ====================================================== clear_associations
1029 void clear_associations (std::vector<Edge*>& table)
1031 int nbelts = table.size();
1032 for (int nro=0 ; nro<nbelts ; nro++)
1034 Edge* elt = table [nro];
1035 if (elt != NULL && elt->isValid())
1036 elt->clearAssociation ();
1039 // ====================================================== clear_associations
1040 void clear_associations (std::vector<Vertex*>& table)
1042 int nbelts = table.size();
1043 for (int nro=0 ; nro<nbelts ; nro++)
1045 Vertex* elt = table [nro];
1046 if (elt != NULL && elt->isValid())
1047 elt->clearAssociation ();
1050 // ====================================================== clearAssociation
1051 void Elements::clearAssociation ()
1053 clear_associations (tab_quad);
1054 clear_associations (tab_edge);
1055 clear_associations (tab_vertex);
1058 clear_associations (ker_hquad);
1059 clear_associations (ker_vquad);
1060 clear_associations (ker_hedge);
1061 clear_associations (ker_vedge);
1063 // ============================================================ findVertex
1064 int Elements::findVertex (double vx, double vy, double vz)
1066 double tol = el_root->getTolerance ();
1067 double xmin = vx - tol;
1068 double xmax = vx + tol;
1069 double ymin = vy - tol;
1070 double ymax = vy + tol;
1071 double zmin = vz - tol;
1072 double zmax = vz + tol;
1074 int nbre = tab_vertex.size();
1075 for (int nro=0 ; nro<nbre ; nro++)
1077 Vertex* node = tab_vertex [nro];
1078 if (node != NULL && node->isHere ()
1079 && node->isin (xmin, xmax, ymin, ymax, zmin, zmax))