2 // C++ : Table d'hexaedres
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 "HexDocument.hxx"
26 #include "HexVector.hxx"
27 #include "HexVertex.hxx"
28 #include "HexHexa.hxx"
29 #include "HexEdge.hxx"
31 #include "HexGlobale.hxx"
32 #include "HexCylinder.hxx"
33 #include "HexShape.hxx"
39 void geom_dump_asso (Edge* edge);
40 void geom_create_circle (double* milieu, double rayon, double* normale,
41 double* base, string& brep);
43 // ====================================================== getHexaIJK
44 Hexa* Elements::getHexaIJK (int nx, int ny, int nz)
46 if (nx<0 || nx>=size_hx || ny<0 || ny>=size_hy || nz<0 || nz>=size_hz)
51 int nro = nx + size_hx*ny + size_hx*size_hy*nz;
53 return tab_hexa [nro];
55 // ====================================================== getQuadIJ
56 Quad* Elements::getQuadIJ (int nx, int ny, int nz)
58 if (nx<0 || nx>=size_qx || ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
63 int nro = nx + size_qx*ny + size_qx*size_qy*nz
64 + size_qx*size_qy*size_qz*dir_z;
65 return tab_quad [nro];
67 // ====================================================== getQuadJK
68 Quad* Elements::getQuadJK (int nx, int ny, int nz)
70 if (nx<0 || nx>=size_qx || ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
75 int nro = nx + size_qx*ny + size_qx*size_qy*nz; // + dir_x*...
77 return tab_quad [nro];
79 // ====================================================== getQuadIK
80 Quad* Elements::getQuadIK (int nx, int ny, int nz)
82 if (nx<0 || nx>=size_qx || ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
87 int nro = nx + size_qx*ny + size_qx*size_qy*nz + size_qx*size_qy*size_qz;
89 return tab_quad [nro];
91 // ====================================================== getEdgeI
92 Edge* Elements::getEdgeI (int nx, int ny, int nz)
94 if (nx<0 || nx>=size_ex || ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
99 int nro = nx + size_ex*ny + size_ex*size_ey*nz;
101 return tab_edge [nro];
103 // ====================================================== getEdgeJ
104 Edge* Elements::getEdgeJ (int nx, int ny, int nz)
106 if (nx<0 || nx>=size_ex || ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
108 else if (grid_nocart)
111 int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez;
113 return tab_edge [nro];
115 // ====================================================== getEdgeK
116 Edge* Elements::getEdgeK (int nx, int ny, int nz)
118 if (nx<0 || nx>=size_ex || ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
120 else if (grid_nocart)
123 int nro = nx + size_ex*ny + size_ex*size_ey*nz
124 + size_ex*size_ey*size_ez*dir_z;
125 return tab_edge [nro];
127 // ====================================================== getVertexIJK
128 Vertex* Elements::getVertexIJK (int nx, int ny, int nz)
130 if (nx<0 || nx>=size_vx || ny<0 || ny>=size_vy || nz<0 || nz>=size_vz)
132 else if (grid_nocart)
135 int nro = nx + size_vx*ny + size_vx*size_vy*nz;
137 return tab_vertex [nro];
139 // ====================================================== setVertex
140 void Elements::setVertex (Vertex* elt, int nx, int ny, int nz)
142 if ( nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy
143 || nz < 0 || nz >= size_vz) return;
145 int nro = nx + size_vx*ny + size_vx*size_vy*nz;
146 tab_vertex [nro] = elt;
148 // ====================================================== setVertex (2)
149 void Elements::setVertex (int nx, int ny, int nz, double px, double py,
152 if ( nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy
153 || nz < 0 || nz >= size_vz) return;
155 Vertex* node = el_root->addVertex (px, py, pz);
156 setVertex (node, nx, ny, nz);
158 // ====================================================== setEdge
159 void Elements::setEdge (Edge* elt, EnumCoord dir, int nx, int ny, int nz)
161 if (nx<0 || nx>=size_ex || ny<0 || ny>=size_ey || nz<0 || nz>=size_ez
162 || dir < dir_x || dir > dir_z )
165 int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez*dir;
166 tab_edge [nro] = elt;
168 // ====================================================== setQuad
169 void Elements::setQuad (Quad* elt, EnumCoord dir, int nx, int ny, int nz)
171 if (nx<0 || nx>=size_ex || ny<0 || ny>=size_ey || nz<0 || nz>=size_ez
172 || dir < dir_x || dir > dir_z )
175 int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez*dir;
176 tab_quad [nro] = elt;
178 // ====================================================== setHexa
179 void Elements::setHexa (Hexa* elt, int nx, int ny, int nz)
181 if ( nx < 0 || nx >= size_hx || ny < 0 || ny >= size_hy
182 || nz < 0 || nz >= size_hz) return;
184 int nro = nx + size_hx*ny + size_hx*size_hy*nz;
185 tab_hexa [nro] = elt;
187 // ====================================================== remove
188 void Elements::remove ()
190 int nbre=tab_hexa.size ();
192 for (int nh=0 ; nh<nbre ; nh++)
193 if (tab_hexa[nh] != NULL)
194 tab_hexa[nh]->remove();
196 // ====================================================== makeCylinder
197 int Elements::makeCylinder (Cylinder* cyl, Vector* vx, int nr, int na, int nl)
199 Vertex* orig = cyl->getBase ();
200 Vector* dir = cyl->getDirection ();
201 double ray = cyl->getRadius ();
202 double haut = cyl->getHeight ();
204 resize (GR_CYLINDRIC, nr, na, nl);
206 makeCylindricalNodes (orig, vx, dir, ray/(nr+1), 360, haut/nl,
209 assoCylinder (orig, dir, 360);
212 // ====================================================== makePipe
213 int Elements::makePipe (Cylinder* cyl, Vector* vx, int nr, int na, int nl)
215 Vertex* orig = cyl->getBase ();
216 Vector* dir = cyl->getDirection ();
217 double ray = cyl->getRadius ();
218 double haut = cyl->getHeight ();
220 resize (GR_CYLINDRIC, nr, na, nl);
222 makeCylindricalNodes (orig, vx, dir, ray, 360, haut, nr, na, nl, false);
224 assoCylinder (orig, dir, 360);
228 // ---------------------------------------- prism Quads
230 // ====================================================== prismQuads
231 int Elements::prismQuads (Quads& tstart, Vector* dir, int nbiter)
233 el_root->markAll (NO_USED);
234 int nbcells = tstart.size ();
238 nbr_hexas = nbcells*nbiter;
240 tab_hexa.resize (nbr_hexas);
241 tab_quad.clear (); // verticaux
242 ker_hquad.clear (); // Horizontaux
249 gen_matrix.defTranslation (dir);
251 for (int nro=0 ; nro<nbcells ; nro++)
253 prismHexas (nro, tstart[nro], nbiter);
259 // ====================================================== prismQuadsVec
260 int Elements::prismQuadsVec (Quads& tstart, Vector* dir, RealVector& tlen,
263 int nbiter = tlen.size();
267 el_root->markAll (NO_USED);
268 int nbcells = tstart.size ();
272 nbr_hexas = nbcells*nbiter;
274 tab_hexa.resize (nbr_hexas);
275 tab_quad.clear (); // verticaux
276 ker_hquad.clear (); // Horizontaux
283 dir->getCoord (prism_dir);
284 normer_vecteur (prism_dir);
287 for (int nro=0 ; nro<nbcells ; nro++)
289 prismHexas (nro, tstart[nro], nbiter);
295 // ======================================================== revolutionQuads
296 int Elements::revolutionQuads (Quads& start, Vertex* center, Vector* axis,
299 int nbiter = angles.size();
300 if (center==NULL || axis==NULL || nbiter==0)
303 el_root->markAll (NO_USED);
304 int nbcells = start.size ();
308 nbr_hexas = nbcells*nbiter;
310 tab_hexa.resize (nbr_hexas);
311 tab_quad.clear (); // verticaux
312 ker_hquad.clear (); // Horizontaux
320 revo_center = center;
323 for (int nro=0 ; nro<nbcells ; nro++)
325 prismHexas (nro, start[nro], nbiter);
331 // ====================================================== prismHexas
332 int Elements::prismHexas (int nro, Quad* qbase, int hauteur)
334 int ind_node [QUAD4];
337 // ----------------------------- Vertex + aretes verticales
338 for (int ns=0 ; ns<QUAD4 ; ns++)
340 Vertex* vbase = qbase ->getVertex (ns);
341 int indx = vbase->getMark ();
345 vbase->setMark (indx);
351 Real3 centre, vk, point, om;
352 revo_center->getPoint (centre);
353 vbase ->getPoint (point);
354 revo_axis ->getCoord (vk);
357 calc_vecteur (centre, point, om);
358 double oh = prod_scalaire (om, vk);
361 for (int dd=dir_x; dd<=dir_z ; dd++)
363 ph [dd] = centre [dd] + oh*vk[dd];
364 hm [dd] = point [dd] - ph[dd];
365 rayon += hm[dd] * hm[dd];
367 rayon = sqrt (rayon);
368 /********************************
376 ********************************/
377 geom_create_circle (ph, rayon, vk, hm, c_rep);
380 for (int nh=0 ; nh<hauteur ; nh++)
382 nd1 = el_root->addVertex (nd0->getX(), nd0->getY(), nd0->getZ());
384 gen_matrix.perform (nd1);
385 tab_vertex.push_back (nd1);
386 Edge* pilier = newEdge (nd0, nd1);
387 tab_pilier.push_back (pilier);
391 beta = alpha + gen_values[nh];
392 Shape* shape = new Shape (c_rep);
393 shape->setBounds (alpha/360, beta/360);
394 pilier->addAssociation (shape);
395 // geom_dump_asso (pilier);
400 ind_node [ns] = indx;
402 // ----------------------------- Aretes horizontales
403 // ----------------------------- + face verticales
404 int ind_poutre [QUAD4];
405 for (int ns=0 ; ns<QUAD4 ; ns++)
407 Edge* ebase = qbase->getEdge (ns);
408 int indx = ebase->getMark ();
412 ebase->setMark (indx);
413 int nd1 = ind_node [ns];
414 int nd2 = ind_node [(ns+1) MODULO QUAD4];
416 Edge *ed1, *ed2, *ed3;
417 for (int nh=0 ; nh<hauteur ; nh++)
420 ed0 = newEdge (tab_vertex [nd1*hauteur + nh],
421 tab_vertex [nd2*hauteur + nh]);
422 ed1 = tab_pilier [nd1*hauteur + nh];
423 ed3 = tab_pilier [nd2*hauteur + nh];
425 Quad* mur = newQuad (ed0, ed1, ed2, ed3);
426 tab_edge.push_back (ed0);
427 tab_quad.push_back (mur);
428 prismAssociation (ed2, ed0, nh, ed1);
431 ind_poutre [ns] = indx;
433 // ----------------------------- Faces horizontales
434 // ----------------------------- + Hexas
436 Quad *qb, *qc, *qd, *qe, *qf;
437 int nv0 = hauteur*ind_poutre [0];
438 int nv1 = hauteur*ind_poutre [1];
439 int nv2 = hauteur*ind_poutre [2];
440 int nv3 = hauteur*ind_poutre [3];
441 for (int nh=0 ; nh<hauteur ; nh++)
443 qb = newQuad (tab_edge [nh+nv0], tab_edge [nh+nv1],
444 tab_edge [nh+nv2], tab_edge [nh+nv3]);
445 qc = tab_quad [nh + nv0];
446 qd = tab_quad [nh + nv2];
447 qe = tab_quad [nh + nv1];
448 qf = tab_quad [nh + nv3];
450 // *** tab_hexa [nh*hauteur + nro] = newHexa (qa, qb, qc, qd, qe, qf); Abu
451 tab_hexa [nro*hauteur + nh] = newHexa (qa, qb, qc, qd, qe, qf);
452 ker_hquad.push_back (qb);
457 // ====================================================== updateMatrix
458 void Elements::updateMatrix (int hauteur)
462 gen_matrix.defRotation (revo_center, revo_axis, gen_values[hauteur]);
466 double h0 = hauteur>0 ? gen_values[hauteur-1] : 0;
467 double dh = gen_values[hauteur] - h0;
469 for (int nc=dir_x ; nc<=dir_z ; nc++)
470 decal [nc] = prism_dir [nc]*dh;
471 gen_matrix.defTranslation (decal);
474 // ====================================================== endPrism
475 void Elements::endPrism ()
477 int nbelts = ker_hquad.size();
478 for (int nro=0 ; nro<nbelts ; nro++)
479 tab_quad.push_back (ker_hquad[nro]);
481 nbelts = tab_pilier.size();
482 for (int nro=0 ; nro<nbelts ; nro++)
483 tab_edge.push_back (tab_pilier[nro]);
486 nbr_hexas = tab_hexa.size ();
487 nbr_edges = tab_edge.size ();
488 nbr_quads = tab_quad.size ();
489 nbr_vertex = tab_vertex.size ();