1 // C++ : Table d'hexaedres (Evol Versions 3)
3 // Copyright (C) 2009-2013 CEA/DEN, EDF R&D
5 // This library is free software; you can redistribute it and/or
6 // modify it under the terms of the GNU Lesser General Public
7 // License as published by the Free Software Foundation; either
8 // version 2.1 of the License.
10 // This library is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 // Lesser General Public License for more details.
15 // You should have received a copy of the GNU Lesser General Public
16 // License along with this library; if not, write to the Free Software
17 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 #include "HexElements.hxx"
24 #include "HexDocument.hxx"
25 #include "HexVector.hxx"
26 #include "HexVertex.hxx"
27 #include "HexEdge.hxx"
28 #include "HexHexa.hxx"
29 #include "HexMatrix.hxx"
30 #include "HexGlobale.hxx"
32 #include "HexNewShape.hxx"
33 #include "HexEdgeShape.hxx"
39 void translate_brep (string& brep, double vdir[], string& trep);
40 void transfo_brep (string& brep, Matrix* matrice, string& trep);
44 // ====================================================== makeRind
45 int Elements::makeRind (EnumGrid type, Vertex* center, Vector* vx, Vector* vz,
46 double radext, double radint, double radhole,
47 Vertex* plorig, double angle, int nr, int na, int nl)
50 int ier = controlRind (type, center, vx, vz, radext, radint, radhole,
51 plorig, angle, nr, na, nl, cyl_phi0, phi1);
55 resize (type, nr, na, nl);
59 cyl_radhole = radhole;
60 cyl_closed = type==GR_HEMISPHERIC || type==GR_RIND;
62 double theta = cyl_closed ? 2*M_PI : M_PI*angle/180;
63 cyl_dphi = (phi1-cyl_phi0)/ size_hz;
64 cyl_dtheta = theta / size_hy;
66 int nb_secteurs = cyl_closed ? size_vy-1 : size_vy;
68 for (int ny=0 ; ny<nb_secteurs ; ny++)
69 for (int nx=0 ; nx<size_vx ; nx++)
70 for (int nz=0 ; nz<size_vz ; nz++)
73 getCylPoint (nx, ny, nz, px, py, pz);
74 Vertex* node = el_root->addVertex (px, py, pz);
75 setVertex (node, nx, ny, nz);
78 for (int nx=0 ; nx<size_vx ; nx++)
79 for (int nz=0 ; nz<size_vz ; nz++)
81 Vertex* node = getVertexIJK (nx, 0, nz);
82 setVertex (node, nx, size_vy-1, nz);
85 transfoVertices (center, vx, vz);
87 assoCylinder (center, vz, angle);
90 // ====================================================== getCylPoint
91 int Elements::getCylPoint (int nr, int na, int nh, double& px, double& py,
94 if (grid_type == GR_CYLINDRIC)
96 px = cyl_radext * cos (na*cyl_dtheta);
97 py = cyl_radext * sin (na*cyl_dtheta);
102 bool rind = (grid_type == GR_RIND || grid_type == GR_PART_RIND);
103 double sinphi = sin (cyl_phi0 + nh*cyl_dphi);
104 double cosphi = cos (cyl_phi0 + nh*cyl_dphi);
109 rayon = cyl_radint + nr*(cyl_radext-cyl_radint)/size_hx;
111 rayon = rayon*cosphi;
115 pz = cyl_radext*sinphi;
116 rayon = cyl_radhole + nr*(cyl_radext*cosphi - cyl_radhole)/size_hx;
117 rayon = std::max (cyl_radhole, rayon);
120 px = rayon * cos (na*cyl_dtheta);
121 py = rayon * sin (na*cyl_dtheta);
125 // ====================================================== controlRind
126 int Elements::controlRind (EnumGrid type, Vertex* cx, Vector* vx, Vector* vz,
127 double rext, double rint, double rhole,
128 Vertex* px, double angle,
129 int nrad, int nang, int nhaut,
130 double &phi0, double &phi1)
132 const double Epsil1 = 1e-6;
135 if (cx == NULL || vx == NULL || vz == NULL)
138 if (nrad<=0 || nang<=0 || nhaut<=0)
150 double nvx = vx->getNorm();
151 double nvz = vz->getNorm();
153 if (nvx < Epsil1 || nvz < Epsil1)
156 double alpha = asin (rhole/rext);
157 double beta = -M_PI*DEMI;
158 if (type==GR_HEMISPHERIC || type==GR_PART_SPHERIC)
164 double oh = ((px->getX() - cx->getX()) * vz->getDx()
165 + (px->getY() - cx->getY()) * vz->getDy()
166 + (px->getZ() - cx->getZ()) * vz->getDz()) / nvz;
170 beta = asin (oh/rext);
173 phi0 = std::max (alpha - M_PI*DEMI, beta);
174 phi1 = M_PI*DEMI - alpha;
177 // ====================================================== getHexas
178 int Elements::getHexas (Hexas& liste)
181 for (int nro = 0 ; nro<nbr_hexas ; nro++)
183 Hexa* cell = tab_hexa [nro];
184 if (cell!=NULL && cell->isValid())
185 liste.push_back (cell);
189 // ====================================================== assoCylinder
190 void Elements::assoCylinder (Vertex* ori, Vector* normal, double angle)
193 assoCylinders (ori, normal, angle, tangles);
195 // ====================================================== assoCylinders
196 void Elements::assoCylinders (Vertex* ori, Vector* normal, double angle,
197 RealVector& t_angles)
200 sprintf (name, "grid_%02d", el_id);
201 NewShape* geom = el_root->addShape (name, SH_CYLINDER);
204 int na = t_angles.size();
205 bool regul = na == 0;
206 double alpha = angle/size_hy;
209 Real3 vk = { normal->getDx(), normal->getDy(), normal->getDz() };
212 for (int nz=0 ; nz<size_vz ; nz++)
214 for (int nx=0 ; nx<size_vx ; nx++)
216 Vertex* pm = getVertexIJK (nx, 0, nz);
217 Real3 om = { pm->getX() - ori->getX(),
218 pm->getY() - ori->getY(),
219 pm->getZ() - ori->getZ() };
221 double oh = prod_scalaire (om, vk);
224 for (int dd=dir_x; dd<=dir_z ; dd++)
226 ph [dd] = ori->getCoord(dd) + oh*vk[dd];
227 hm [dd] = pm ->getCoord(dd) - ph[dd];
228 rayon += hm[dd] * hm[dd];
231 rayon = sqrt (rayon);
232 int subid = geom->addCircle (ph, rayon, vk, hm);
235 for (int ny=0 ; ny<size_hy ; ny++)
237 double beta = regul ? alpha : t_angles [ny];
240 Edge* edge = getEdgeJ (nx, ny, nz);
241 geom->addAssociation (edge, subid, pmin, pmax);
242 // Shape* shape = new Shape (brep);
243 // shape-> setBounds (pmin, pmax);
244 // edge->addAssociation (shape);
248 // Association automatique des vertex non associes -> bph
249 // Traitement des faces spheriques
251 Real3 vi = { -vk[dir_x], -vk[dir_y], -vk[dir_z] };
252 Real3 po = { ori->getX(), ori->getY(), ori->getZ() };
256 case GR_HEMISPHERIC : // Pour l'exterieur
257 case GR_PART_SPHERIC :
258 assoRind (po, vi, size_vx-1, geom);
260 case GR_PART_RIND : // Exterieur + interieur
262 assoRind (po, vi, 0, geom);
263 assoRind (po, vi, size_vx-1, geom);
268 // assoResiduelle (); // Association des sommets residuels
271 // ====================================================== assoRind
272 // Association des meridiennes
273 // Creation sphere geometrique + association faces
274 void Elements::assoRind (double* ori, double* vi, int nx, NewShape* geom)
277 double radius = nx==0 ? cyl_radint : cyl_radext;
278 double paramin = std::max ((cyl_phi0 + M_PI/2) / (2*M_PI), 0.0);
279 int sphid = geom->addSphere (ori, radius);
282 int nb_secteurs = cyl_closed ? size_vy-1 : size_vy;
284 for (int ny=0 ; ny<nb_secteurs ; ny++)
286 Vertex* pm = getVertexIJK (nx, ny, nz1);
287 Real3 vj = { pm->getX(), pm->getY(), pm->getZ() };
288 prod_vectoriel (vi, vj, vk);
289 double rayon = cyl_radint + nx*(cyl_radext-cyl_radint)/size_hx;
290 int subid = geom->addCircle (ori, rayon, vk, vi);
291 double pmax = paramin;
293 for (int nz=0 ; nz<size_hz ; nz++)
295 Quad* quad = getQuadJK (nx, ny, nz);
296 Edge* edge = getEdgeK (nx, ny, nz);
298 pmax = pmin + cyl_dphi/(2*M_PI);
299 geom->addAssociation (edge, subid, pmin, pmax);
300 geom->addAssociation (quad, sphid);
304 // ====================================================== assoCircle
305 // ==== utilise pour les spheres carrees
306 void Elements::assoCircle (double* center, Edge* ed1, Edge* ed2, NewShape* geom)
308 Real3 oa, ob, normal;
309 Real3 pta, ptb, ptc, ptd;
312 // Les 2 edges dont les petits cotes d'un rectangle de rapport L/l=sqrt(2)
313 // Soit le cercle circonscrit a ce rectangle.
314 // * la largeur est balayee par l'angle alpha
315 // * la longueur par l'angle beta = pi -alpha
317 ed1->getVertex(V_AMONT)->getPoint (pta);
318 ed1->getVertex(V_AVAL )->getPoint (ptb);
319 ed2->getVertex(V_AMONT)->getPoint (ptc);
320 ed2->getVertex(V_AVAL )->getPoint (ptd);
322 double d1 = calc_distance (pta, ptc);
323 double d2 = calc_distance (pta, ptd);
327 ed2->getVertex(V_AMONT)->getPoint (ptd);
328 ed2->getVertex(V_AVAL )->getPoint (ptc);
331 calc_vecteur (center, pta, oa);
332 calc_vecteur (center, ptb, ob);
333 prod_vectoriel (oa, ob, normal);
335 double rayon = calc_norme (oa);
336 int subid = geom->addCircle (center, rayon, normal, oa);
338 // Shape* asso1 = new Shape (brep);
339 // Shape* asso2 = new Shape (brep);
341 const double alpha = atan (sqrt(2)/2)/M_PI; // angle proche de 70.528 degres
342 // asso1->setBounds (0, alpha);
343 // asso2->setBounds (0.5, alpha + 0.5);
345 geom->addAssociation (ed1, subid, 0.0, alpha);
346 geom->addAssociation (ed2, subid, 0.5, alpha+0.5);
348 // ====================================================== assoSphere
349 void Elements::assoSphere (Vertex* ori, Edge* t_edge[], Quad* t_quad[])
352 sprintf (name, "grid_%02d", el_id);
353 NewShape* geom = el_root->addShape (name, SH_SPHERE);
354 geom -> openShape ();
356 Real3 center, sommet;
357 ori->getPoint(center);
359 assoCircle (center, t_edge [E_AC], t_edge [E_BD], geom);
360 assoCircle (center, t_edge [E_AD], t_edge [E_BC], geom);
361 assoCircle (center, t_edge [E_AE], t_edge [E_BF], geom);
362 assoCircle (center, t_edge [E_AF], t_edge [E_BE], geom);
363 assoCircle (center, t_edge [E_CE], t_edge [E_DF], geom);
364 assoCircle (center, t_edge [E_CF], t_edge [E_DE], geom);
366 t_edge[E_AC]->getVertex(V_AMONT)->getPoint (sommet);
367 double radius = calc_distance (center, sommet);;
369 int sphid = geom -> addSphere (center, radius);
372 for (int nf=0 ; nf < HQ_MAXI ; nf++)
373 t_quad[nf]->addAssociation (geom, sphid);
375 // assoResiduelle (); // Association des sommets residuels
377 // ====================================================== makeSphericalGrid
378 int Elements::makeSphericalGrid (Vertex* c, double rayon, int nb, double k)
380 if (nb<=0 || rayon <=Epsil || k <= Epsil || BadElement (c))
386 resize (GR_SPHERIC, nb);
388 Vertex* i_node [HV_MAXI]; // Les noeuds de l'hexa englobant
389 Edge* i_edge [HE_MAXI]; // Les noeuds de l'hexa englobant
390 Quad* i_quad [HQ_MAXI]; // Les noeuds de l'hexa englobant
392 for (int nro=0 ; nro<HV_MAXI; nro++)
394 double dx = glob->CoordVertex (nro, dir_x) * rayon;
395 double dy = glob->CoordVertex (nro, dir_y) * rayon;
396 double dz = glob->CoordVertex (nro, dir_z) * rayon;
398 i_node [nro] = el_root->addVertex (c->getX ()+dx, c->getY ()+dy,
402 for (int nro=0 ; nro<HE_MAXI; nro++)
404 int v1 = glob->EdgeVertex (nro, V_AMONT);
405 int v2 = glob->EdgeVertex (nro, V_AVAL);
406 i_edge[nro] = newEdge (i_node[v1], i_node[v2]);
410 char nm0[8], nm1 [8], nm2 [8];
411 printf (" %2d : %s = %s = [%s, %s] = [%d,%d] = [%s,%s]\n", nro,
412 glob->namofHexaEdge(nro), i_edge[nro]->getName(nm0),
413 glob->namofHexaVertex(v1), glob->namofHexaVertex(v2), v1, v2,
414 i_node[v1]->getName(nm1), i_node[v2]->getName(nm2));
418 for (int nro=0 ; nro<HQ_MAXI; nro++)
419 i_quad[nro] = newQuad (i_edge[glob->QuadEdge (nro, E_A)],
420 i_edge[glob->QuadEdge (nro, E_B)],
421 i_edge[glob->QuadEdge (nro, E_C)],
422 i_edge[glob->QuadEdge (nro, E_D)]);
424 tab_hexa.push_back (newHexa (i_quad[Q_A], i_quad[Q_B], i_quad[Q_C],
425 i_quad[Q_D], i_quad[Q_E], i_quad[Q_F]));
428 for (int niv=0; niv<gr_rayon ; niv++)
430 double lambda0 = lambda;
433 addStrate (i_quad, i_edge, i_node, c, lambda/lambda0);
436 assoSphere (c, i_edge, i_quad);
439 // ==================================================== propagateAssociation
440 int Elements::propagateAssociation (Edge* orig, Edge* dest, Edge* dir)
444 if (revo_lution || orig==NULL || dest==NULL || dir==NULL)
447 const Shapes& tab_shapes = orig->getAssociations ();
448 const Shapes& tab_dest = dest->getAssociations ();
449 int nbdest = tab_dest.size();
450 int nbshapes = tab_shapes.size();
451 bool on_edge = nbshapes!=0 && nbdest==0;
453 Vertex* vo1 = orig->commonVertex (dir);
454 Vertex* vd1 = dest->commonVertex (dir);
455 Vertex* vo2 = orig->opposedVertex (vo1);
456 Vertex* vd2 = dest->opposedVertex (vd1);
458 if (vo1==NULL || vd1==NULL)
462 Real3 pa, pb, vdir1, vdir2;
463 calc_vecteur (vo1->getPoint (pa), vd1->getPoint (pb), vdir1);
464 calc_vecteur (vo2->getPoint (pa), vd2->getPoint (pb), vdir2);
466 double dd = calc_distance (vdir1, vdir2);
467 bool para = dd < 1.0e-3;
471 for (int nro=0 ; nro<nbshapes ; nro++)
473 Shape* shape = tab_shapes[nro];
476 string brep = shape->getBrep();
477 translate_brep (brep, vdir1, trep);
478 // Shape* tshape = new Shape (trep);
479 // tshape->setBounds (shape->getStart(), shape->getEnd());
480 // dest->addAssociation (tshape);
485 double* vdir = vdir1;
486 for (int nro=V_AMONT ; nro<=V_AVAL ; nro++)
488 Shape* shape = vo1->getAssociation ();
489 if (shape!=NULL && vd1->getAssociation ()==NULL)
491 string brep = shape->getBrep();
492 translate_brep (brep, vdir, trep);
493 // Shape* tshape = new Shape (trep);
494 // vd1->setAssociation (tshape);
504 // ==================================================== prismAssociation
505 int Elements::prismAssociation (Edge* orig, Edge* dest, int nh, Edge* dir)
509 if (orig==NULL || dest==NULL || dir==NULL)
514 const Shapes& tab_shapes = orig->getAssociations ();
515 const Shapes& tab_dest = dest->getAssociations ();
516 int nbdest = tab_dest.size();
517 int nbshapes = tab_shapes.size();
518 bool on_edge = nbshapes>0 && nbdest==0;
520 Vertex* vo1 = orig->commonVertex (dir);
521 Vertex* vd1 = dest->commonVertex (dir);
523 if (vo1==NULL || vd1==NULL)
527 Real3 porig, pdest, vdir;
528 vo1->getPoint (porig);
529 vd1->getPoint (pdest);
530 calc_vecteur (porig, pdest, vdir);
534 for (int nro=0 ; nro<nbshapes ; nro++)
536 Shape* shape = tab_shapes[nro];
539 string brep = shape->getBrep();
540 // translate_brep (brep, vdir, trep);
541 transfo_brep (brep, &gen_matrix, trep);
542 // Shape* tshape = new Shape (trep);
543 // tshape->setBounds (shape->getStart(), shape->getEnd());
544 // dest->addAssociation (tshape);
549 for (int nro=V_AMONT ; nro<=V_AVAL ; nro++)
551 Shape* shape = vo1->getAssociation ();
552 if (shape!=NULL && vd1->getAssociation ()==NULL)
554 string brep = shape->getBrep();
555 // translate_brep (brep, vdir, trep);
556 transfo_brep (brep, &gen_matrix, trep);
557 // Shape* tshape = new Shape (trep);
558 // vd1->setAssociation (tshape);
560 vo1 = orig->opposedVertex (vo1);
561 vd1 = dest->opposedVertex (vd1);
566 // ====================================================== assoResiduelle
567 void Elements::assoResiduelle ()
569 // int nbre = tab_vertex.size();
570 // for (int nv=0 ; nv<nbre ; nv++)
572 // geom_asso_point (tab_vertex [nv]);
575 // ====================================================== moveDisco
576 void Elements::moveDisco (Hexa* hexa)
580 hexa->getCenter (center);
581 matrix.defScale (center, 0.55);
583 int nbre = tab_vertex.size();
584 for (int nv=0 ; nv<nbre ; nv++)
586 matrix.perform (tab_vertex [nv]);