2 // C++ : Grilles hexa 6
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"
24 #include "HexGlobale.hxx"
25 #include "HexDocument.hxx"
26 #include "HexVector.hxx"
27 #include "HexVertex.hxx"
28 #include "HexHexa.hxx"
29 #include "HexEdge.hxx"
30 #include "HexPropagation.hxx"
38 double calcul_phimax (double radhole, double radext, bool sphere);
40 // ====================================================== makeCartesian
41 int Elements::makeCartesian (Vertex* orig, double* ux, double* uy, double* uz,
42 RealVector& tx, RealVector& ty, RealVector& tz)
44 resize (GR_CARTESIAN, tx.size(), ty.size(), tz.size());
48 double px0 = orig->getX ();
49 double py0 = orig->getY ();
50 double pz0 = orig->getZ ();
52 setVertex (orig, 0, 0, 0);
53 for (int nz=0 ; nz<size_vz ; nz++)
55 kz = nz==0 ? 0 : tz [nz-1];
56 for (int ny=0 ; ny<size_vy ; ny++)
58 ky = ny==0 ? 0 : ty [ny-1];
59 for (int nx=0 ; nx<size_vx ; nx++)
61 kx = nx==0 ? 0 : tx [nx-1];
62 dx = px0 + kx*ux[dir_x] + ky*uy[dir_x] + kz*uz[dir_x];
63 dy = py0 + kx*ux[dir_y] + ky*uy[dir_y] + kz*uz[dir_y];
64 dz = pz0 + kx*ux[dir_z] + ky*uy[dir_z] + kz*uz[dir_z];
68 node = el_root->addVertex (dx, dy, dz);
69 setVertex (node, nx, ny, nz);
76 // ====================================================== makeCylinder
77 int Elements::makeCylinder (Vertex* vorig, double* base, double* haut,
78 RealVector& tray, RealVector& tang, RealVector& thaut,
81 Real3 orig = { 0, 0, 0 };
83 vorig->getPoint (orig);
85 int nr = tray .size() - 1;
86 int na = tang .size() - 1;
87 int nl = thaut.size() - 1;
88 double angle = tang[na];
90 resize (GR_CYLINDRIC, nr, na, nl);
91 cyl_closed = angle >= 359.9;
94 int ier = makeBasicCylinder (tray, tang, thaut, fill);
98 transfoVertices (orig, base, haut);
101 assoCylinders (orig, haut, angle, tang);
104 // ====================================================== transfoVertices
105 void Elements::transfoVertices (double* omega, double* base, double* haut)
110 prod_vectoriel (haut, base, jprim);
111 prod_vectoriel (jprim, haut, iprim);
113 matrice.define (omega, iprim, jprim, haut);
114 int nbre = tab_vertex.size ();
115 for (int nro=0 ; nro<nbre ; nro++)
117 if (tab_vertex[nro] != NULL)
118 tab_vertex[nro]->setMark (NO_USED);
121 for (int nro=0 ; nro<nbre ; nro++)
123 Vertex* node = tab_vertex[nro];
124 if (node != NULL && node->getMark() == NO_USED)
126 matrice.perform (node);
127 node->setMark (IS_USED);
131 // ====================================================== makeSpherical
132 int Elements::makeSpherical (double* center, double* base, double* haut,
133 RealVector& trayon, int critere)
136 prod_vectoriel (haut, base, jprim);
137 prod_vectoriel (jprim, haut, iprim);
139 int nbre = trayon.size();
140 resize (GR_SPHERIC, nbre);
142 Vertex* i_node [HV_MAXI]; // Les noeuds de l'hexa englobant
143 Edge* i_edge [HE_MAXI]; // Les noeuds de l'hexa englobant
144 Quad* i_quad [HQ_MAXI]; // Les noeuds de l'hexa englobant
146 double rayon = trayon [0];
147 for (int nro=0 ; nro<HV_MAXI; nro++)
149 double dx = glob->CoordVertex (nro, dir_x) * rayon;
150 double dy = glob->CoordVertex (nro, dir_y) * rayon;
151 double dz = glob->CoordVertex (nro, dir_z) * rayon;
153 double px = center[0] + dx*iprim[0] + dy*jprim[0] + dz*haut[0] ;
154 double py = center[1] + dx*iprim[1] + dy*jprim[1] + dz*haut[1] ;
155 double pz = center[2] + dx*iprim[2] + dy*jprim[2] + dz*haut[2] ;
157 i_node [nro] = el_root->addVertex (px, py, pz);
160 for (int nro=0 ; nro<HE_MAXI; nro++)
162 int v1 = glob->EdgeVertex (nro, V_AMONT);
163 int v2 = glob->EdgeVertex (nro, V_AVAL);
164 i_edge[nro] = newEdge (i_node[v1], i_node[v2]);
168 char nm0[8], nm1 [8], nm2 [8];
169 printf (" %2d : %s = %s = [%s, %s] = [%d,%d] = [%s,%s]\n", nro,
170 glob->namofHexaEdge(nro), i_edge[nro]->getName(nm0),
171 glob->namofHexaVertex(v1), glob->namofHexaVertex(v2), v1, v2,
172 i_node[v1]->getName(nm1), i_node[v2]->getName(nm2));
176 for (int nro=0 ; nro<HQ_MAXI; nro++)
177 i_quad[nro] = newQuad (i_edge[glob->QuadEdge (nro, E_A)],
178 i_edge[glob->QuadEdge (nro, E_B)],
179 i_edge[glob->QuadEdge (nro, E_C)],
180 i_edge[glob->QuadEdge (nro, E_D)]);
182 tab_hexa.push_back (newHexa (i_quad[Q_A], i_quad[Q_B], i_quad[Q_C],
183 i_quad[Q_D], i_quad[Q_E], i_quad[Q_F]));
185 for (int niv=1; niv<gr_rayon ; niv++)
187 double rayon0 = rayon;
188 rayon = trayon [niv];
189 addStrate (i_quad, i_edge, i_node, center, rayon/rayon0);
192 nbr_hexas = tab_hexa.size();
193 assoSphere (center, i_edge, i_quad);
196 // ====================================================== addStrate
197 int Elements::addStrate (Quad* i_quad[], Edge* i_edge[], Vertex* i_node[],
198 Real* center, double lambda)
200 Vertex* e_node [HV_MAXI]; // Les noeuds de l'hexa englobant
201 Edge* e_edge [HE_MAXI]; // Les noeuds de l'hexa englobant
202 Quad* e_quad [HQ_MAXI]; // Les noeuds de l'hexa englobant
204 Edge* d_edge [HV_MAXI]; // Les aretes diagonales (1 par sommet)
205 Quad* d_quad [HE_MAXI]; // Les faces diagonales (1 par arete)
208 // + les aretes diagonales
209 double px0 = center [dir_x];
210 double py0 = center [dir_y];
211 double pz0 = center [dir_z];
213 for (int nv=0 ; nv<HV_MAXI ; nv++)
215 e_node[nv] = el_root->addVertex (px0+lambda*(i_node[nv]->getX()-px0),
216 py0+lambda*(i_node[nv]->getY()-py0),
217 pz0+lambda*(i_node[nv]->getZ()-pz0));
219 d_edge[nv] = newEdge (i_node[nv], e_node[nv]);
221 // Les aretes exterieures
222 // + les faces diagonales
223 for (int nro=0 ; nro<HE_MAXI ; nro++)
225 int nv0 = glob->EdgeVertex (nro, V_AMONT);
226 int nv1 = glob->EdgeVertex (nro, V_AVAL );
227 e_edge[nro] = newEdge (e_node [nv0], e_node [nv1]);
228 d_quad[nro] = newQuad (i_edge [nro], d_edge [nv0],
229 e_edge [nro], d_edge [nv1]);
231 // Les faces exterieures
234 for (int nro=0 ; nro<HQ_MAXI ; nro++)
236 int ne0 = glob->QuadEdge (nro, E_A);
237 int ne1 = glob->QuadEdge (nro, E_B);
238 int ne2 = glob->QuadEdge (nro, E_C);
239 int ne3 = glob->QuadEdge (nro, E_D);
241 e_quad[nro] = newQuad (e_edge[ne0], e_edge[ne1],
242 e_edge[ne2], e_edge[ne3]);
243 strate = newHexa (i_quad[nro], e_quad[nro], d_quad[ne0],
244 d_quad[ne2], d_quad[ne1], d_quad[ne3]);
245 tab_hexa.push_back (strate);
248 for (int nv=0 ; nv<HV_MAXI ; nv++) i_node [nv] = e_node [nv];
249 for (int ns=0 ; ns<HE_MAXI ; ns++) i_edge [ns] = e_edge [ns];
250 for (int nq=0 ; nq<HQ_MAXI ; nq++) i_quad [nq] = e_quad [nq];
253 // ====================================================== makeRind
254 int Elements::makeRind (EnumGrid type, double* center, double* vx, double* vz,
255 RealVector& tray, RealVector& tang, RealVector& tphi)
257 int na = tang.size() - 1;
258 int nr = tray.size() - 1;
259 int nh = tphi.size() - 1;
264 if (type==GR_PART_RIND || type==GR_RIND)
272 resize (type, nr, na, nh);
274 cyl_radhole = tray[0];
275 cyl_radint = tray[nint];
276 cyl_radext = tray[rad0+nr];
277 cyl_closed = type==GR_HEMISPHERIC || type==GR_RIND;
278 cyl_phimax = calcul_phimax (cyl_radhole, cyl_radext, sphere);
279 PutData (cyl_phimax);
280 int nb_secteurs = cyl_closed ? size_vy-1 : size_vy;
282 for (int ny=0 ; ny<nb_secteurs ; ny++)
284 double theta = tang[ny];
285 for (int nx=0 ; nx<size_vx ; nx++)
287 double rayon = tray[rad0+nx];
288 if (sphere && ny==0 && nx==0)
289 printf (" rayon = %g, rad_int=%g\n", rayon, cyl_radint);
291 for (int nz=0 ; nz<size_vz ; nz++)
294 double phi = tphi [nz];
295 makeRindPoint (rayon, theta, phi, px, py, pz);
296 Vertex* node = el_root->addVertex (px, py, pz);
297 setVertex (node, nx, ny, nz);
304 for (int nx=0 ; nx<size_vx ; nx++)
305 for (int nz=0 ; nz<size_vz ; nz++)
307 Vertex* node = getVertexIJK (nx, 0, nz);
308 setVertex (node, nx, size_vy-1, nz);
312 transfoVertices (center, vx, vz);
314 double angle = tang[na];
315 assoCylinders (center, vz, angle, tang);
318 // ====================================================== makeRindPoint
319 int Elements::makeRindPoint (double ray, double theta, double phi,
320 double& px, double& py, double& pz)
322 bool rind = (grid_type == GR_RIND || grid_type == GR_PART_RIND);
323 double rphi = deg2radians (phi);
324 rphi = std::min (cyl_phimax, std::max(-cyl_phimax, rphi));
325 double rtheta = deg2radians (theta);
326 double sinphi = sin (rphi);
327 double cosphi = cos (rphi);
332 // rayon = cyl_radint + nr*(cyl_radext-cyl_radint)/size_hx;
333 // pz = rayon*sinphi;
334 // rayon = rayon*cosphi;
340 // rayon = cyl_radhole + nr*(cyl_radext*cosphi - cyl_radhole)/size_hx;
341 pz = cyl_radext*sinphi;
342 rayon = cyl_radhole + (ray-cyl_radhole)*(cyl_radext*cosphi-cyl_radhole)
343 / (cyl_radext-cyl_radhole);
346 rayon = std::max (cyl_radhole, rayon);
347 px = rayon * cos (rtheta);
348 py = rayon * sin (rtheta);
351 // ====================================================== extrudeQuads
352 int Elements::extrudeQuads (Quads& tstart, double* axe, RealVector& tabval,
353 bool revol, Vertex* center)
356 int nbiter = tabval.size ();
357 int nbcells = tstart.size ();
359 el_root->markAll (NO_USED);
363 nbr_hexas = nbcells*nbiter;
365 tab_hexa.resize (nbr_hexas, NULL);
366 tab_quad.clear (); // verticaux
367 ker_hquad.clear (); // Horizontaux
371 copy_vecteur (axe, prism_dir);
374 prism_vec = NOT revol;
375 copy_vecteur (axe, revo_axe);
376 revo_center = center;
378 cum_values = gen_values;
380 for (int nro=1 ; nro<nbiter ; nro++)
381 gen_values [nro] = cum_values [nro] - cum_values [nro-1];
386 for (int nro=0 ; nro<nbcells ; nro++)
388 prismHexas (nro, tstart[nro], nbiter);
394 // ============================================================ joinQuads
395 int Elements::joinQuads (Quads& tstart, Quad* cible, Vertex* va1, Vertex* vb1,
396 Vertex* va2, Vertex* vb2, RealVector& tlen)
399 int nbiter = tlen .size();
400 int nbquads = tstart.size();
401 resize (GR_JOINT, nbquads, nbiter);
403 el_root->markAll (IS_NONE);
411 for (int nq=0 ; nq<nbquads ; nq++)
412 tab_orig [nq]->setMark (nq);
414 StrOrient orient (va1, va2, vb1, vb2);
415 int ier = this ->coupler (0, cible, &orient);
421 ier = tstart[0]->coupler (cible, &orient, this);
424 // ============================================================ cutEdge
425 int Elements::cutEdge (Edge* edge, RealVector& tlen)
429 cum_values = gen_values;
431 Propagation* prop = el_root->findPropagation (edge);
432 const Edges& t_edges = prop->getEdges ();
434 int nbcuts = tlen.size ();
435 cutHexas (t_edges, nbcuts);
437 // el_root->majPropagation ();
440 // ============================================================ nearestVertex
441 Vertex* Elements::nearestVertex (Vertex* other)
443 if (BadElement (other))
448 int nbre = countVertex ();
449 for (int nro=0 ; nro<nbre ; nro++)
451 Vertex* vert = getVertex (nro);
452 double dist = other->dist2 (vert);