1 // C++ : Table d'hexaedres (Evol Versions 3)
3 // Copyright (C) 2009-2015 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, or (at your option) any later version.
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"
42 // ====================================================== getCylPoint
43 int Elements::getCylPoint (int nr, int na, int nh, double& px, double& py,
46 if (grid_type == GR_CYLINDRIC)
48 px = cyl_radext * cos (na*cyl_dtheta);
49 py = cyl_radext * sin (na*cyl_dtheta);
54 bool rind = (grid_type == GR_RIND || grid_type == GR_PART_RIND);
55 double sinphi = sin (cyl_phi0 + nh*cyl_dphi);
56 double cosphi = cos (cyl_phi0 + nh*cyl_dphi);
61 rayon = cyl_radint + nr*(cyl_radext-cyl_radint)/size_hx;
67 pz = cyl_radext*sinphi;
68 rayon = cyl_radhole + nr*(cyl_radext*cosphi - cyl_radhole)/size_hx;
69 rayon = std::max (cyl_radhole, rayon);
72 px = rayon * cos (na*cyl_dtheta);
73 py = rayon * sin (na*cyl_dtheta);
77 // ====================================================== controlRind
78 int Elements::controlRind (EnumGrid type, Vertex* cx, Vector* vx, Vector* vz,
79 double rext, double rint, double rhole,
80 Vertex* px, double angle,
81 int nrad, int nang, int nhaut,
82 double &phi0, double &phi1)
84 const double Epsil1 = 1e-6;
87 if (cx == NULL || vx == NULL || vz == NULL)
90 if (nrad<=0 || nang<=0 || nhaut<=0)
102 double nvx = vx->getNorm();
103 double nvz = vz->getNorm();
105 if (nvx < Epsil1 || nvz < Epsil1)
108 double alpha = asin (rhole/rext);
109 double beta = -M_PI*DEMI;
110 if (type==GR_HEMISPHERIC || type==GR_PART_SPHERIC)
116 double oh = ((px->getX() - cx->getX()) * vz->getDx()
117 + (px->getY() - cx->getY()) * vz->getDy()
118 + (px->getZ() - cx->getZ()) * vz->getDz()) / nvz;
122 beta = asin (oh/rext);
125 phi0 = std::max (alpha - M_PI*DEMI, beta);
126 phi1 = M_PI*DEMI - alpha;
129 // ====================================================== getHexas
130 int Elements::getHexas (Hexas& liste)
133 for (int nro = 0 ; nro<nbr_hexas ; nro++)
135 Hexa* cell = tab_hexa [nro];
136 if (cell!=NULL && cell->isValid())
137 liste.push_back (cell);
141 // ====================================================== assoCylinder
142 void Elements::assoCylinder (Vertex* ori, Vector* normal, double angle)
144 if (normal==NULL || ori==NULL)
149 normal->getCoord (vz);
150 ori ->getPoint (center);
151 assoCylinders (center, vz, angle, tangles);
153 // ====================================================== assoCylinders
154 void Elements::assoCylinders (Vertex* ori, Vector* normal, double angle,
155 RealVector& t_angles)
157 if (normal==NULL || ori==NULL)
161 normal->getCoord (vz);
162 ori ->getPoint (center);
163 assoCylinders (center, vz, angle, t_angles);
165 // ====================================================== assoCylinders
166 void Elements::assoCylinders (double* ori, double* vk, double angle,
167 RealVector& t_angles)
170 sprintf (name, "grid_%02d", el_id);
171 NewShape* geom = el_root->addShape (name, SH_CYLINDER);
175 // Real3 vk = { normal->getDx(), normal->getDy(), normal->getDz() };
176 // normer_vecteur (vk);
178 for (int nz=0 ; nz<size_vz ; nz++)
180 for (int nx=0 ; nx<size_vx ; nx++)
182 Vertex* pm = getVertexIJK (nx, 0, nz);
183 Real3 om = { pm->getX() - ori[dir_x],
184 pm->getY() - ori[dir_y],
185 pm->getZ() - ori[dir_z] };
187 double oh = prod_scalaire (om, vk);
190 for (int dd=dir_x; dd<=dir_z ; dd++)
192 ph [dd] = ori[dd] + oh*vk[dd];
193 hm [dd] = pm ->getCoord(dd) - ph[dd];
194 rayon += hm[dd] * hm[dd];
197 rayon = sqrt (rayon);
198 int subid = geom->addCircle (ph, rayon, vk, hm);
200 for (int ny=0 ; ny<size_hy ; ny++)
202 double pmin = t_angles [ny]/360;
203 double pmax = t_angles [ny+1]/360;
204 Edge* edge = getEdgeJ (nx, ny, nz);
205 geom->addAssociation (edge, subid, pmin, pmax);
206 if (db) cout << " assoCylinders : ny= " << ny
207 << ", nz= " << nz << " param = ("
208 << pmin << ", " << pmax << ")\n";
212 // Association automatique des vertex non associes -> bph
213 // Traitement des faces spheriques
215 Real3 vi = { -vk[dir_x], -vk[dir_y], -vk[dir_z] };
219 case GR_HEMISPHERIC : // Pour l'exterieur
220 case GR_PART_SPHERIC :
221 assoRind (ori, vi, size_vx-1, geom);
223 case GR_PART_RIND : // Exterieur + interieur
225 assoRind (ori, vi, 0, geom);
226 assoRind (ori, vi, size_vx-1, geom);
233 // ====================================================== calcul_param
234 double calcul_param (double* ori, double* vi, Vertex* vert)
237 vert->getPoint (pnt);
238 calc_vecteur (ori, pnt, vj);
240 double kos = prod_scalaire (vi, vj);
241 double alpha = acos (kos) / (2*M_PI);
244 // ====================================================== assoRind
245 // Association des meridiennes
246 // Creation sphere geometrique + association faces
247 void Elements::assoRind (double* ori, double* vi, int nx, NewShape* geom)
250 double radius = nx==0 ? cyl_radint : cyl_radext;
251 int sphid = geom->addSphere (ori, radius);
254 int nb_secteurs = cyl_closed ? size_vy-1 : size_vy;
256 for (int ny=0 ; ny<nb_secteurs ; ny++)
258 Vertex* pm = getVertexIJK (nx, ny, nz1);
259 Real3 vj = { pm->getX(), pm->getY(), pm->getZ() };
260 prod_vectoriel (vi, vj, vk);
261 double rayon = cyl_radint + nx*(cyl_radext-cyl_radint)/size_hx;
262 int subid = geom->addCircle (ori, rayon, vk, vi);
264 for (int nz=0 ; nz<size_hz ; nz++)
266 Quad* quad = getQuadJK (nx, ny, nz);
267 Edge* edge = getEdgeK (nx, ny, nz);
268 Vertex* nd1 = edge->getVertex (V_AMONT);
269 Vertex* nd2 = edge->getVertex (V_AVAL);
270 double pmin = calcul_param (ori, vi, nd1);
271 double pmax = calcul_param (ori, vi, nd2);
272 cout << " assoRind : ny= " << ny << ", nz= " << nz
273 << " param = (" << pmin << ", " << pmax << ")\n";
275 geom->addAssociation (edge, subid, pmin, pmax);
276 geom->addAssociation (quad, sphid);
280 // ====================================================== assoCircle
281 // ==== utilise pour les spheres carrees
282 void Elements::assoCircle (double* center, Edge* ed1, Edge* ed2, NewShape* geom)
284 Real3 oa, ob, normal;
285 Real3 pta, ptb, ptc, ptd;
288 // Les 2 edges dont les petits cotes d'un rectangle de rapport L/l=sqrt(2)
289 // Soit le cercle circonscrit a ce rectangle.
290 // * la largeur est balayee par l'angle alpha
291 // * la longueur par l'angle beta = pi -alpha
293 ed1->getVertex(V_AMONT)->getPoint (pta);
294 ed1->getVertex(V_AVAL )->getPoint (ptb);
295 ed2->getVertex(V_AMONT)->getPoint (ptc);
296 ed2->getVertex(V_AVAL )->getPoint (ptd);
298 double d1 = calc_distance (pta, ptc);
299 double d2 = calc_distance (pta, ptd);
303 ed2->getVertex(V_AMONT)->getPoint (ptd);
304 ed2->getVertex(V_AVAL )->getPoint (ptc);
307 calc_vecteur (center, pta, oa);
308 calc_vecteur (center, ptb, ob);
309 prod_vectoriel (oa, ob, normal);
311 double rayon = calc_norme (oa);
312 int subid = geom->addCircle (center, rayon, normal, oa);
314 const double alpha = atan (sqrt(2.)/2)/M_PI; // angle proche de 70.528 degres
315 // asso1->setBounds (0, alpha);
316 // asso2->setBounds (0.5, alpha + 0.5);
318 geom->addAssociation (ed1, subid, 0.0, alpha);
319 geom->addAssociation (ed2, subid, 0.5, alpha+0.5);
321 // ====================================================== assoSphere
322 void Elements::assoSphere (double* center, Edge* t_edge[], Quad* t_quad[])
325 sprintf (name, "grid_%02d", el_id);
326 NewShape* geom = el_root->addShape (name, SH_SPHERE);
327 geom -> openShape ();
331 assoCircle (center, t_edge [E_AC], t_edge [E_BD], geom);
332 assoCircle (center, t_edge [E_AD], t_edge [E_BC], geom);
333 assoCircle (center, t_edge [E_AE], t_edge [E_BF], geom);
334 assoCircle (center, t_edge [E_AF], t_edge [E_BE], geom);
335 assoCircle (center, t_edge [E_CE], t_edge [E_DF], geom);
336 assoCircle (center, t_edge [E_CF], t_edge [E_DE], geom);
338 t_edge[E_AC]->getVertex(V_AMONT)->getPoint (sommet);
339 double radius = calc_distance (center, sommet);
341 int sphid = geom -> addSphere (center, radius);
344 for (int nf=0 ; nf < HQ_MAXI ; nf++)
345 t_quad[nf]->addAssociation (geom, sphid);
347 // ==================================================== propagateAssociation
348 int Elements::propagateAssociation (Edge* orig, Edge* dest, Edge* dir)
352 if (revo_lution || orig==NULL || dest==NULL || dir==NULL)
355 const Shapes& tab_shapes = orig->getAssociations ();
356 const Shapes& tab_dest = dest->getAssociations ();
357 int nbdest = tab_dest.size();
358 int nbshapes = tab_shapes.size();
359 bool on_edge = nbshapes!=0 && nbdest==0;
361 Vertex* vo1 = orig->commonVertex (dir);
362 Vertex* vd1 = dest->commonVertex (dir);
363 Vertex* vo2 = orig->opposedVertex (vo1);
364 Vertex* vd2 = dest->opposedVertex (vd1);
366 if (vo1==NULL || vd1==NULL)
370 Real3 pa, pb, vdir1, vdir2;
371 calc_vecteur (vo1->getPoint (pa), vd1->getPoint (pb), vdir1);
372 calc_vecteur (vo2->getPoint (pa), vd2->getPoint (pb), vdir2);
374 double dd = calc_distance (vdir1, vdir2);
375 bool para = dd < 1.0e-3;
379 for (int nro=0 ; nro<nbshapes ; nro++)
381 Shape* shape = tab_shapes[nro];
384 string brep = shape->getBrep();
385 translate_brep (brep, vdir1, trep);
386 // Shape* tshape = new Shape (trep);
387 // tshape->setBounds (shape->getStart(), shape->getEnd());
388 // dest->addAssociation (tshape);
393 double* vdir = vdir1;
394 for (int nro=V_AMONT ; nro<=V_AVAL ; nro++)
396 Shape* shape = vo1->getAssociation ();
397 if (shape!=NULL && vd1->getAssociation ()==NULL)
399 string brep = shape->getBrep();
400 translate_brep (brep, vdir, trep);
401 // Shape* tshape = new Shape (trep);
402 // vd1->setAssociation (tshape);
411 // ==================================================== prismAssociation
412 int Elements::prismAssociation (Edge* lorig, Edge* ldest, int nh)
414 if (lorig==NULL || ldest==NULL)
418 int nbshapes = lorig->countAssociation ();
422 NewShape* new_shape = getShape ();
423 for (int nro=0 ; nro<nbshapes ; nro++)
425 AssoEdge* asso = lorig->getAssociation (nro);
426 EdgeShape* shape = asso ->getEdgeShape();
427 double p1 = asso ->getStart();
428 double p2 = asso ->getEnd ();
432 sprintf (name, "0x%lx#%d", (unsigned long) shape, nh);
433 map<string,int>::iterator iter = map_shape.find (name);
434 if (iter != map_shape.end())
435 subid = iter->second;
437 subid = map_shape[name] = new_shape->transfoShape (cum_matrix, shape);
439 new_shape->addAssociation (ldest, subid, p1, p2);
440 printf (" prismAsso : %s -> %s(%d) = %s [ %g , %g]\n",
441 lorig->getName(), ldest->getName(), nh, name, p1, p2);
445 // ====================================================== assoResiduelle
446 void Elements::assoResiduelle ()
448 // int nbre = tab_vertex.size();
449 // for (int nv=0 ; nv<nbre ; nv++)
451 // geom_asso_point (tab_vertex [nv]);
454 // ====================================================== moveDisco
455 void Elements::moveDisco (Hexa* hexa)
459 hexa->getCenter (center);
460 matrix.defScale (center, 0.55);
462 int nbre = tab_vertex.size();
463 for (int nv=0 ; nv<nbre ; nv++)
465 matrix.perform (tab_vertex [nv]);
468 // =================================================== getStrate
469 Hexa* Elements::getStrate (int couche, EnumHQuad nroface)
472 int nro = couche <= 0 ? 0 : (couche-1)*HQ_MAXI + nroface + 1;
474 if (nbr_hexas==0 || nro >= nbr_hexas)
477 cell = tab_hexa [nro];
481 // ============================================================ setHexa
482 void Elements::setHexa (Hexa* elt, int nro)
484 if (nro >=0 && nro < nbr_hexas)
485 tab_hexa [nro] = elt;
487 // ============================================================ setQuad
488 void Elements::setQuad (Quad* elt, int nro)
490 if (nro >=0 && nro < nbr_quads)
491 tab_quad [nro] = elt;
493 // ============================================================ setEdge
494 void Elements::setEdge (Edge* elt, int nro)
496 if (nro >=0 && nro < nbr_edges)
497 tab_edge [nro] = elt;
499 // ============================================================ setVertex
500 void Elements::setVertex (Vertex* elt, int nro)
502 if (nro >=0 && nro < nbr_vertex)
503 tab_vertex [nro] = elt;
505 // -----------------------------------------------------------------------
506 // -----------------------------------------------------------------------
507 // ============================================================ getHexa
508 Hexa* Elements::getHexa (int nro)
510 DumpStart ("getHexa", nro);
513 int nombre=tab_hexa.size();
514 // if (nro >=0 && nro < nbr_hexas && el_status == HOK Abu 2010/05/06
515 if (nro >=0 && nro < nombre && el_status == HOK
516 && tab_hexa [nro] != NULL && tab_hexa [nro]->isValid())
517 elt = tab_hexa [nro];
522 // ============================================================ getQuad
523 Quad* Elements::getQuad (int nro)
525 DumpStart ("getQuad", nro);
528 if (nro >=0 && nro < nbr_quads && el_status == HOK
529 && tab_quad [nro] != NULL && tab_quad [nro]->isValid())
530 elt = tab_quad [nro];
535 // ============================================================ getEdge
536 Edge* Elements::getEdge (int nro)
538 DumpStart ("getEdge", nro);
541 if (nro >=0 && nro < nbr_edges && el_status == HOK
542 && tab_edge [nro] != NULL && tab_edge [nro]->isValid())
543 elt = tab_edge [nro];
548 // ============================================================ getVertex
549 Vertex* Elements::getVertex (int nro)
551 DumpStart ("getVertex", nro);
554 if (nro >=0 && nro < nbr_vertex && el_status == HOK
555 && tab_vertex [nro] != NULL && tab_vertex [nro]->isValid())
556 elt = tab_vertex [nro];
561 // ============================================================ indVertex
562 int Elements::indVertex (int nquad, int nsommet, int nh)
564 int nro = nsommet + QUAD4*nquad + nbr_orig*QUAD4*(nh+1);
567 // ============================================================ nroVertex
568 int Elements::nroVertex (int nquad, int nsommet, int nh)
570 int nro = nsommet + QUAD4*nquad + nbr_orig*QUAD4*nh;
573 // ============================================================ indVertex
574 int Elements::indVertex (int ref_edge, int nh)
576 int nro = ref_edge + nbr_orig*QUAD4*nh;
579 // ============================================================ nroEdgeH
580 int Elements::nroEdgeH (int nvertex)
582 return QUAD4*nbr_orig*gr_hauteur + nvertex;
584 // ============================================================ nroEdgeH
585 int Elements::nroEdgeH (int nquad, int nsommet, int nh)
587 return QUAD4*nbr_orig*gr_hauteur + indVertex (nquad, nsommet, nh);
589 // ============================================================ nroHexa
590 int Elements::nroHexa (int nquad, int nh)
592 int nro = gr_hauteur*nquad + nh;
596 // ============================================================ addHexa
597 void Elements::addHexa (Hexa* element)
599 tab_hexa.push_back (element);
602 // ============================================================ addQuad
603 void Elements::addQuad (Quad* element)
605 tab_quad.push_back (element);
608 // ============================================================ addEdge
609 void Elements::addEdge (Edge* element)
611 tab_edge.push_back (element);
614 // ============================================================ addVertex
615 void Elements::addVertex (Vertex* element)
617 tab_vertex.push_back (element);
620 // ============================================================ findHexa
621 int Elements::findHexa (Hexa* element)
623 int nbre = tab_hexa.size();
624 for (int nro=0 ; nro<nbre ; nro++)
625 if (tab_hexa [nro] == element)
629 // ============================================================ findQuad
630 int Elements::findQuad (Quad* element)
632 int nbre = tab_quad.size();
633 for (int nro=0 ; nro<nbre ; nro++)
634 if (tab_quad [nro] == element)
638 // ============================================================ findEdge
639 int Elements::findEdge (Edge* element)
641 int nbre = tab_edge.size();
642 for (int nro=0 ; nro<nbre ; nro++)
643 if (tab_edge [nro] == element)
647 // ============================================================ findVertex
648 int Elements::findVertex (Vertex* element)
650 int nbre = tab_vertex.size();
651 for (int nro=0 ; nro<nbre ; nro++)
652 if (tab_vertex [nro] == element)
656 // ========================================================= saveVtk (avec nro)
657 int Elements::saveVtk (cpchar radical, int &nro)
660 sprintf (num, "%d", nro);
663 string filename = radical;
666 int ier = saveVtk (filename.c_str());