2 // C++ : Table d'hexaedres
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 "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 "HexOldShape.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)
48 if (nx<0 || nx>=size_hx || ny<0 || ny>=size_hy || nz<0 || nz>=size_hz)
53 int nro = nx + size_hx*ny + size_hx*size_hy*nz;
55 DumpStart ("getHexaIJK", nx << ny << nz);
56 DumpReturn (tab_hexa [nro]);
57 return tab_hexa [nro];
59 // ====================================================== getQuadIJ
60 Quad* Elements::getQuadIJ (int nx, int ny, int nz)
64 if (nx<0 || nx>=size_qx || ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
69 int nro = nx + size_qx*ny + size_qx*size_qy*nz
70 + size_qx*size_qy*size_qz*dir_z;
72 DumpStart ("getQuadIJ", nx << ny << nz);
73 DumpReturn (tab_quad [nro]);
74 return tab_quad [nro];
76 // ====================================================== getQuadJK
77 Quad* Elements::getQuadJK (int nx, int ny, int nz)
81 if (nx<0 || nx>=size_qx || ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
86 int nro = nx + size_qx*ny + size_qx*size_qy*nz; // + dir_x*...
88 DumpStart ("getQuadJK", nx << ny << nz);
89 DumpReturn (tab_quad [nro]);
90 return tab_quad [nro];
92 // ====================================================== getQuadIK
93 Quad* Elements::getQuadIK (int nx, int ny, int nz)
97 if (nx<0 || nx>=size_qx || ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
102 int nro = nx + size_qx*ny + size_qx*size_qy*nz + size_qx*size_qy*size_qz;
104 DumpStart ("getQuadIK", nx << ny << nz);
105 DumpReturn (tab_quad [nro]);
106 return tab_quad [nro];
108 // ====================================================== getEdgeI
109 Edge* Elements::getEdgeI (int nx, int ny, int nz)
113 if (nx<0 || nx>=size_ex || ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
115 else if (grid_nocart)
118 int nro = nx + size_ex*ny + size_ex*size_ey*nz;
120 DumpStart ("getEdgeI", nx << ny << nz);
121 DumpReturn (tab_edge [nro]);
122 return tab_edge [nro];
124 // ====================================================== getEdgeJ
125 Edge* Elements::getEdgeJ (int nx, int ny, int nz)
129 if (nx<0 || nx>=size_ex || ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
131 else if (grid_nocart)
134 int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez;
136 DumpStart ("getEdgeJ", nx << ny << nz);
137 DumpReturn (tab_edge [nro]);
138 return tab_edge [nro];
140 // ====================================================== getEdgeK
141 Edge* Elements::getEdgeK (int nx, int ny, int nz)
145 if (nx<0 || nx>=size_ex || ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
147 else if (grid_nocart)
150 int nro = nx + size_ex*ny + size_ex*size_ey*nz
151 + size_ex*size_ey*size_ez*dir_z;
153 DumpStart ("getEdgeK", nx << ny << nz);
154 DumpReturn (tab_edge [nro]);
155 return tab_edge [nro];
157 // ====================================================== getVertexIJK
158 Vertex* Elements::getVertexIJK (int nx, int ny, int nz)
162 if (nx<0 || nx>=size_vx || ny<0 || ny>=size_vy || nz<0 || nz>=size_vz)
164 else if (grid_nocart)
167 int nro = nx + size_vx*ny + size_vx*size_vy*nz;
169 DumpStart ("getVertexIJK", nx << ny << nz);
170 DumpReturn (tab_vertex [nro]);
171 return tab_vertex [nro];
173 // ====================================================== setVertex
174 void Elements::setVertex (Vertex* elt, int nx, int ny, int nz)
176 if ( nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy
177 || nz < 0 || nz >= size_vz) return;
179 int nro = nx + size_vx*ny + size_vx*size_vy*nz;
180 tab_vertex [nro] = elt;
182 // ====================================================== setVertex (2)
183 void Elements::setVertex (int nx, int ny, int nz, double px, double py,
186 if ( nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy
187 || nz < 0 || nz >= size_vz) return;
189 Vertex* node = el_root->addVertex (px, py, pz);
190 setVertex (node, nx, ny, nz);
192 // ====================================================== setEdge
193 void Elements::setEdge (Edge* elt, EnumCoord dir, int nx, int ny, int nz)
195 if (nx<0 || nx>=size_ex || ny<0 || ny>=size_ey || nz<0 || nz>=size_ez
196 || dir < dir_x || dir > dir_z )
199 int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez*dir;
200 tab_edge [nro] = elt;
202 // ====================================================== setQuad
203 void Elements::setQuad (Quad* elt, EnumCoord dir, int nx, int ny, int nz)
205 if (nx<0 || nx>=size_ex || ny<0 || ny>=size_ey || nz<0 || nz>=size_ez
206 || dir < dir_x || dir > dir_z )
209 int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez*dir;
210 tab_quad [nro] = elt;
212 // ====================================================== setHexa
213 void Elements::setHexa (Hexa* elt, int nx, int ny, int nz)
215 if ( nx < 0 || nx >= size_hx || ny < 0 || ny >= size_hy
216 || nz < 0 || nz >= size_hz) return;
218 int nro = nx + size_hx*ny + size_hx*size_hy*nz;
219 tab_hexa [nro] = elt;
221 // ====================================================== remove
222 void Elements::remove ()
224 int nbre=tab_hexa.size ();
226 for (int nh=0 ; nh<nbre ; nh++)
227 if (tab_hexa[nh] != NULL)
228 tab_hexa[nh]->remove();
230 // ====================================================== makeCylinder
231 int Elements::makeCylinder (Cylinder* cyl, Vector* vx, int nr, int na, int nl)
233 if (BadElement (cyl) || BadElement (vx) || nr<=0 || na <=3 || nl <=0
234 || vx->getNorm () <= Epsil)
240 Vertex* orig = cyl->getBase ();
241 Vector* dir = cyl->getDirection ();
242 double ray = cyl->getRadius ();
243 double haut = cyl->getHeight ();
245 resize (GR_CYLINDRIC, nr, na, nl);
247 makeCylindricalNodes (orig, vx, dir, ray/(nr+1), 360, haut/nl,
250 assoCylinder (orig, dir, 360);
253 // ====================================================== makePipe
254 int Elements::makePipe (Cylinder* cyl, Vector* vx, int nr, int na, int nl)
256 if (BadElement (cyl) || BadElement (vx) || nr<=0 || na <=3 || nl <=0
257 || vx->getNorm () <= Epsil)
263 Vertex* orig = cyl->getBase ();
264 Vector* dir = cyl->getDirection ();
265 double ray = cyl->getRadius ();
266 double haut = cyl->getHeight ();
268 resize (GR_CYLINDRIC, nr, na, nl);
270 makeCylindricalNodes (orig, vx, dir, ray, 360, haut, nr, na, nl, false);
272 assoCylinder (orig, dir, 360);
276 // ---------------------------------------- prism Quads
278 // ====================================================== prismQuads
279 int Elements::prismQuads (Quads& tstart, Vector* dir, int nbiter)
281 if (BadElement (dir) || dir->getNorm () <= Epsil || nbiter <= 0)
287 el_root->markAll (NO_USED);
288 int nbcells = tstart.size ();
292 nbr_hexas = nbcells*nbiter;
294 tab_hexa.resize (nbr_hexas);
295 tab_quad.clear (); // verticaux
296 ker_hquad.clear (); // Horizontaux
303 gen_matrix.defTranslation (dir);
305 for (int nro=0 ; nro<nbcells ; nro++)
307 prismHexas (nro, tstart[nro], nbiter);
313 // ====================================================== prismQuadsVec
314 int Elements::prismQuadsVec (Quads& tstart, Vector* dir, RealVector& tlen,
317 int nbiter = tlen.size();
318 if (BadElement (dir) || dir->getNorm () <= Epsil || nbiter <= 0)
324 el_root->markAll (NO_USED);
325 int nbcells = tstart.size ();
329 nbr_hexas = nbcells*nbiter;
331 tab_hexa.resize (nbr_hexas);
332 tab_quad.clear (); // verticaux
333 ker_hquad.clear (); // Horizontaux
340 dir->getCoord (prism_dir);
341 normer_vecteur (prism_dir);
344 for (int nro=0 ; nro<nbcells ; nro++)
346 prismHexas (nro, tstart[nro], nbiter);
352 // ======================================================== revolutionQuads
353 int Elements::revolutionQuads (Quads& start, Vertex* center, Vector* axis,
356 int nbiter = angles.size();
357 int nbcells = start.size ();
358 if (BadElement (center) || BadElement(axis) || nbiter==0 || nbcells==0
359 || axis->getNorm () <= Epsil)
365 el_root->markAll (NO_USED);
369 nbr_hexas = nbcells*nbiter;
371 tab_hexa.resize (nbr_hexas);
372 tab_quad.clear (); // verticaux
373 ker_hquad.clear (); // Horizontaux
381 revo_center = center;
384 for (int nro=0 ; nro<nbcells ; nro++)
386 prismHexas (nro, start[nro], nbiter);
392 // ====================================================== prismHexas
393 int Elements::prismHexas (int nro, Quad* qbase, int hauteur)
395 int ind_node [QUAD4];
398 // ----------------------------- Vertex + aretes verticales
399 for (int ns=0 ; ns<QUAD4 ; ns++)
401 Vertex* vbase = qbase ->getVertex (ns);
402 int indx = vbase->getMark ();
406 vbase->setMark (indx);
412 Real3 centre, vk, point, om;
413 revo_center->getPoint (centre);
414 vbase ->getPoint (point);
415 revo_axis ->getCoord (vk);
418 calc_vecteur (centre, point, om);
419 double oh = prod_scalaire (om, vk);
422 for (int dd=dir_x; dd<=dir_z ; dd++)
424 ph [dd] = centre [dd] + oh*vk[dd];
425 hm [dd] = point [dd] - ph[dd];
426 rayon += hm[dd] * hm[dd];
428 rayon = sqrt (rayon);
429 /********************************
437 ********************************/
438 geom_create_circle (ph, rayon, vk, hm, c_rep);
441 for (int nh=0 ; nh<hauteur ; nh++)
443 nd1 = el_root->addVertex (nd0->getX(), nd0->getY(), nd0->getZ());
445 gen_matrix.perform (nd1);
446 tab_vertex.push_back (nd1);
447 Edge* pilier = newEdge (nd0, nd1);
448 tab_pilier.push_back (pilier);
452 beta = alpha + gen_values[nh];
453 Shape* shape = new Shape (c_rep);
454 shape->setBounds (alpha/360, beta/360);
455 pilier->addAssociation (shape);
456 // geom_dump_asso (pilier);
461 ind_node [ns] = indx;
463 // ----------------------------- Aretes horizontales
464 // ----------------------------- + face verticales
465 int ind_poutre [QUAD4];
466 for (int ns=0 ; ns<QUAD4 ; ns++)
468 Edge* ebase = qbase->getEdge (ns);
469 int indx = ebase->getMark ();
473 ebase->setMark (indx);
474 int nd1 = ind_node [ns];
475 int nd2 = ind_node [(ns+1) MODULO QUAD4];
477 Edge *ed1, *ed2, *ed3;
478 for (int nh=0 ; nh<hauteur ; nh++)
481 ed0 = newEdge (tab_vertex [nd1*hauteur + nh],
482 tab_vertex [nd2*hauteur + nh]);
483 ed1 = tab_pilier [nd1*hauteur + nh];
484 ed3 = tab_pilier [nd2*hauteur + nh];
486 Quad* mur = newQuad (ed0, ed1, ed2, ed3);
487 tab_edge.push_back (ed0);
488 tab_quad.push_back (mur);
489 prismAssociation (ed2, ed0, nh, ed1);
492 ind_poutre [ns] = indx;
494 // ----------------------------- Faces horizontales
495 // ----------------------------- + Hexas
497 Quad *qb, *qc, *qd, *qe, *qf;
498 int nv0 = hauteur*ind_poutre [0];
499 int nv1 = hauteur*ind_poutre [1];
500 int nv2 = hauteur*ind_poutre [2];
501 int nv3 = hauteur*ind_poutre [3];
502 for (int nh=0 ; nh<hauteur ; nh++)
504 qb = newQuad (tab_edge [nh+nv0], tab_edge [nh+nv1],
505 tab_edge [nh+nv2], tab_edge [nh+nv3]);
506 qc = tab_quad [nh + nv0];
507 qd = tab_quad [nh + nv2];
508 qe = tab_quad [nh + nv1];
509 qf = tab_quad [nh + nv3];
511 // *** tab_hexa [nh*hauteur + nro] = newHexa (qa, qb, qc, qd, qe, qf); Abu
512 tab_hexa [nro*hauteur + nh] = newHexa (qa, qb, qc, qd, qe, qf);
513 ker_hquad.push_back (qb);
518 // ====================================================== updateMatrix
519 void Elements::updateMatrix (int hauteur)
523 gen_matrix.defRotation (revo_center, revo_axis, gen_values[hauteur]);
527 double h0 = hauteur>0 ? gen_values[hauteur-1] : 0;
528 double dh = gen_values[hauteur] - h0;
530 for (int nc=dir_x ; nc<=dir_z ; nc++)
531 decal [nc] = prism_dir [nc]*dh;
532 gen_matrix.defTranslation (decal);
535 // ====================================================== endPrism
536 void Elements::endPrism ()
538 int nbelts = ker_hquad.size();
539 for (int nro=0 ; nro<nbelts ; nro++)
540 tab_quad.push_back (ker_hquad[nro]);
542 nbelts = tab_pilier.size();
543 for (int nro=0 ; nro<nbelts ; nro++)
544 tab_edge.push_back (tab_pilier[nro]);
547 nbr_hexas = tab_hexa.size ();
548 nbr_edges = tab_edge.size ();
549 nbr_quads = tab_quad.size ();
550 nbr_vertex = tab_vertex.size ();