2 // C++ : Table d'hexaedres
4 #include "HexElements.hxx"
6 #include "HexDocument.hxx"
7 #include "HexVector.hxx"
8 #include "HexVertex.hxx"
10 #include "HexEdge.hxx"
12 #include "HexGlobale.hxx"
13 #include "HexCylinder.hxx"
21 // ====================================================== getHexaIJK
22 Hexa* Elements::getHexaIJK (int nx, int ny, int nz)
24 if (nx<0 || nx>=size_hx || ny<0 || ny>=size_hy || nz<0 || nz>=size_hz)
26 else if (grid_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
29 int nro = nx + size_hx*ny + size_hx*size_hy*nz;
31 return tab_hexa [nro];
33 // ====================================================== getQuadIJ
34 Quad* Elements::getQuadIJ (int nx, int ny, int nz)
36 if (nx<0 || nx>=size_qx || ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
38 else if (grid_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
41 int nro = nx + size_qx*ny + size_qx*size_qy*nz
42 + size_qx*size_qy*size_qz*dir_z;
43 return tab_quad [nro];
45 // ====================================================== getQuadJK
46 Quad* Elements::getQuadJK (int nx, int ny, int nz)
48 if (nx<0 || nx>=size_qx || ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
50 else if (grid_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
53 int nro = nx + size_qx*ny + size_qx*size_qy*nz; // + dir_x*...
55 return tab_quad [nro];
57 // ====================================================== getQuadIK
58 Quad* Elements::getQuadIK (int nx, int ny, int nz)
60 if (nx<0 || nx>=size_qx || ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
62 else if (grid_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
65 int nro = nx + size_qx*ny + size_qx*size_qy*nz + size_qx*size_qy*size_qz;
67 return tab_quad [nro];
69 // ====================================================== getEdgeI
70 Edge* Elements::getEdgeI (int nx, int ny, int nz)
72 if (nx<0 || nx>=size_ex || ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
74 else if (grid_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
77 int nro = nx + size_ex*ny + size_ex*size_ey*nz;
79 return tab_edge [nro];
81 // ====================================================== getEdgeJ
82 Edge* Elements::getEdgeJ (int nx, int ny, int nz)
84 if (nx<0 || nx>=size_ex || ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
86 else if (grid_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
89 int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez;
91 return tab_edge [nro];
93 // ====================================================== getEdgeK
94 Edge* Elements::getEdgeK (int nx, int ny, int nz)
96 if (nx<0 || nx>=size_ex || ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
98 else if (grid_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
101 int nro = nx + size_ex*ny + size_ex*size_ey*nz
102 + size_ex*size_ey*size_ez*dir_z;
103 return tab_edge [nro];
105 // ====================================================== getVertexIJK
106 Vertex* Elements::getVertexIJK (int nx, int ny, int nz)
108 if (nx<0 || nx>=size_vx || ny<0 || ny>=size_vy || nz<0 || nz>=size_vz)
110 else if (grid_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
113 int nro = nx + size_vx*ny + size_vx*size_vy*nz;
115 return tab_vertex [nro];
117 // ====================================================== setVertex
118 void Elements::setVertex (Vertex* elt, int nx, int ny, int nz)
120 if ( nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy
121 || nz < 0 || nz >= size_vz) return;
123 int nro = nx + size_vx*ny + size_vx*size_vy*nz;
124 tab_vertex [nro] = elt;
126 // ====================================================== setVertex (2)
127 void Elements::setVertex (int nx, int ny, int nz, double px, double py,
130 if ( nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy
131 || nz < 0 || nz >= size_vz) return;
133 Vertex* node = el_root->addVertex (px, py, pz);
134 setVertex (node, nx, ny, nz);
136 // ====================================================== setEdge
137 void Elements::setEdge (Edge* elt, EnumCoord dir, int nx, int ny, int nz)
139 if (nx<0 || nx>=size_ex || ny<0 || ny>=size_ey || nz<0 || nz>=size_ez
140 || dir < dir_x || dir > dir_z )
143 int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez*dir;
144 tab_edge [nro] = elt;
146 // ====================================================== setQuad
147 void Elements::setQuad (Quad* elt, EnumCoord dir, int nx, int ny, int nz)
149 if (nx<0 || nx>=size_ex || ny<0 || ny>=size_ey || nz<0 || nz>=size_ez
150 || dir < dir_x || dir > dir_z )
153 int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez*dir;
154 tab_quad [nro] = elt;
156 // ====================================================== setHexa
157 void Elements::setHexa (Hexa* elt, int nx, int ny, int nz)
159 if ( nx < 0 || nx >= size_hx || ny < 0 || ny >= size_hy
160 || nz < 0 || nz >= size_hz) return;
162 int nro = nx + size_hx*ny + size_hx*size_hy*nz;
163 tab_hexa [nro] = elt;
165 // ====================================================== remove
166 void Elements::remove ()
168 int nbre=tab_hexa.size ();
170 for (int nh=0 ; nh<nbre ; nh++)
171 if (tab_hexa[nh] != NULL)
172 tab_hexa[nh]->remove();
174 // ====================================================== prismQuads
175 int Elements::prismQuads (Quads& tstart, Vector* dir, int hauteur)
177 el_root->markAll (NO_USED);
178 int nbcells = tstart.size ();
182 nbr_hexas = nbcells*hauteur;
184 tab_hexa.resize (nbr_hexas);
185 tab_quad.clear (); // verticaux
190 for (int nro=0 ; nro<nbcells ; nro++)
192 pushHexas (nro, tstart[nro], dir->getDx(), dir->getDy(), dir->getDz(),
197 // ====================================================== pushHexas
198 int Elements::pushHexas (int nro, Quad* qbase, double dx, double dy,
199 double dz, int hauteur, Quad* toit, int tcorr[])
201 int ind_node [QUAD4];
203 // ----------------------------- Vertex + aretes verticales
204 for (int ns=0 ; ns<QUAD4 ; ns++)
206 Vertex* vbase = qbase ->getVertex (ns);
207 int indx = vbase->getMark ();
211 vbase->setMark (indx);
214 for (int nh=0 ; nh<hauteur ; nh++)
216 nd1 = el_root->addVertex (nd0->getX() + dx, nd0->getY() + dy,
218 tab_vertex.push_back (nd1);
219 tab_pilier.push_back (newEdge (nd0, nd1));
223 ind_node [ns] = indx;
225 // ----------------------------- Aretes horizontales
226 // ----------------------------- + face verticales
227 int ind_poutre [QUAD4];
228 for (int ns=0 ; ns<QUAD4 ; ns++)
230 Edge* ebase = qbase ->getEdge (ns);
231 int indx = ebase->getMark ();
235 ebase->setMark (indx);
236 int nd1 = ind_node [ns];
237 int nd2 = ind_node [(ns+1) MODULO QUAD4];
239 Edge *ed1, *ed2, *ed3;
240 for (int nh=0 ; nh<hauteur ; nh++)
243 ed0 = newEdge (tab_vertex [nd1*hauteur + nh],
244 tab_vertex [nd2*hauteur + nh]);
245 ed1 = tab_pilier [nd1*hauteur + nh];
246 ed3 = tab_pilier [nd2*hauteur + nh];
248 Quad* mur = newQuad (ed0, ed1, ed2, ed3);
249 tab_edge.push_back (ed0);
250 tab_quad.push_back (mur);
253 ind_poutre [ns] = indx;
255 // ----------------------------- Faces horizontales
256 // ----------------------------- + Hexas
258 Quad *qb, *qc, *qd, *qe, *qf;
259 int nv0 = hauteur*ind_poutre [0];
260 int nv1 = hauteur*ind_poutre [1];
261 int nv2 = hauteur*ind_poutre [2];
262 int nv3 = hauteur*ind_poutre [3];
263 for (int nh=0 ; nh<hauteur ; nh++)
265 qb = newQuad (tab_edge [nh+nv0], tab_edge [nh+nv1],
266 tab_edge [nh+nv2], tab_edge [nh+nv3]);
267 qc = tab_quad [nh + nv0];
268 qd = tab_quad [nh + nv2];
269 qe = tab_quad [nh + nv1];
270 qf = tab_quad [nh + nv3];
272 // *** tab_hexa [nh*hauteur + nro] = newHexa (qa, qb, qc, qd, qe, qf); Abu
273 tab_hexa [nro*hauteur + nh] = newHexa (qa, qb, qc, qd, qe, qf);
278 // ====================================================== cutHexas
279 int Elements::cutHexas0 (const Edges& t_edges, int nbcuts)
281 // 1) marquage des hexas
282 el_root->markAll (NO_USED);
284 vector <Quad*> q_amont;
285 vector <Quad*> q_aval;
287 int nbnodes = t_edges.size();
288 vector <Vertex*> v_amont (nbnodes);
289 vector <Vertex*> v_aval (nbnodes);
291 std::map <Vertex*, Vertex*> vertex_aval;
293 for (int nro=0; nro<nbnodes ; nro++)
295 Edge* arete = t_edges [nro];
296 Vertex* amont = arete->getAmont ();
297 v_amont [nro] = amont;
298 vertex_aval [amont] = arete->getAval ();
299 int nbcells = arete->getNbrParents ();
301 for (int nq=0 ; nq<nbcells ; nq++)
303 Quad* quad = arete->getParent (nq);
304 if (quad->getMark () != IS_USED)
306 quad->setMark (IS_USED);
307 int nbcubes = quad->getNbrParents ();
308 for (int nh=0 ; nh<nbcubes ; nh++)
310 Hexa* hexa = quad->getParent (nh);
311 if (hexa->getMark () != IS_USED)
313 hexa->setMark (IS_USED);
314 int namont = hexa->getBase (v_amont[nro], arete);
315 int naval = glob->getOpposedQuad (namont);
316 q_amont.push_back (hexa->getQuad (namont));
317 q_aval .push_back (hexa->getQuad (naval ));
322 hexa->printName (", ");
323 printf (" Faces = (");
324 hexa->getQuad (namont)->printName (", ");
325 hexa->getQuad (naval )->printName (")\n");
332 // ------------------- Dimensionnement
333 int hauteur = nbcuts + 1;
334 int nbcells = q_amont.size ();
335 nbr_hexas = nbcells * hauteur;
337 tab_hexa.resize (nbr_hexas);
338 tab_quad.clear (); // verticaux
343 for (int ned=0; ned<nbnodes ; ned++)
344 t_edges [ned]->remove ();
347 for (int nro=0; nro<nbcells ; nro++)
350 Quad* sol = q_amont[nro];
351 Quad* toit = q_aval [nro];
353 for (int nv=0; nv<QUAD4 ; nv++)
355 Vertex* dest = vertex_aval [sol->getVertex(nv)];
356 tcorr [nv] = toit->indexVertex (dest);
359 pushHexas (nro, sol, 0,0,0, hauteur, toit, tcorr);
363 // ====================================================== makeCylinder
364 int Elements::makeCylinder (Cylinder* cyl, Vector* vx, int nr, int na, int nl)
366 Vertex* orig = cyl->getBase ();
367 Vector* dir = cyl->getDirection ();
368 double ray = cyl->getRadius ();
369 double haut = cyl->getHeight ();
371 resize (GR_CYLINDRIC, nr, na, nl);
373 makeCylindricalNodes (orig, vx, dir, ray, 360, haut, nr, na, nl, true);
377 // ====================================================== makePipe
378 int Elements::makePipe (Cylinder* cyl, Vector* vx, int nr, int na, int nl)
380 Vertex* orig = cyl->getBase ();
381 Vector* dir = cyl->getDirection ();
382 double ray = cyl->getRadius ();
383 double haut = cyl->getHeight ();
385 resize (GR_CYLINDRIC, nr, na, nl);
387 makeCylindricalNodes (orig, vx, dir, ray, 360, haut, nr, na, nl, false);