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"
34 #include "HexAssoEdge.hxx"
40 // static bool db=false;
42 // ====================================================== makeRind
43 int Elements::makeRind (EnumGrid type, Vertex* center, Vector* vx, Vector* vz,
44 double radext, double radint, double radhole,
45 Vertex* plorig, double angle, int nr, int na, int nl)
48 int ier = controlRind (type, center, vx, vz, radext, radint, radhole,
49 plorig, angle, nr, na, nl, cyl_phi0, phi1);
53 resize (type, nr, na, nl);
57 cyl_radhole = radhole;
58 cyl_closed = type==GR_HEMISPHERIC || type==GR_RIND;
60 double theta = cyl_closed ? 2*M_PI : M_PI*angle/180;
61 cyl_dphi = (phi1-cyl_phi0)/ size_hz;
62 cyl_dtheta = theta / size_hy;
64 int nb_secteurs = cyl_closed ? size_vy-1 : size_vy;
66 for (int ny=0 ; ny<nb_secteurs ; ny++)
67 for (int nx=0 ; nx<size_vx ; nx++)
68 for (int nz=0 ; nz<size_vz ; nz++)
71 getCylPoint (nx, ny, nz, px, py, pz);
72 Vertex* node = el_root->addVertex (px, py, pz);
73 setVertex (node, nx, ny, nz);
76 for (int nx=0 ; nx<size_vx ; nx++)
77 for (int nz=0 ; nz<size_vz ; nz++)
79 Vertex* node = getVertexIJK (nx, 0, nz);
80 setVertex (node, nx, size_vy-1, nz);
83 transfoVertices (center, vx, vz);
85 assoCylinder (center, vz, angle);
88 // ====================================================== getCylPoint
89 int Elements::getCylPoint (int nr, int na, int nh, double& px, double& py,
92 if (grid_type == GR_CYLINDRIC)
94 px = cyl_radext * cos (na*cyl_dtheta);
95 py = cyl_radext * sin (na*cyl_dtheta);
100 bool rind = (grid_type == GR_RIND || grid_type == GR_PART_RIND);
101 double sinphi = sin (cyl_phi0 + nh*cyl_dphi);
102 double cosphi = cos (cyl_phi0 + nh*cyl_dphi);
107 rayon = cyl_radint + nr*(cyl_radext-cyl_radint)/size_hx;
109 rayon = rayon*cosphi;
113 pz = cyl_radext*sinphi;
114 rayon = cyl_radhole + nr*(cyl_radext*cosphi - cyl_radhole)/size_hx;
115 rayon = std::max (cyl_radhole, rayon);
118 px = rayon * cos (na*cyl_dtheta);
119 py = rayon * sin (na*cyl_dtheta);
123 // ====================================================== controlRind
124 int Elements::controlRind (EnumGrid type, Vertex* cx, Vector* vx, Vector* vz,
125 double rext, double rint, double rhole,
126 Vertex* px, double angle,
127 int nrad, int nang, int nhaut,
128 double &phi0, double &phi1)
130 const double Epsil1 = 1e-6;
133 if (cx == NULL || vx == NULL || vz == NULL)
136 if (nrad<=0 || nang<=0 || nhaut<=0)
148 double nvx = vx->getNorm();
149 double nvz = vz->getNorm();
151 if (nvx < Epsil1 || nvz < Epsil1)
154 double alpha = asin (rhole/rext);
155 double beta = -M_PI*DEMI;
156 if (type==GR_HEMISPHERIC || type==GR_PART_SPHERIC)
162 double oh = ((px->getX() - cx->getX()) * vz->getDx()
163 + (px->getY() - cx->getY()) * vz->getDy()
164 + (px->getZ() - cx->getZ()) * vz->getDz()) / nvz;
168 beta = asin (oh/rext);
171 phi0 = std::max (alpha - M_PI*DEMI, beta);
172 phi1 = M_PI*DEMI - alpha;
175 // ====================================================== getHexas
176 int Elements::getHexas (Hexas& liste)
179 for (int nro = 0 ; nro<nbr_hexas ; nro++)
181 Hexa* cell = tab_hexa [nro];
182 if (cell!=NULL && cell->isValid())
183 liste.push_back (cell);
187 // ====================================================== assoCylinder
188 void Elements::assoCylinder (Vertex* ori, Vector* normal, double angle)
190 if (normal==NULL || ori==NULL)
195 normal->getCoord (vz);
196 ori ->getPoint (center);
197 assoCylinders (center, vz, angle, tangles);
199 // ====================================================== assoCylinders
200 void Elements::assoCylinders (Vertex* ori, Vector* normal, double angle,
201 RealVector& t_angles)
203 if (normal==NULL || ori==NULL)
207 normal->getCoord (vz);
208 ori ->getPoint (center);
209 assoCylinders (center, vz, angle, t_angles);
211 // ====================================================== assoCylinders
212 void Elements::assoCylinders (double* ori, double* vk, double angle,
213 RealVector& t_angles)
216 sprintf (name, "grid_%02d", el_id);
217 NewShape* geom = el_root->addShape (name, SH_CYLINDER);
221 // Real3 vk = { normal->getDx(), normal->getDy(), normal->getDz() };
222 // normer_vecteur (vk);
224 for (int nz=0 ; nz<size_vz ; nz++)
226 for (int nx=0 ; nx<size_vx ; nx++)
228 Vertex* pm = getVertexIJK (nx, 0, nz);
229 Real3 om = { pm->getX() - ori[dir_x],
230 pm->getY() - ori[dir_y],
231 pm->getZ() - ori[dir_z] };
233 double oh = prod_scalaire (om, vk);
236 for (int dd=dir_x; dd<=dir_z ; dd++)
238 ph [dd] = ori[dd] + oh*vk[dd];
239 hm [dd] = pm ->getCoord(dd) - ph[dd];
240 rayon += hm[dd] * hm[dd];
243 rayon = sqrt (rayon);
244 int subid = geom->addCircle (ph, rayon, vk, hm);
246 for (int ny=0 ; ny<size_hy ; ny++)
248 double pmin = t_angles [ny]/360;
249 double pmax = t_angles [ny+1]/360;
250 Edge* edge = getEdgeJ (nx, ny, nz);
251 geom->addAssociation (edge, subid, pmin, pmax);
252 cout << " assoCylinders : ny= " << ny << ", nz= " << nz
253 << " param = (" << pmin << ", " << pmax << ")\n";
257 // Association automatique des vertex non associes -> bph
258 // Traitement des faces spheriques
260 Real3 vi = { -vk[dir_x], -vk[dir_y], -vk[dir_z] };
264 case GR_HEMISPHERIC : // Pour l'exterieur
265 case GR_PART_SPHERIC :
266 assoRind (ori, vi, size_vx-1, geom);
268 case GR_PART_RIND : // Exterieur + interieur
270 assoRind (ori, vi, 0, geom);
271 assoRind (ori, vi, size_vx-1, geom);
278 // ====================================================== calcul_param
279 double calcul_param (double* ori, double* vi, Vertex* vert)
282 vert->getPoint (pnt);
283 calc_vecteur (ori, pnt, vj);
285 double kos = prod_scalaire (vi, vj);
286 double alpha = acos (kos) / (2*M_PI);
289 // ====================================================== assoRind
290 // Association des meridiennes
291 // Creation sphere geometrique + association faces
292 void Elements::assoRind (double* ori, double* vi, int nx, NewShape* geom)
295 double radius = nx==0 ? cyl_radint : cyl_radext;
296 int sphid = geom->addSphere (ori, radius);
299 int nb_secteurs = cyl_closed ? size_vy-1 : size_vy;
301 for (int ny=0 ; ny<nb_secteurs ; ny++)
303 Vertex* pm = getVertexIJK (nx, ny, nz1);
304 Real3 vj = { pm->getX(), pm->getY(), pm->getZ() };
305 prod_vectoriel (vi, vj, vk);
306 double rayon = cyl_radint + nx*(cyl_radext-cyl_radint)/size_hx;
307 int subid = geom->addCircle (ori, rayon, vk, vi);
309 for (int nz=0 ; nz<size_hz ; nz++)
311 Quad* quad = getQuadJK (nx, ny, nz);
312 Edge* edge = getEdgeK (nx, ny, nz);
313 Vertex* nd1 = edge->getVertex (V_AMONT);
314 Vertex* nd2 = edge->getVertex (V_AVAL);
315 double pmin = calcul_param (ori, vi, nd1);
316 double pmax = calcul_param (ori, vi, nd2);
317 cout << " assoRind : ny= " << ny << ", nz= " << nz
318 << " param = (" << pmin << ", " << pmax << ")\n";
320 geom->addAssociation (edge, subid, pmin, pmax);
321 geom->addAssociation (quad, sphid);
325 // ====================================================== assoCircle
326 // ==== utilise pour les spheres carrees
327 void Elements::assoCircle (double* center, Edge* ed1, Edge* ed2, NewShape* geom)
329 Real3 oa, ob, normal;
330 Real3 pta, ptb, ptc, ptd;
333 // Les 2 edges dont les petits cotes d'un rectangle de rapport L/l=sqrt(2)
334 // Soit le cercle circonscrit a ce rectangle.
335 // * la largeur est balayee par l'angle alpha
336 // * la longueur par l'angle beta = pi -alpha
338 ed1->getVertex(V_AMONT)->getPoint (pta);
339 ed1->getVertex(V_AVAL )->getPoint (ptb);
340 ed2->getVertex(V_AMONT)->getPoint (ptc);
341 ed2->getVertex(V_AVAL )->getPoint (ptd);
343 double d1 = calc_distance (pta, ptc);
344 double d2 = calc_distance (pta, ptd);
348 ed2->getVertex(V_AMONT)->getPoint (ptd);
349 ed2->getVertex(V_AVAL )->getPoint (ptc);
352 calc_vecteur (center, pta, oa);
353 calc_vecteur (center, ptb, ob);
354 prod_vectoriel (oa, ob, normal);
356 double rayon = calc_norme (oa);
357 int subid = geom->addCircle (center, rayon, normal, oa);
359 const double alpha = atan (sqrt(2)/2)/M_PI; // angle proche de 70.528 degres
360 // asso1->setBounds (0, alpha);
361 // asso2->setBounds (0.5, alpha + 0.5);
363 geom->addAssociation (ed1, subid, 0.0, alpha);
364 geom->addAssociation (ed2, subid, 0.5, alpha+0.5);
366 // ====================================================== assoSphere
367 void Elements::assoSphere (double* center, Edge* t_edge[], Quad* t_quad[])
370 sprintf (name, "grid_%02d", el_id);
371 NewShape* geom = el_root->addShape (name, SH_SPHERE);
372 geom -> openShape ();
376 assoCircle (center, t_edge [E_AC], t_edge [E_BD], geom);
377 assoCircle (center, t_edge [E_AD], t_edge [E_BC], geom);
378 assoCircle (center, t_edge [E_AE], t_edge [E_BF], geom);
379 assoCircle (center, t_edge [E_AF], t_edge [E_BE], geom);
380 assoCircle (center, t_edge [E_CE], t_edge [E_DF], geom);
381 assoCircle (center, t_edge [E_CF], t_edge [E_DE], geom);
383 t_edge[E_AC]->getVertex(V_AMONT)->getPoint (sommet);
384 double radius = calc_distance (center, sommet);
386 int sphid = geom -> addSphere (center, radius);
389 for (int nf=0 ; nf < HQ_MAXI ; nf++)
390 t_quad[nf]->addAssociation (geom, sphid);
392 // ==================================================== propagateAssociation
393 int Elements::propagateAssociation (Edge* orig, Edge* dest, Edge* dir)
397 if (revo_lution || orig==NULL || dest==NULL || dir==NULL)
400 const Shapes& tab_shapes = orig->getAssociations ();
401 const Shapes& tab_dest = dest->getAssociations ();
402 int nbdest = tab_dest.size();
403 int nbshapes = tab_shapes.size();
404 bool on_edge = nbshapes!=0 && nbdest==0;
406 Vertex* vo1 = orig->commonVertex (dir);
407 Vertex* vd1 = dest->commonVertex (dir);
408 Vertex* vo2 = orig->opposedVertex (vo1);
409 Vertex* vd2 = dest->opposedVertex (vd1);
411 if (vo1==NULL || vd1==NULL)
415 Real3 pa, pb, vdir1, vdir2;
416 calc_vecteur (vo1->getPoint (pa), vd1->getPoint (pb), vdir1);
417 calc_vecteur (vo2->getPoint (pa), vd2->getPoint (pb), vdir2);
419 double dd = calc_distance (vdir1, vdir2);
420 bool para = dd < 1.0e-3;
424 for (int nro=0 ; nro<nbshapes ; nro++)
426 Shape* shape = tab_shapes[nro];
429 string brep = shape->getBrep();
430 translate_brep (brep, vdir1, trep);
431 // Shape* tshape = new Shape (trep);
432 // tshape->setBounds (shape->getStart(), shape->getEnd());
433 // dest->addAssociation (tshape);
438 double* vdir = vdir1;
439 for (int nro=V_AMONT ; nro<=V_AVAL ; nro++)
441 Shape* shape = vo1->getAssociation ();
442 if (shape!=NULL && vd1->getAssociation ()==NULL)
444 string brep = shape->getBrep();
445 translate_brep (brep, vdir, trep);
446 // Shape* tshape = new Shape (trep);
447 // vd1->setAssociation (tshape);
456 // ==================================================== prismAssociation
457 int Elements::prismAssociation (Edge* lorig, Edge* ldest, int nh)
459 if (lorig==NULL || ldest==NULL)
463 int nbshapes = lorig->countAssociation ();
467 NewShape* new_shape = getShape ();
468 for (int nro=0 ; nro<nbshapes ; nro++)
470 AssoEdge* asso = lorig->getAssociation (nro);
471 EdgeShape* shape = asso ->getEdgeShape();
472 double p1 = asso ->getStart();
473 double p2 = asso ->getEnd ();
477 sprintf (name, "0x%lx#%d", (unsigned long) shape, nh);
478 map<string,int>::iterator iter = map_shape.find (name);
479 if (iter != map_shape.end())
480 subid = iter->second;
482 subid = map_shape[name] = new_shape->transfoShape (cum_matrix, shape);
484 new_shape->addAssociation (ldest, subid, p1, p2);
485 printf (" prismAsso : %s -> %s(%d) = %s [ %g , %g]\n",
486 lorig->getName(), ldest->getName(), nh, name, p1, p2);
490 // ====================================================== assoResiduelle
491 void Elements::assoResiduelle ()
493 // int nbre = tab_vertex.size();
494 // for (int nv=0 ; nv<nbre ; nv++)
496 // geom_asso_point (tab_vertex [nv]);
499 // ====================================================== moveDisco
500 void Elements::moveDisco (Hexa* hexa)
504 hexa->getCenter (center);
505 matrix.defScale (center, 0.55);
507 int nbre = tab_vertex.size();
508 for (int nv=0 ; nv<nbre ; nv++)
510 matrix.perform (tab_vertex [nv]);
513 // =================================================== getStrate
514 Hexa* Elements::getStrate (int couche, EnumHQuad nroface)
517 int nro = couche <= 0 ? 0 : (couche-1)*HQ_MAXI + nroface + 1;
519 if (nbr_hexas==0 || nro >= nbr_hexas)
522 cell = tab_hexa [nro];
526 // ============================================================ setHexa
527 void Elements::setHexa (Hexa* elt, int nro)
529 if (nro >=0 && nro < nbr_hexas)
530 tab_hexa [nro] = elt;
532 // ============================================================ setQuad
533 void Elements::setQuad (Quad* elt, int nro)
535 if (nro >=0 && nro < nbr_quads)
536 tab_quad [nro] = elt;
538 // ============================================================ setEdge
539 void Elements::setEdge (Edge* elt, int nro)
541 if (nro >=0 && nro < nbr_edges)
542 tab_edge [nro] = elt;
544 // ============================================================ setVertex
545 void Elements::setVertex (Vertex* elt, int nro)
547 if (nro >=0 && nro < nbr_vertex)
548 tab_vertex [nro] = elt;
550 // -----------------------------------------------------------------------
551 // -----------------------------------------------------------------------
552 // ============================================================ getHexa
553 Hexa* Elements::getHexa (int nro)
555 DumpStart ("getHexa", nro);
558 int nombre=tab_hexa.size();
559 // if (nro >=0 && nro < nbr_hexas && el_status == HOK Abu 2010/05/06
560 if (nro >=0 && nro < nombre && el_status == HOK
561 && tab_hexa [nro] != NULL && tab_hexa [nro]->isValid())
562 elt = tab_hexa [nro];
567 // ============================================================ getQuad
568 Quad* Elements::getQuad (int nro)
570 DumpStart ("getQuad", nro);
573 if (nro >=0 && nro < nbr_quads && el_status == HOK
574 && tab_quad [nro] != NULL && tab_quad [nro]->isValid())
575 elt = tab_quad [nro];
580 // ============================================================ getEdge
581 Edge* Elements::getEdge (int nro)
583 DumpStart ("getEdge", nro);
586 if (nro >=0 && nro < nbr_edges && el_status == HOK
587 && tab_edge [nro] != NULL && tab_edge [nro]->isValid())
588 elt = tab_edge [nro];
593 // ============================================================ getVertex
594 Vertex* Elements::getVertex (int nro)
596 DumpStart ("getVertex", nro);
599 if (nro >=0 && nro < nbr_vertex && el_status == HOK
600 && tab_vertex [nro] != NULL && tab_vertex [nro]->isValid())
601 elt = tab_vertex [nro];
606 // ============================================================ indVertex
607 int Elements::indVertex (int nquad, int nsommet, int nh)
609 int nro = nsommet + QUAD4*nquad + nbr_orig*QUAD4*(nh+1);
612 // ============================================================ nroVertex
613 int Elements::nroVertex (int nquad, int nsommet, int nh)
615 int nro = nsommet + QUAD4*nquad + nbr_orig*QUAD4*nh;
618 // ============================================================ indVertex
619 int Elements::indVertex (int ref_edge, int nh)
621 int nro = ref_edge + nbr_orig*QUAD4*nh;
624 // ============================================================ nroEdgeH
625 int Elements::nroEdgeH (int nvertex)
627 return QUAD4*nbr_orig*gr_hauteur + nvertex;
629 // ============================================================ nroEdgeH
630 int Elements::nroEdgeH (int nquad, int nsommet, int nh)
632 return QUAD4*nbr_orig*gr_hauteur + indVertex (nquad, nsommet, nh);
634 // ============================================================ nroHexa
635 int Elements::nroHexa (int nquad, int nh)
637 int nro = gr_hauteur*nquad + nh;
641 // ============================================================ addHexa
642 void Elements::addHexa (Hexa* element)
644 tab_hexa.push_back (element);
647 // ============================================================ addQuad
648 void Elements::addQuad (Quad* element)
650 tab_quad.push_back (element);
653 // ============================================================ addEdge
654 void Elements::addEdge (Edge* element)
656 tab_edge.push_back (element);
659 // ============================================================ addVertex
660 void Elements::addVertex (Vertex* element)
662 tab_vertex.push_back (element);
665 // ============================================================ findHexa
666 int Elements::findHexa (Hexa* element)
668 int nbre = tab_hexa.size();
669 for (int nro=0 ; nro<nbre ; nro++)
670 if (tab_hexa [nro] == element)
674 // ============================================================ findQuad
675 int Elements::findQuad (Quad* element)
677 int nbre = tab_quad.size();
678 for (int nro=0 ; nro<nbre ; nro++)
679 if (tab_quad [nro] == element)
683 // ============================================================ findEdge
684 int Elements::findEdge (Edge* element)
686 int nbre = tab_edge.size();
687 for (int nro=0 ; nro<nbre ; nro++)
688 if (tab_edge [nro] == element)
692 // ============================================================ findVertex
693 int Elements::findVertex (Vertex* element)
695 int nbre = tab_vertex.size();
696 for (int nro=0 ; nro<nbre ; nro++)
697 if (tab_vertex [nro] == element)
701 // ========================================================= saveVtk (avec nro)
702 int Elements::saveVtk (cpchar radical, int &nro)
705 sprintf (num, "%d", nro);
708 string filename = radical;
711 int ier = saveVtk (filename.c_str());