TrioCFD 1.9.9_beta
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 Does nothing.
23 *
24 * @param (Sortie& s) an output stream
25 * @return (Sortie&) the output stream
26 */
27template <typename _SIZE_>
28Sortie& Hexaedre_VEF_32_64<_SIZE_>::printOn(Sortie& s ) const
29{
30 return s;
31}
32
33
34/*! @brief Does nothing.
35 *
36 * @param (Entree& s) an input stream
37 * @return (Entree&) the input stream
38 */
39template <typename _SIZE_>
40Entree& Hexaedre_VEF_32_64<_SIZE_>::readOn(Entree& s )
41{
42 return s;
43}
44
45/*! @brief Returns the LML name of a Hexaedre_VEF = "HEXA8".
46 *
47 * @return (Nom&) always equal to "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 // for 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 // for 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 // check if we are not at the boundary
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: the mixed products (som0, som1, som2,som3) and (pos, som1, som2,som3) must have the same sign
172 //construction of vectors
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 Returns 1 if element ielem of the domain associated with the geometric element contains the point
245 *
246 * whose coordinates are specified by parameter "pos".
247 * Returns 0 otherwise.
248 *
249 * @param (DoubleVect& pos) coordinates of the point to locate
250 * @param (int ielem) the index of the domain element in which to search for the point.
251 * @return (int) 1 if the point with the specified coordinates belongs to element ielem, 0 otherwise
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 check if the point is indeed between the 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 Returns 1 if the vertices specified by parameter "pos" are the vertices of element "element" of the domain associated with
298 *
299 * the geometric element.
300 *
301 * @param (IntVect& pos) the vertex indices to compare with those of element "element"
302 * @param (int element) the index of the domain element whose vertices are to be compared
303 * @return (int) 1 if the vertices passed as parameter are those of the specified element, 0 otherwise
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 Computes the volumes of the elements of the associated domain.
323 *
324 * @param (DoubleVect& volumes) the vector containing the volume values of the domain elements
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 // initialize the volume
348 // volumes[num_poly] =0;
349 // for each element, split it into tetrahedra.
350 // Compute the volume as the sum of the volumes of the tetrahedra.
351
352 // define the vertices of the element
353 // retrieve the center of the 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 // build the working array face_sommet_global
369 // faces will be ordered as Xmin, Xmax, Ymin, Ymax, Zmin, Zmax
370 // always preserving the TRUST numbering
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 // vertex 0 will always be the center of gravity
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 << "the processed poly index " << num_poly << finl;
407 // Cerr << "the center of gravity " << x0 << " " << y0 << " " << z0 << finl;
408
409 volume =0;
410
411
412 // Split into tetrahedra.
413 // for each face of the hexahedron (0 to 7)
414 for ( int k=0 ; k<nb_faces(); k++)
415 {
416 // retrieve the vertices of the current face
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 << " the vertices of face " << k << " = " << s1 << " " << s2 << " " << s3 << " " << s4 << finl;
423 // Cerr << "coord of s1 " << dom.coord(s1,0) << " " << dom.coord(s1,1) << " " << dom.coord(s1,2) << finl;
424 // Cerr << "coord of s2 " << dom.coord(s2,0) << " " << dom.coord(s2,1) << " " << dom.coord(s2,2) << finl;
425 // Cerr << "coord of s3 " << dom.coord(s3,0) << " " << dom.coord(s3,1) << " " << dom.coord(s3,2) << finl;
426 // Cerr << "coord of 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 // first check that we actually have a plane
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 // vertex closest to the plane
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 // not a plane
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 // should we test if the 4 other points are on the same side ???
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 // we have two negative products, do nothing
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: Generalization of the method (but slower!)
605 Numbering of a VEF hexahedron:
606 6--------7
607 | \ |\
608 4 2-----3 5
609 \ / \ |
610 0---------1
611 Loop over each face and compute the dot product 03.12/|03||12|
612 and compare it with 02.13/|02||13|
613 If it is larger, swap vertices 2 and 3.
614 The algorithm works if 0123 and 4567 are already on the same 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 // array containing valid permutations of a 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 // start by correctly repositioning the 2 opposite planes
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 // test the 8 possibilities
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 // check if the 6 planes are valid
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
Returns the (global) number of the j-th vertex of the i-th element.
Definition Domaine.h:136
Elem_geom_base_32_64
Class Elem_geom_base This class is the base class for the definition of 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
Class Hexaedre_VEF: represents the geometric element 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
Returns the number of faces of the specified type that the geometric element has.
Definition Hexaedre_VEF.h:107
Hexaedre_VEF_32_64::contient
int contient(const ArrOfDouble &pos, int_t elem) const override
Returns 1 if element ielem of the domain associated with the geometric element contains the 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
Returns the LML name of a 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
Computes the volumes of the elements of the associated domain.
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: a character string for naming TRUST objects.
Definition Nom.h:31
Objet_U::readOn
virtual Entree & readOn(Entree &)
Reads an Objet_U from an input stream. Virtual method to override.
Definition Objet_U.cpp:289
Objet_U::precision_geom
static double precision_geom
Definition Objet_U.h:81
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Writes the object to an output stream. Virtual method to override.
Definition Objet_U.cpp:278
Process::exit
static void exit(int exit_code=-1)
Exit routine for TRUST within a Kokkos region.
Definition Process.cpp:466
Sortie
Base class for output streams.
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