2 // C++ : Table d'hexaedres (Evol Versions 3)
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 "HexVector.hxx"
26 #include "HexVertex.hxx"
27 #include "HexEdge.hxx"
28 #include "HexHexa.hxx"
29 #include "HexMatrix.hxx"
30 #include "HexShape.hxx"
31 #include "HexGlobale.hxx"
37 void geom_create_circle (double* milieu, double rayon, double* normale,
38 double* base, string& brep);
39 void geom_create_sphere (double* milieu, double radius, string& brep);
41 void translate_brep (string& brep, double vdir[], string& trep);
42 void transfo_brep (string& brep, Matrix* matrice, string& trep);
43 void geom_asso_point (Vertex* node);
47 // ====================================================== makeRind
48 int Elements::makeRind (EnumGrid type, Vertex* center, Vector* vx, Vector* vz,
49 double radext, double radint, double radhole,
50 Vertex* plorig, double angle, int nr, int na, int nl)
53 int ier = controlRind (type, center, vx, vz, radext, radint, radhole,
54 plorig, angle, nr, na, nl, cyl_phi0, phi1);
58 resize (type, nr, na, nl);
62 cyl_radhole = radhole;
63 cyl_closed = type==GR_HEMISPHERIC || type==GR_RIND;
65 double theta = cyl_closed ? 2*M_PI : M_PI*angle/180;
66 cyl_dphi = (phi1-cyl_phi0)/ size_hz;
67 cyl_dtheta = theta / size_hy;
69 int nb_secteurs = cyl_closed ? size_vy-1 : size_vy;
71 for (int ny=0 ; ny<nb_secteurs ; ny++)
72 for (int nx=0 ; nx<size_vx ; nx++)
73 for (int nz=0 ; nz<size_vz ; nz++)
76 getCylPoint (nx, ny, nz, px, py, pz);
77 Vertex* node = el_root->addVertex (px, py, pz);
78 setVertex (node, nx, ny, nz);
81 for (int nx=0 ; nx<size_vx ; nx++)
82 for (int nz=0 ; nz<size_vz ; nz++)
84 Vertex* node = getVertexIJK (nx, 0, nz);
85 setVertex (node, nx, size_vy-1, nz);
88 transfoVertices (center, vx, vz);
90 assoCylinder (center, vz, angle);
93 // ====================================================== getCylPoint
94 int Elements::getCylPoint (int nr, int na, int nh, double& px, double& py,
97 if (grid_type == GR_CYLINDRIC)
99 px = cyl_radext * cos (na*cyl_dtheta);
100 py = cyl_radext * sin (na*cyl_dtheta);
101 pz = cyl_length * nh;
105 bool rind = (grid_type == GR_RIND || grid_type == GR_PART_RIND);
106 double sinphi = sin (cyl_phi0 + nh*cyl_dphi);
107 double cosphi = cos (cyl_phi0 + nh*cyl_dphi);
112 rayon = cyl_radint + nr*(cyl_radext-cyl_radint)/size_hx;
114 rayon = rayon*cosphi;
118 pz = cyl_radext*sinphi;
119 rayon = cyl_radhole + nr*(cyl_radext*cosphi - cyl_radhole)/size_hx;
120 rayon = std::max (cyl_radhole, rayon);
123 px = rayon * cos (na*cyl_dtheta);
124 py = rayon * sin (na*cyl_dtheta);
128 // ====================================================== controlRind
129 int Elements::controlRind (EnumGrid type, Vertex* cx, Vector* vx, Vector* vz,
130 double rext, double rint, double rhole,
131 Vertex* px, double angle,
132 int nrad, int nang, int nhaut,
133 double &phi0, double &phi1)
135 const double Epsil1 = 1e-6;
138 if (cx == NULL || vx == NULL || vz == NULL)
141 if (nrad<=0 || nang<=0 || nhaut<=0)
153 double nvx = vx->getNorm();
154 double nvz = vz->getNorm();
156 if (nvx < Epsil1 || nvz < Epsil1)
159 double alpha = asin (rhole/rext);
160 double beta = -M_PI*DEMI;
161 if (type==GR_HEMISPHERIC || type==GR_PART_SPHERIC)
167 double oh = ((px->getX() - cx->getX()) * vz->getDx()
168 + (px->getY() - cx->getY()) * vz->getDy()
169 + (px->getZ() - cx->getZ()) * vz->getDz()) / nvz;
173 beta = asin (oh/rext);
176 phi0 = std::max (alpha - M_PI*DEMI, beta);
177 phi1 = M_PI*DEMI - alpha;
180 // ====================================================== getHexas
181 int Elements::getHexas (Hexas& liste)
184 for (int nro = 0 ; nro<nbr_hexas ; nro++)
186 Hexa* cell = tab_hexa [nro];
187 if (cell!=NULL && cell->isValid())
188 liste.push_back (cell);
192 // ====================================================== assoCylinder
193 void Elements::assoCylinder (Vertex* ori, Vector* normal, double angle)
196 assoCylinders (ori, normal, angle, tangles);
198 // ====================================================== assoCylinders
199 void Elements::assoCylinders (Vertex* ori, Vector* normal, double angle,
200 RealVector& t_angles)
202 int na = t_angles.size();
203 bool regul = na == 0;
204 double alpha = angle/size_hy;
207 Real3 vk = { normal->getDx(), normal->getDy(), normal->getDz() };
210 for (int nz=0 ; nz<size_vz ; nz++)
212 for (int nx=0 ; nx<size_vx ; nx++)
214 Vertex* pm = getVertexIJK (nx, 0, nz);
215 Real3 om = { pm->getX() - ori->getX(),
216 pm->getY() - ori->getY(),
217 pm->getZ() - ori->getZ() };
219 double oh = prod_scalaire (om, vk);
222 for (int dd=dir_x; dd<=dir_z ; dd++)
224 ph [dd] = ori->getCoord(dd) + oh*vk[dd];
225 hm [dd] = pm ->getCoord(dd) - ph[dd];
226 rayon += hm[dd] * hm[dd];
229 rayon = sqrt (rayon);
230 geom_create_circle (ph, rayon, vk, hm, brep);
233 for (int ny=0 ; ny<size_hy ; ny++)
235 double beta = regul ? alpha : t_angles [ny];
238 Edge* edge = getEdgeJ (nx, ny, nz);
239 Shape* shape = new Shape (brep);
240 shape-> setBounds (pmin, pmax);
241 edge->addAssociation (shape);
246 // Association automatique des vertex non associes -> bph
247 // Traitement des faces spheriques
249 Real3 vi = { -vk[dir_x], -vk[dir_y], -vk[dir_z] };
250 Real3 po = { ori->getX(), ori->getY(), ori->getZ() };
254 case GR_HEMISPHERIC : // Pour l'exterieur
255 case GR_PART_SPHERIC :
256 assoRind (po, vi, size_vx-1);
258 case GR_PART_RIND : // Exterieur + interieur
260 assoRind (po, vi, 0);
261 assoRind (po, vi, size_vx-1);
266 assoResiduelle (); // Association des sommets residuels
268 // ====================================================== assoRind
269 // Association des meridiennes
270 // Creation sphere geometrique + association faces
271 void Elements::assoRind (double* ori, double* vi, int nx)
275 string c_brep, s_brep;
276 Shape shape (c_brep);
278 t_shape.push_back (&shape);
280 double radius = nx==0 ? cyl_radint : cyl_radext;
281 double paramin = (cyl_phi0 + M_PI/2) / (2*M_PI);
282 double paramax = paramin + size_hz*cyl_dphi/(2*M_PI);
284 paramin = std::max (paramin, 0.0);
285 paramax = std::min (paramax, 1.0);
286 geom_create_sphere (ori, radius, s_brep);
289 int nb_secteurs = cyl_closed ? size_vy-1 : size_vy;
291 for (int ny=0 ; ny<nb_secteurs ; ny++)
293 Vertex* pm = getVertexIJK (nx, ny, nz1);
294 Real3 vj = { pm->getX(), pm->getY(), pm->getZ() };
295 prod_vectoriel (vi, vj, vk);
296 double rayon = cyl_radint + nx*(cyl_radext-cyl_radint)/size_hx;
297 geom_create_circle (ori, rayon, vk, vi, c_brep);
298 shape.setBrep (c_brep);
299 shape.setBounds (paramin, paramax);
302 for (int nz=0 ; nz<size_hz ; nz++)
304 contour.push_back (getEdgeK (nx, ny, nz));
305 Quad* quad = getQuadJK (nx, ny, nz);
308 Shape* sshape = new Shape (s_brep);
309 quad->addAssociation (sshape);
312 cutAssociation (t_shape, contour, false);
315 // ====================================================== assoCircle
316 // ==== utilise pour les spheres carrees
317 void Elements::assoCircle (double* center, Edge* ed1, Edge* ed2)
319 Real3 oa, ob, normal;
320 Real3 pta, ptb, ptc, ptd;
323 // Les 2 edges dont les petits cotes d'un rectangle de rapport L/l=sqrt(2)
324 // Soit le cercle circonscrit a ce rectangle.
325 // * la largeur est balayee par l'angle alpha
326 // * la longueur par l'angle beta = pi -alpha
328 ed1->getVertex(V_AMONT)->getPoint (pta);
329 ed1->getVertex(V_AVAL )->getPoint (ptb);
330 ed2->getVertex(V_AMONT)->getPoint (ptc);
331 ed2->getVertex(V_AVAL )->getPoint (ptd);
333 double d1 = calc_distance (pta, ptc);
334 double d2 = calc_distance (pta, ptd);
338 ed2->getVertex(V_AMONT)->getPoint (ptd);
339 ed2->getVertex(V_AVAL )->getPoint (ptc);
342 calc_vecteur (center, pta, oa);
343 calc_vecteur (center, ptb, ob);
344 prod_vectoriel (oa, ob, normal);
345 double rayon = calc_norme (oa);
347 geom_create_circle (center, rayon, normal, oa, brep);
349 Shape* asso1 = new Shape (brep);
350 Shape* asso2 = new Shape (brep);
352 const double alpha = atan (sqrt(2)/2)/M_PI; // angle proche de 70.528 degres
353 asso1->setBounds (0, alpha);
354 asso2->setBounds (0.5, alpha + 0.5);
356 ed1->addAssociation (asso1);
357 ed2->addAssociation (asso2);
359 // ====================================================== assoSphere
360 void Elements::assoSphere (Vertex* ori, Edge* t_edge[], Quad* t_quad[])
362 Real3 center, sommet;
363 ori->getPoint(center);
365 assoCircle (center, t_edge [E_AC], t_edge [E_BD]);
366 assoCircle (center, t_edge [E_AD], t_edge [E_BC]);
367 assoCircle (center, t_edge [E_AE], t_edge [E_BF]);
368 assoCircle (center, t_edge [E_AF], t_edge [E_BE]);
369 assoCircle (center, t_edge [E_CE], t_edge [E_DF]);
370 assoCircle (center, t_edge [E_CF], t_edge [E_DE]);
372 t_edge[E_AC]->getVertex(V_AMONT)->getPoint (sommet);
373 double radius = calc_distance (center, sommet);;
376 geom_create_sphere (center, radius, brep);
378 for (int nf=0 ; nf < HQ_MAXI ; nf++)
380 Shape* shape = new Shape (brep);
381 t_quad [nf]->addAssociation (shape);
384 assoResiduelle (); // Association des sommets residuels
386 // ====================================================== makeSphericalGrid
387 int Elements::makeSphericalGrid (Vertex* c, double rayon, int nb, double k)
389 resize (GR_SPHERIC, nb);
393 else if (rayon <=ZEROR)
396 Vertex* i_node [HV_MAXI]; // Les noeuds de l'hexa englobant
397 Edge* i_edge [HE_MAXI]; // Les noeuds de l'hexa englobant
398 Quad* i_quad [HQ_MAXI]; // Les noeuds de l'hexa englobant
400 for (int nro=0 ; nro<HV_MAXI; nro++)
402 double dx = glob->CoordVertex (nro, dir_x) * rayon;
403 double dy = glob->CoordVertex (nro, dir_y) * rayon;
404 double dz = glob->CoordVertex (nro, dir_z) * rayon;
406 i_node [nro] = el_root->addVertex (c->getX ()+dx, c->getY ()+dy,
410 for (int nro=0 ; nro<HE_MAXI; nro++)
412 int v1 = glob->EdgeVertex (nro, V_AMONT);
413 int v2 = glob->EdgeVertex (nro, V_AVAL);
414 i_edge[nro] = newEdge (i_node[v1], i_node[v2]);
418 char nm0[8], nm1 [8], nm2 [8];
419 printf (" %2d : %s = %s = [%s, %s] = [%d,%d] = [%s,%s]\n", nro,
420 glob->namofHexaEdge(nro), i_edge[nro]->getName(nm0),
421 glob->namofHexaVertex(v1), glob->namofHexaVertex(v2), v1, v2,
422 i_node[v1]->getName(nm1), i_node[v2]->getName(nm2));
426 for (int nro=0 ; nro<HQ_MAXI; nro++)
427 i_quad[nro] = newQuad (i_edge[glob->QuadEdge (nro, E_A)],
428 i_edge[glob->QuadEdge (nro, E_B)],
429 i_edge[glob->QuadEdge (nro, E_C)],
430 i_edge[glob->QuadEdge (nro, E_D)]);
432 tab_hexa.push_back (newHexa (i_quad[Q_A], i_quad[Q_B], i_quad[Q_C],
433 i_quad[Q_D], i_quad[Q_E], i_quad[Q_F]));
436 for (int niv=0; niv<gr_rayon ; niv++)
438 double lambda0 = lambda;
441 addStrate (i_quad, i_edge, i_node, c, lambda/lambda0);
444 assoSphere (c, i_edge, i_quad);
447 // ==================================================== propagateAssociation
448 int Elements::propagateAssociation (Edge* orig, Edge* dest, Edge* dir)
450 if (revo_lution || orig==NULL || dest==NULL || dir==NULL)
453 const Shapes& tab_shapes = orig->getAssociations ();
454 const Shapes& tab_dest = dest->getAssociations ();
455 int nbdest = tab_dest.size();
456 int nbshapes = tab_shapes.size();
457 bool on_edge = nbshapes!=0 && nbdest==0;
459 Vertex* vo1 = orig->commonVertex (dir);
460 Vertex* vd1 = dest->commonVertex (dir);
461 Vertex* vo2 = orig->opposedVertex (vo1);
462 Vertex* vd2 = dest->opposedVertex (vd1);
464 if (vo1==NULL || vd1==NULL)
468 Real3 pa, pb, vdir1, vdir2;
469 calc_vecteur (vo1->getPoint (pa), vd1->getPoint (pb), vdir1);
470 calc_vecteur (vo2->getPoint (pa), vd2->getPoint (pb), vdir2);
472 double dd = calc_distance (vdir1, vdir2);
473 bool para = dd < 1.0e-3;
477 for (int nro=0 ; nro<nbshapes ; nro++)
479 Shape* shape = tab_shapes[nro];
482 string brep = shape->getBrep();
483 translate_brep (brep, vdir1, trep);
484 Shape* tshape = new Shape (trep);
485 tshape->setBounds (shape->debut, shape->fin);
486 dest->addAssociation (tshape);
491 double* vdir = vdir1;
492 for (int nro=V_AMONT ; nro<=V_AVAL ; nro++)
494 Shape* shape = vo1->getAssociation ();
495 if (shape!=NULL && vd1->getAssociation ()==NULL)
497 string brep = shape->getBrep();
498 translate_brep (brep, vdir, trep);
499 Shape* tshape = new Shape (trep);
500 vd1->setAssociation (tshape);
509 // ==================================================== prismAssociation
510 int Elements::prismAssociation (Edge* orig, Edge* dest, int nh, Edge* dir)
512 if (orig==NULL || dest==NULL || dir==NULL)
517 const Shapes& tab_shapes = orig->getAssociations ();
518 const Shapes& tab_dest = dest->getAssociations ();
519 int nbdest = tab_dest.size();
520 int nbshapes = tab_shapes.size();
521 bool on_edge = nbshapes>0 && nbdest==0;
523 Vertex* vo1 = orig->commonVertex (dir);
524 Vertex* vd1 = dest->commonVertex (dir);
526 if (vo1==NULL || vd1==NULL)
530 Real3 porig, pdest, vdir;
531 vo1->getPoint (porig);
532 vd1->getPoint (pdest);
533 calc_vecteur (porig, pdest, vdir);
537 for (int nro=0 ; nro<nbshapes ; nro++)
539 Shape* shape = tab_shapes[nro];
542 string brep = shape->getBrep();
543 // translate_brep (brep, vdir, trep);
544 transfo_brep (brep, &gen_matrix, trep);
545 Shape* tshape = new Shape (trep);
546 tshape->setBounds (shape->debut, shape->fin);
547 dest->addAssociation (tshape);
552 for (int nro=V_AMONT ; nro<=V_AVAL ; nro++)
554 Shape* shape = vo1->getAssociation ();
555 if (shape!=NULL && vd1->getAssociation ()==NULL)
557 string brep = shape->getBrep();
558 // translate_brep (brep, vdir, trep);
559 transfo_brep (brep, &gen_matrix, trep);
560 Shape* tshape = new Shape (trep);
561 vd1->setAssociation (tshape);
563 vo1 = orig->opposedVertex (vo1);
564 vd1 = dest->opposedVertex (vd1);
568 // ====================================================== assoResiduelle
569 void Elements::assoResiduelle ()
571 int nbre = tab_vertex.size();
572 for (int nv=0 ; nv<nbre ; nv++)
574 geom_asso_point (tab_vertex [nv]);
577 // ====================================================== moveDisco
578 void Elements::moveDisco (Hexa* hexa)
582 hexa->getCenter (center);
583 matrix.defScale (center, 0.55);
585 int nbre = tab_vertex.size();
586 for (int nv=0 ; nv<nbre ; nv++)
588 matrix.perform (tab_vertex [nv]);