TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Hexaedre_VEF.cpp
1/****************************************************************************
2* Copyright (c) 2024, CEA
3* All rights reserved.
4*
5* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
6* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
7* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
8* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
9*
10* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
11* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
12* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
13*
14*****************************************************************************/
15
16#include <Hexaedre_VEF.h>
17#include <Domaine.h>
18
19Implemente_instanciable_32_64(Hexaedre_VEF_32_64,"Hexaedre_VEF",Elem_geom_base_32_64<_T_>);
20
21
22/*! @brief NE FAIT RIEN
23 *
24 * @param (Sortie& s) un flot de sortie
25 * @return (Sortie&) le flot de sortie
26 */
27template <typename _SIZE_>
28Sortie& Hexaedre_VEF_32_64<_SIZE_>::printOn(Sortie& s ) const
29{
30 return s;
31}
32
33
34/*! @brief NE FAIT RIEN
35 *
36 * @param (Entree& s) un flot d'entree
37 * @return (Entree&) le flot d'entree
38 */
39template <typename _SIZE_>
40Entree& Hexaedre_VEF_32_64<_SIZE_>::readOn(Entree& s )
41{
42 return s;
43}
44
45/*! @brief Renvoie le nom LML d'un Hexaedre_VEF = "HEXA8".
46 *
47 * @return (Nom&) toujours egal a "HEXA8"
48 */
49template <typename _SIZE_>
50const Nom& Hexaedre_VEF_32_64<_SIZE_>::nom_lml() const
51{
52 static Nom nom="VOXEL8";
53 return nom;
54}
55
56/*! @brief teste si une position pos est un tetraedre de sommets som(i)
57 */
58template <typename _SIZE_>
59bool Hexaedre_VEF_32_64<_SIZE_>::entre_faces(const ArrOfDouble& pos, int_t Asom0_, int_t Asom1_, int_t Asom2_,
60 int_t Bsom0_, int_t Bsom1_, int_t Bsom2_) const
61{
62 const Domaine_t& dom = this->mon_dom.valeur();
63 double prodA,prodB;
64 const DoubleTab_t& coord=dom.les_sommets();
65
66 // pour A
67 double vecA1[3], vecA2[3], AP[3];
68 double normaleA[3];
69 for (int i=0; i<3; i++)
70 {
71 double ref=coord(Asom0_,i);
72 vecA1[i]=coord(Asom1_,i)-ref;
73 vecA2[i]=coord(Asom2_,i)-ref;
74 AP[i]=pos[i]-ref;
75 }
76 normaleA[0]=(vecA1[1]*vecA2[2]-vecA1[2]*vecA2[1]);
77 normaleA[1]=(vecA1[2]*vecA2[0]-vecA1[0]*vecA2[2]);
78 normaleA[2]=(vecA1[0]*vecA2[1]-vecA1[1]*vecA2[0]);
79 prodA=0;
80 for (int i=0; i<3; i++)
81 prodA+=normaleA[i]*AP[i];
82
83 // pour B
84 double vecB1[3], vecB2[3], BP[3];
85 double normaleB[3];
86 for (int i=0; i<3; i++)
87 {
88 double ref=coord(Bsom0_,i);
89 vecB1[i]=coord(Bsom1_,i)-ref;
90 vecB2[i]=coord(Bsom2_,i)-ref;
91 BP[i]=pos[i]-ref;
92 }
93 normaleB[0]=(vecB1[1]*vecB2[2]-vecB1[2]*vecB2[1]);
94 normaleB[1]=(vecB1[2]*vecB2[0]-vecB1[0]*vecB2[2]);
95 normaleB[2]=(vecB1[0]*vecB2[1]-vecB1[1]*vecB2[0]);
96 prodB=0;
97 for (int i=0; i<3; i++)
98 prodB+=normaleB[i]*BP[i];
99
100 if (prodA*prodB<=0)
101 return true;
102 else
103 {
104 // on regarde si on n'est pas limite
105 double rap=prodA/prodB;
106 if (rap>1) rap=prodB/prodA;
107 if (rap< Objet_U::precision_geom) return true;
108 }
109 return false;
110}
111
112/*! @brief teste si une position pos est un tetraedre de sommets som(i)
113 *
114 * @return (Nom&) toujours egal a "HEXA8"
115 */
116template <typename _SIZE_>
117bool Hexaedre_VEF_32_64<_SIZE_>::contient_Tetra(const ArrOfDouble& pos, int_t som0_, int_t som1_, int_t som2_, int_t som3_, int aff) const
118{
119 const Domaine_t& dom = this->mon_dom.valeur();
120 double prod1,prod2;
121 int_t som0, som1, som2, som3;
122 som0 = som0_;
123 som1 = som1_;
124 som2 = som2_;
125 som3 = som3_;
126 double vec0[3], vec1[3], vec2[3];
127 double vecx0[3], vecx1[3], vecx2[3];
128
129 if (aff)
130 {
131 Cerr<<"Test contient_Tetra nodes="<<som0_<<" "<<som1_<<" "<<som2_<<" "<<som3_<<finl;
132 Cerr<<" position="<<pos[0]<<" "<<pos[1]<<" "<<pos[2]<<finl;
133 }
134
135 if( ( est_egal(dom.coord(som0,0),pos[0]) && est_egal(dom.coord(som0,1),pos[1]) && est_egal(dom.coord(som0,2),pos[2]) )
136 || (est_egal(dom.coord(som1,0),pos[0]) && est_egal(dom.coord(som1,1),pos[1]) && est_egal(dom.coord(som1,2),pos[2]))
137 || (est_egal(dom.coord(som2,0),pos[0]) && est_egal(dom.coord(som2,1),pos[1]) && est_egal(dom.coord(som2,2),pos[2]))
138 || (est_egal(dom.coord(som3,0),pos[0]) && est_egal(dom.coord(som3,1),pos[1]) && est_egal(dom.coord(som3,2),pos[2])) )
139 return 1;
140
141 for (int j=0; j<4; j++)
142 {
143 switch(j)
144 {
145 case 0 :
146 som0 = som0_;
147 som1 = som1_;
148 som2 = som2_;
149 som3 = som3_;
150 break;
151 case 1 :
152 som0 = som1_;
153 som1 = som2_;
154 som2 = som3_;
155 som3 = som0_;
156 break;
157 case 2 :
158 som0 = som2_;
159 som1 = som3_;
160 som2 = som0_;
161 som3 = som1_;
162 break;
163 case 3 :
164 som0 = som3_;
165 som1 = som0_;
166 som2 = som1_;
167 som3 = som2_;
168 break;
169 }
170
171 //algo : les produits mixtes (som0, som1, som2,som3) et (pos, som1, som2,som3) doivent a voir le meme signe
172 //construction des vecteurs
173 vec0[0] = dom.coord(som1,0) - dom.coord(som0,0);
174 vec0[1] = dom.coord(som1,1) - dom.coord(som0,1);
175 vec0[2] = dom.coord(som1,2) - dom.coord(som0,2);
176 vec1[0] = dom.coord(som2,0) - dom.coord(som0,0);
177 vec1[1] = dom.coord(som2,1) - dom.coord(som0,1);
178 vec1[2] = dom.coord(som2,2) - dom.coord(som0,2);
179 vec2[0] = dom.coord(som3,0) - dom.coord(som0,0);
180 vec2[1] = dom.coord(som3,1) - dom.coord(som0,1);
181 vec2[2] = dom.coord(som3,2) - dom.coord(som0,2);
182 vecx0[0] = dom.coord(som1,0) - pos[0];
183 vecx0[1] = dom.coord(som1,1) - pos[1];
184 vecx0[2] = dom.coord(som1,2) - pos[2];
185 vecx1[0] = dom.coord(som2,0) - pos[0];
186 vecx1[1] = dom.coord(som2,1) - pos[1];
187 vecx1[2] = dom.coord(som2,2) - pos[2];
188 vecx2[0] = dom.coord(som3,0) - pos[0];
189 vecx2[1] = dom.coord(som3,1) - pos[1];
190 vecx2[2] = dom.coord(som3,2) - pos[2];
191 prod1 = vec0[0]*vec1[1]*vec2[2] + vec0[1]*vec1[2]*vec2[0] + vec0[2]*vec1[0]*vec2[1]
192 - vec0[0]*vec1[2]*vec2[1] - vec0[1]*vec1[0]*vec2[2] - vec0[2]*vec1[1]*vec2[0];
193 prod2 = vecx0[0]*vecx1[1]*vecx2[2] + vecx0[1]*vecx1[2]*vecx2[0] + vecx0[2]*vecx1[0]*vecx2[1]
194 - vecx0[0]*vecx1[2]*vecx2[1] - vecx0[1]*vecx1[0]*vecx2[2] - vecx0[2]*vecx1[1]*vecx2[0];
195
196 if (aff)
197 {
198 Cerr<<" test "<<som0<<" "<<som1<<" "<<som2<<" "<<som3<<finl;
199 Cerr<<" node0="<<som0<<" : "<<dom.coord(som0,0)<<" "<<dom.coord(som0,1)<<" "<<dom.coord(som0,2)<<finl;
200 Cerr<<" node1="<<som1<<" : "<<dom.coord(som1,0)<<" "<<dom.coord(som1,1)<<" "<<dom.coord(som1,2)<<finl;
201 Cerr<<" node2="<<som2<<" : "<<dom.coord(som2,0)<<" "<<dom.coord(som2,1)<<" "<<dom.coord(som2,2)<<finl;
202 Cerr<<" node3="<<som3<<" : "<<dom.coord(som3,0)<<" "<<dom.coord(som3,1)<<" "<<dom.coord(som3,2)<<finl;
203 Cerr<<" position : "<<pos[0]<<" "<<pos[1]<<" "<<pos[2]<<finl;
204 Cerr<<" prod1="<<prod1<<" prod2="<<prod2<<finl;
205 }
206
207 if (std::fabs(prod1)<Objet_U::precision_geom)
208 {
209 double norm_v0=0,norm_v1=0,norm_v2=0;
210 for (int ii=0; ii<3; ii++)
211 {
212 norm_v0+=vec0[ii]*vec0[ii];
213 norm_v1+=vec1[ii]*vec1[ii];
214 norm_v2+=vec2[ii]*vec2[ii];
215 }
216 if (std::fabs(prod1)<Objet_U::precision_geom*sqrt(norm_v0*norm_v1*norm_v2))
217 {
218 Cerr<<"Test contient_Tetra ERROR : coplanar nodes"<<finl;
219 Cerr<<" node0="<<som0<<" : "<<dom.coord(som0,0)<<" "<<dom.coord(som0,1)<<" "<<dom.coord(som0,2)<<finl;
220 Cerr<<" node1="<<som1<<" : "<<dom.coord(som1,0)<<" "<<dom.coord(som1,1)<<" "<<dom.coord(som1,2)<<finl;
221 Cerr<<" node2="<<som2<<" : "<<dom.coord(som2,0)<<" "<<dom.coord(som2,1)<<" "<<dom.coord(som2,2)<<finl;
222 Cerr<<" node3="<<som3<<" : "<<dom.coord(som3,0)<<" "<<dom.coord(som3,1)<<" "<<dom.coord(som3,2)<<finl;
223 Cerr<<" position : "<<pos[0]<<" "<<pos[1]<<" "<<pos[2]<<finl;
224 Process::Process::exit();
225 }
226 }
227
228 if (prod1*prod2<0 && std::fabs(prod2)>std::fabs(prod1)*Objet_U::precision_geom)
229 {
230 if (aff)
231 {
232 Cerr<<" CONTAINS=0"<<finl;
233 }
234 return 0;
235 }
236 }
237 if (aff)
238 {
239 Cerr<<" CONTAINS=1"<<finl;
240 }
241 return 1;
242}
243
244/*! @brief Renvoie 1 si l'element ielem du domaine associe a l'element geometrique contient le point
245 *
246 * de coordonnees specifiees par le parametre "pos".
247 * Renvoie 0 sinon.
248 *
249 * @param (DoubleVect& pos) coordonnees du point que l'on cherche a localiser
250 * @param (int ielem) le numero de l'element du domaine dans lequel on cherche le point.
251 * @return (int) 1 si le point de coordonnees specifiees appartient a l'element ielem 0 sinon
252 */
253template <typename _SIZE_>
254int Hexaedre_VEF_32_64<_SIZE_>::contient(const ArrOfDouble& pos, int_t element ) const
255{
256 assert(pos.size_array()==3);
257 const Domaine_t& dom=this->mon_dom.valeur();
258 int_t som0 = dom.sommet_elem(element,0);
259 int_t som1 = dom.sommet_elem(element,1);
260 int_t som2 = dom.sommet_elem(element,2);
261 int_t som3 = dom.sommet_elem(element,3);
262 int_t som4 = dom.sommet_elem(element,4);
263 int_t som5 = dom.sommet_elem(element,5);
264 int_t som6 = dom.sommet_elem(element,6);
265 int_t som7 = dom.sommet_elem(element,7);
266
267 // GF on teste si le point est boen entre les 2 faces
268 bool new_algo = true;
269 if (new_algo)
270 {
271 if (entre_faces(pos, som0, som1, som2, som4, som5, som6))
272 if (entre_faces(pos, som0, som1, som4, som2, som3, som6))
273 if (entre_faces(pos, som0, som2, som4, som1, som3, som5))
274 return true;
275 return false;
276 }
277 else
278 {
279 // PQ : 04/03 : optimisation
280 int aff = 0;
281 if (element == -701)
282 aff = 1;
283 if (contient_Tetra(pos, som0, som1, som2, som4, aff) ||
284 contient_Tetra(pos, som1, som2, som3, som6, aff) ||
285 contient_Tetra(pos, som1, som2, som4, som6, aff) ||
286 contient_Tetra(pos, som3, som5, som6, som7, aff) ||
287 contient_Tetra(pos, som1, som4, som5, som6, aff) ||
288 contient_Tetra(pos, som1, som3, som5, som6, aff))
289 {
290 return true;
291 }
292 return false;
293 }
294}
295
296
297/*! @brief Renvoie 1 si les sommets specifies par le parametre "pos" sont les sommets de l'element "element" du domaine associe a
298 *
299 * l'element geometrique.
300 *
301 * @param (IntVect& pos) les numeros des sommets a comparer avec ceux de l'elements "element"
302 * @param (int element) le numero de l'element du domaine dont on veut comparer les sommets
303 * @return (int) 1 si les sommets passes en parametre sont ceux de l'element specifie, 0 sinon
304 */
305template <typename _SIZE_>
306int Hexaedre_VEF_32_64<_SIZE_>::contient(const SmallArrOfTID_t& som, int_t element ) const
307{
308 const Domaine_t& dom=this->mon_dom.valeur();
309 if((dom.sommet_elem(element,0)==som[0])&&
310 (dom.sommet_elem(element,1)==som[1])&&
311 (dom.sommet_elem(element,2)==som[2])&&
312 (dom.sommet_elem(element,3)==som[3])&&
313 (dom.sommet_elem(element,4)==som[4])&&
314 (dom.sommet_elem(element,5)==som[5])&&
315 (dom.sommet_elem(element,6)==som[6])&&
316 (dom.sommet_elem(element,7)==som[7]))
317 return 1;
318 else
319 return 0;
320}
321
322/*! @brief Calcule les volumes des elements du domaine associe.
323 *
324 * @param (DoubleVect& volumes) le vecteur contenant les valeurs des des volumes des elements du domaine
325 */
326template <typename _SIZE_>
327void Hexaedre_VEF_32_64<_SIZE_>::calculer_volumes(DoubleVect_t& volumes) const
328{
329 const Domaine_t& dom=this->mon_dom.valeur();
330 IntTab_t face_sommet_global;
331 face_sommet_global.resize(6,4);
332 double x0,y0,z0;
333 double x1,y1,z1;
334 double x2,y2,z2;
335 double x3,y3,z3;
336 double volume;
337 // double volume_global;
338 int_t s1,s2,s3,s4;
339 int_t som0,som1,som2,som3,som4,som5,som6,som7;
340
341 // volume_global =0;
342
343 int_t size_tot = dom.nb_elem_tot();
344 assert(volumes.size_totale()==size_tot);
345 for (int_t num_poly=0; num_poly<size_tot; num_poly++)
346 {
347 // on initialise le volume
348 // volumes[num_poly] =0;
349 // pour chaque elem on le decoupe en tetra.
350 // On calcul le volume par somme des volumes des tetras.
351
352 // on definit les sommets de l'elem
353 // on recupere le centre de l'element.
354
355 som0 = dom.sommet_elem(num_poly,0);
356 som1 = dom.sommet_elem(num_poly,1);
357 som2 = dom.sommet_elem(num_poly,2);
358 som3 = dom.sommet_elem(num_poly,3);
359 som4 = dom.sommet_elem(num_poly,4);
360 som5 = dom.sommet_elem(num_poly,5);
361 som6 = dom.sommet_elem(num_poly,6);
362 som7 = dom.sommet_elem(num_poly,7);
363
364 // Cerr << " som0 " << som0 << " som1 " << som1 << " som2 " << som2 << " som3 " << som3 << finl;
365 // Cerr << " som4 " << som4 << " som5 " << som5 << " som6 " << som6 << " som7 " << som7 << finl;
366
367
368 // on construit le tableau de travail face_sommet_global
369 // les faces seront ordonnee suivant Xmin, Xmax, Ymin, Ymax, Zmin,Zmax
370 // toujours en gardant la numerotation de TRUST
371 face_sommet_global(0,0) =som0;
372 face_sommet_global(0,1) =som1;
373 face_sommet_global(0,2) =som2;
374 face_sommet_global(0,3) =som3;
375
376 face_sommet_global(1,0) =som4;
377 face_sommet_global(1,1) =som5;
378 face_sommet_global(1,2) =som6;
379 face_sommet_global(1,3) =som7;
380
381 face_sommet_global(2,0) =som0;
382 face_sommet_global(2,1) =som4;
383 face_sommet_global(2,2) =som2;
384 face_sommet_global(2,3) =som6;
385
386 face_sommet_global(3,0) =som1;
387 face_sommet_global(3,1) =som5;
388 face_sommet_global(3,2) =som3;
389 face_sommet_global(3,3) =som7;
390
391 face_sommet_global(4,0) =som0;
392 face_sommet_global(4,1) =som4;
393 face_sommet_global(4,2) =som1;
394 face_sommet_global(4,3) =som5;
395
396 face_sommet_global(5,0) =som2;
397 face_sommet_global(5,1) =som6;
398 face_sommet_global(5,2) =som3;
399 face_sommet_global(5,3) =som7;
400
401 // le sommet 0 sera toujours le centre de gravite
402 x0 = ( dom.coord(som0,0) + dom.coord(som1,0) + dom.coord(som2,0) + dom.coord(som3,0) + dom.coord(som4,0) + dom.coord(som5,0) + dom.coord(som6,0) + dom.coord(som7,0) )*0.125;
403 y0 = ( dom.coord(som0,1) + dom.coord(som1,1) + dom.coord(som2,1) + dom.coord(som3,1) + dom.coord(som4,1) + dom.coord(som5,1) + dom.coord(som6,1) + dom.coord(som7,1) )*0.125;
404 z0 = ( dom.coord(som0,2) + dom.coord(som1,2) + dom.coord(som2,2) + dom.coord(som3,2) + dom.coord(som4,2) + dom.coord(som5,2) + dom.coord(som6,2) + dom.coord(som7,2) )*0.125;
405
406 // Cerr << "le num poly traite " << num_poly << finl;
407 // Cerr << "le centre de gravite " << x0 << " " << y0 << " " << z0 << finl;
408
409 volume =0;
410
411
412 // On decoupe en tetra.
413 // pour chaque face de l'hexa ( 0 a 7 )
414 for ( int k=0 ; k<nb_faces(); k++)
415 {
416 // on recupere les sommets de la face consideree
417 s1 = face_sommet_global(k,0);
418 s2 = face_sommet_global(k,1);
419 s3 = face_sommet_global(k,2);
420 s4 = face_sommet_global(k,3);
421
422 // Cerr << " les sommets de la face " << k << " = " << s1 << " " << s2 << " " << s3 << " " << s4 << finl;
423 // Cerr << "coord de s1 " << dom.coord(s1,0) << " " << dom.coord(s1,1) << " " << dom.coord(s1,2) << finl;
424 // Cerr << "coord de s2 " << dom.coord(s2,0) << " " << dom.coord(s2,1) << " " << dom.coord(s2,2) << finl;
425 // Cerr << "coord de s3 " << dom.coord(s3,0) << " " << dom.coord(s3,1) << " " << dom.coord(s3,2) << finl;
426 // Cerr << "coord de s4 " << dom.coord(s4,0) << " " << dom.coord(s4,1) << " " << dom.coord(s4,2) << finl;
427
428 x1 = dom.coord(s1,0);
429 y1 = dom.coord(s1,1);
430 z1 = dom.coord(s1,2);
431
432 x2 = dom.coord(s2,0);
433 y2 = dom.coord(s2,1);
434 z2 = dom.coord(s2,2);
435
436 x3 = dom.coord(s4,0);
437 y3 = dom.coord(s4,1);
438 z3 = dom.coord(s4,2);
439
440 //volume_tetra1 = std::fabs((x1-x0)*((y2-y0)*(z3-z0)-(y3-y0)*(z2-z0))-
441 // (x2-x0)*((y1-y0)*(z3-z0)-(y3-y0)*(z1-z0))+
442 // (x3-x0)*((y1-y0)*(z2-z0)-(y2-y0)*(z1-z0)))/6;
443
444 volume += std::fabs((x1-x0)*((y2-y0)*(z3-z0)-(y3-y0)*(z2-z0))-
445 (x2-x0)*((y1-y0)*(z3-z0)-(y3-y0)*(z1-z0))+
446 (x3-x0)*((y1-y0)*(z2-z0)-(y2-y0)*(z1-z0)))/6;
447
448 x1 = dom.coord(s1,0);
449 y1 = dom.coord(s1,1);
450 z1 = dom.coord(s1,2);
451
452 x2 = dom.coord(s4,0);
453 y2 = dom.coord(s4,1);
454 z2 = dom.coord(s4,2);
455
456 x3 = dom.coord(s3,0);
457 y3 = dom.coord(s3,1);
458 z3 = dom.coord(s3,2);
459
460 // volume_tetra2 = std::fabs((x1-x0)*((y2-y0)*(z3-z0)-(y3-y0)*(z2-z0))-
461 // (x2-x0)*((y1-y0)*(z3-z0)-(y3-y0)*(z1-z0))+
462 // (x3-x0)*((y1-y0)*(z2-z0)-(y2-y0)*(z1-z0)))/6;
463
464 volume += std::fabs((x1-x0)*((y2-y0)*(z3-z0)-(y3-y0)*(z2-z0))-
465 (x2-x0)*((y1-y0)*(z3-z0)-(y3-y0)*(z1-z0))+
466 (x3-x0)*((y1-y0)*(z2-z0)-(y2-y0)*(z1-z0)))/6;
467 }
468
469 volumes[num_poly] = volume;
470 }
471
472}
473
474template <typename _SIZE_>
475int Hexaedre_VEF_32_64<_SIZE_>::reordonne(int i0,int i1,int i2,int i3,const DoubleTab_t& coord,IntTab_t& elem, int_t& num_poly,int test,
476 DoubleTab& v, ArrOfDouble& prod_,ArrOfDouble& prod_v, ArrOfDouble& dist, DoubleTab& prod_v2)
477{
478 int_t s[4];
479 s[0]=elem(num_poly,i0);
480 s[1]=elem(num_poly,i1);
481 s[2]=elem(num_poly,i2);
482 s[3]=elem(num_poly,i3);
483 // DoubleTab v(3,3);
484 for (int i=1; i<4; i++)
485 for (int dir=0; dir<3; dir++)
486 v(i-1,dir)=coord(s[i],dir)-coord(s[0],dir);
487 //ArrOfDouble prod_(3);
488 int opp=-1;
489 // avant tout on teste si on a bien un plan
490 {
491
492 // ArrOfDouble prod_v(3);
493 int op=0,i=1;
494 // prod_v=O1^O2
495 prod_v[0]=v(op,1)*v(i,2)-v(op,2)*v(i,1);
496 prod_v[1]=v(op,2)*v(i,0)-v(op,0)*v(i,2);
497 prod_v[2]=v(op,0)*v(i,1)-v(op,1)*v(i,0);
498 // ArrOfDouble dist(8);
499 dist=0;
500 for (int j=0; j<8; j++)
501 {
502 for (int dir=0; dir<3; dir++)
503 dist[j]+=prod_v[dir]*(coord(elem(num_poly,j),dir)-coord(s[0],dir));
504 }
505 //Cerr<<"ici "<<i0 <<" i1 " <<i1<<" i2 "<<i2 << " i3 "<<i3 <<" "<<dist<<finl;
506 // sommet le plus proche du plan
507 int s_plan=0;
508 double min_dist=1e25;
509 for (int j=0; j<8; j++)
510 {
511 if ((j!=i0)&&(j!=i1)&&(j!=i2))
512 if (std::fabs(dist[j])<min_dist)
513 {
514 s_plan=j;
515 min_dist=std::fabs(dist[j]);
516 }
517 }
518 if (s_plan!=i3)
519 {
520 // pas un plan
521 if (test)
522 {
523 //Cerr<<"ici pb"<<finl;
524 return 1;
525 }
526 else
527 {
528 Cerr<<"we do not have a plan "<<finl;
529 Process::Process::exit();
530 }
531 }
532 // doit on tester si les 4 autres points sont du meme cote ???
533
534 }
535 for (int op=0; op<3; op++)
536 {
537 //DoubleTab prod_v2(3,3);
538 for (int i=0; i<3; i++)
539 {
540 prod_v2(i,0)=v(op,1)*v(i,2)-v(op,2)*v(i,1);
541 prod_v2(i,1)=v(op,2)*v(i,0)-v(op,0)*v(i,2);
542 prod_v2(i,2)=v(op,0)*v(i,1)-v(op,1)*v(i,0);
543 }
544 double prod=0;
545 int i1bis=-1,j2bis=-1;
546 if (op==0)
547 {
548 i1bis=1;
549 j2bis=2;
550 }
551 else if (op==1)
552 {
553 i1bis=0;
554 j2bis=2;
555 }
556 else if (op==2)
557 {
558 i1bis=0;
559 j2bis=1;
560 }
561 for (int dir=0; dir<3; dir++)
562 prod+=(prod_v2(i1bis,dir)*prod_v2(j2bis,dir));
563 prod_[op]=prod;
564 if (prod<0)
565 {
566 if (opp!=-1)
567 {
568 op=3;
569 if (test) return 1;
570 // on a deux produits negatifs on ne fait rien
571 }
572 opp=op;
573
574 }
575 }
576 if ((test==0)&& ((opp==3)||(opp==-1)))
577 {
578 Cerr<<" this does not seem to be a plan "<<finl;
579 Process::Process::exit();
580 }
581 if ((opp!=2)&&(opp!=-2))
582 {
583 int i2b=i2;
584 if (opp==0) i2b=i1;
585 if (test==0)
586 {
587 int_t tmp=elem(num_poly,i2b);
588 elem(num_poly,i2b)=elem(num_poly,i3);
589 elem(num_poly,i3)=tmp;
590
591
592 Cerr << "Permutation of local nodes "<<i2b<<" and "<<i3<<" on the element " <<num_poly<<" prod "<<prod_<<finl;
593
594 }
595 return 1;
596 }
597 return 0;
598}
599
600template <typename _SIZE_>
601void Hexaedre_VEF_32_64<_SIZE_>::reordonner()
602{
603 /*
604 22/12/04 : Generalisation de la methode (mais c'est plus lent!)
605 Numerotation d'un hexa VEF:
606 6--------7
607 | \ |\
608 4 2-----3 5
609 \ / \ |
610 0---------1
611 On boucle sur chacune des faces et on fait le produit scalaire 03.12/|03||12|
612 et on le compare a 02.13/|02||13|
613 S'il est plus grand on intervertit les sommets 2 et 3.
614 L'algorithme fonctionne si 0123 et 4567 sont deja sur une meme face
615 */
616 Domaine_t& dom=this->mon_dom.valeur();
617 IntTab_t& elem=dom.les_elems();
618 int_t nb_elem=dom.nb_elem();
619 DoubleTab_t& coord=dom.les_sommets();
620 // tableau contenant les permutations valides d'un quadrangle
621 IntTab perm(8,4);
622 perm( 0 , 0 )= 0 ;
623 perm( 0 , 1 )= 1 ;
624 perm( 0 , 2 )= 2 ;
625 perm( 0 , 3 )= 3 ;
626 perm( 1 , 0 )= 0 ;
627 perm( 1 , 1 )= 2 ;
628 perm( 1 , 2 )= 1 ;
629 perm( 1 , 3 )= 3 ;
630 perm( 2 , 0 )= 1 ;
631 perm( 2 , 1 )= 0 ;
632 perm( 2 , 2 )= 3 ;
633 perm( 2 , 3 )= 2 ;
634 perm( 3 , 0 )= 1 ;
635 perm( 3 , 1 )= 3 ;
636 perm( 3 , 2 )= 0 ;
637 perm( 3 , 3 )= 2 ;
638 perm( 4 , 0 )= 2 ;
639 perm( 4 , 1 )= 0 ;
640 perm( 4 , 2 )= 3 ;
641 perm( 4 , 3 )= 1 ;
642 perm( 5 , 0 )= 2 ;
643 perm( 5 , 1 )= 3 ;
644 perm( 5 , 2 )= 0 ;
645 perm( 5 , 3 )= 1 ;
646 perm( 6 , 0 )= 3 ;
647 perm( 6 , 1 )= 1 ;
648 perm( 6 , 2 )= 2 ;
649 perm( 6 , 3 )= 0 ;
650 perm( 7 , 0 )= 3 ;
651 perm( 7 , 1 )= 2 ;
652 perm( 7 , 2 )= 1 ;
653 perm( 7 , 3 )= 0 ;
654
655 DoubleTab v(3,3);
656 ArrOfDouble prod_(3);
657 ArrOfDouble prod_v(3);
658 ArrOfDouble dist(8);
659 DoubleTab prod_v2(3,3);
660
661 SmallArrOfTID_t sa(8);
662 for (int_t num_poly=0; num_poly<nb_elem; num_poly++)
663 {
664 int permutations=-1;
665 // int niter=0;
666 // on commence par remettre correctement les 2 plans opposes
667 permutations+=reordonne(0,1,2,3,coord,elem,num_poly,0, v,prod_, prod_v, dist, prod_v2);
668 permutations+=reordonne(4,5,6,7,coord,elem,num_poly,0, v,prod_, prod_v, dist, prod_v2);
669 for (int i=0; i<8; i++)
670 sa[i]=elem(num_poly,i);
671 // on teste les 8 possibilites
672 int j;
673 int n=0;
674 int perm_valid=-1;
675 for ( j=0; j<8; j++)
676 {
677 // Cerr<<" debut test "<<j<<finl;
678 int test=1;
679 for (int i=4; i<8; i++)
680 elem(num_poly,i)=sa[perm(j,i-4)+4];
681 // on regarde si les 6 plans sont valides
682 permutations=0;
683 permutations+=reordonne(2,0,3,1,coord,elem,num_poly,test, v,prod_, prod_v, dist, prod_v2);
684 permutations+=reordonne(3,1,7,5,coord,elem,num_poly,test, v,prod_, prod_v, dist, prod_v2);
685 permutations+=reordonne(7,5,6,4,coord,elem,num_poly,test, v,prod_, prod_v, dist, prod_v2);
686 permutations+=reordonne(6,4,2,0,coord,elem,num_poly,test, v,prod_, prod_v, dist, prod_v2);
687 permutations+=reordonne(2,3,6,7,coord,elem,num_poly,test, v,prod_, prod_v, dist, prod_v2);
688 permutations+=reordonne(4,5,0,1,coord,elem,num_poly,test, v,prod_, prod_v, dist, prod_v2);
689 if (permutations==0)
690 {
691 if (n>0)
692 Cerr<<" former permutation "<< perm_valid<<" new permutation "<<j<<finl;
693 perm_valid=j;
694 n++;
695 }
696 }
697 if (n==0)
698 {
699 Cerr << "Failure of the algorithm in Hexaedre_VEF_32_64<_SIZE_>::reordonner" << finl;
700 Cerr << "Here is the connectivity element-node of the cell " << num_poly << " that raises problem:" << finl;
701 for (int i=0; i<8; i++)
702 Cerr << sa[i] << " ";
703 Cerr << finl;
704 Process::exit();
705 }
706 else if (n>1)
707 {
708 Cerr<<"Failure of the algorithm in Hexaedre_VEF_32_64<_SIZE_>::reordonner"<<finl;
709 Cerr<<"many permutations possibles"<<finl;
710 Cerr << "Here is the connectivity element-node of the cell " << num_poly << " that raises problem:" << finl;
711 for (int i=0; i<8; i++)
712 Cerr << sa[i] << " ";
713 Cerr << finl;
714 Cerr << "Here is the new connectivity element-node of the cell " << num_poly << " that raises problem:" << finl;
715 for (int i=0; i<8; i++)
716 Cerr << elem(num_poly,i) << " ";
717 Cerr << finl;
718 // Process::exit();
719 }
720 {
721 for (int i=4; i<8; i++)
722 elem(num_poly,i)=sa[perm(perm_valid,i-4)+4];
723
724 if(perm_valid!=0)
725 Cerr<<" we carried out the permutation "<<perm_valid<<" on the element "<<num_poly<<finl;
726 }
727 }
728}
729
730static int faces_sommets_hexa_vef[6][4] =
731{
732 { 0, 2, 4, 6 },
733 { 0, 1, 4, 5 },
734 { 0, 1, 2, 3 },
735 { 1, 3, 5, 7 },
736 { 2, 3, 6, 7 },
737 { 4, 5, 6, 7 }
738};
739
740/*! @brief voir ElemGeomBase::get_tab_faces_sommets_locaux
741 *
742 */
743template <typename _SIZE_>
744int Hexaedre_VEF_32_64<_SIZE_>::get_tab_faces_sommets_locaux(IntTab& faces_som_local) const
745{
746 faces_som_local.resize(6,4);
747 for (int i=0; i<6; i++)
748 for (int j=0; j<4; j++)
749 faces_som_local(i,j) = faces_sommets_hexa_vef[i][j];
750 return 1;
751}
752
753
754template class Hexaedre_VEF_32_64<int>;
755#if INT_is_64_ == 2
756template class Hexaedre_VEF_32_64<trustIdType>;
757#endif
758
Domaine_32_64::nb_elem_tot
int_t nb_elem_tot() const
Definition Domaine.h:132
Domaine_32_64::les_sommets
DoubleTab_t & les_sommets()
Definition Domaine.h:113
Domaine_32_64::les_elems
IntTab_t & les_elems()
Definition Domaine.h:129
Domaine_32_64::nb_elem
int_t nb_elem() const
Definition Domaine.h:131
Domaine_32_64::coord
double coord(int_t i, int j) const
Definition Domaine.h:110
Domaine_32_64::sommet_elem
int_t sommet_elem(int_t i, int j) const
Renvoie le numero (global) du j-ieme sommet du i-ieme element.
Definition Domaine.h:136
Elem_geom_base_32_64
Classe Elem_geom_base Cette classe est la classe de base pour la definition d'elements.
Definition Elem_geom_base.h:35
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Hexaedre_VEF_32_64
Classe Hexaedre_VEF Cette classe represente l'element geometrique Hexaedre_VEF.
Definition Hexaedre_VEF.h:31
Hexaedre_VEF_32_64::get_tab_faces_sommets_locaux
int get_tab_faces_sommets_locaux(IntTab &faces_som_local) const override
voir ElemGeomBase::get_tab_faces_sommets_locaux
Definition Hexaedre_VEF.cpp:744
Hexaedre_VEF_32_64::int_t
_SIZE_ int_t
Definition Hexaedre_VEF.h:36
Hexaedre_VEF_32_64::nb_faces
int nb_faces(int=0) const override
Renvoie le nombre de faces du type specifie que possede l'element geometrique.
Definition Hexaedre_VEF.h:107
Hexaedre_VEF_32_64::contient
int contient(const ArrOfDouble &pos, int_t elem) const override
Renvoie 1 si l'element ielem du domaine associe a l'element geometrique contient le point.
Definition Hexaedre_VEF.cpp:254
Hexaedre_VEF_32_64::DoubleTab_t
DoubleTab_T< _SIZE_ > DoubleTab_t
Definition Hexaedre_VEF.h:40
Hexaedre_VEF_32_64::IntTab_t
IntTab_T< _SIZE_ > IntTab_t
Definition Hexaedre_VEF.h:37
Hexaedre_VEF_32_64::nom_lml
const Nom & nom_lml() const override
Renvoie le nom LML d'un Hexaedre_VEF = "HEXA8".
Definition Hexaedre_VEF.cpp:50
Hexaedre_VEF_32_64::DoubleVect_t
DoubleVect_T< _SIZE_ > DoubleVect_t
Definition Hexaedre_VEF.h:39
Hexaedre_VEF_32_64::reordonner
void reordonner() override
Definition Hexaedre_VEF.cpp:601
Hexaedre_VEF_32_64::calculer_volumes
void calculer_volumes(DoubleVect_t &vols) const override
Calcule les volumes des elements du domaine associe.
Definition Hexaedre_VEF.cpp:327
Hexaedre_VEF_32_64::SmallArrOfTID_t
SmallArrOfTID_T< _SIZE_ > SmallArrOfTID_t
Definition Hexaedre_VEF.h:38
Hexaedre_VEF_32_64::Domaine_t
Domaine_32_64< _SIZE_ > Domaine_t
Definition Hexaedre_VEF.h:41
Nom
class Nom Une chaine de caractere pour nommer les objets de TRUST
Definition Nom.h:31
Objet_U::readOn
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
Definition Objet_U.cpp:293
Objet_U::precision_geom
static double precision_geom
Definition Objet_U.h:86
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52
TRUSTArray::size_array
_SIZE_ size_array() const
Definition TRUSTArray.tpp:187
TRUSTTab::resize
void resize(_SIZE_ n, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTTab.tpp:469
TRUSTVect::size_totale
_SIZE_ size_totale() const
Definition TRUSTVect.tpp:61