TrioCFD 1.9.9_beta
TrioCFD documentation
Loading...
Searching...
No Matches
Hexa_VEF.cpp
1/****************************************************************************
2* Copyright (c) 2026, 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 <Hexa_VEF.h>
17#include <Domaine.h>
18#include <Domaine_VEF.h>
19
20Implemente_instanciable_sans_constructeur(Hexa_VEF,"Hexa_VEF",Elem_VEF_base);
21
22// printOn and readOn
23
24Sortie& Hexa_VEF::printOn(Sortie& s ) const
25{
26 return s << que_suis_je() << finl;
27}
28
29Entree& Hexa_VEF::readOn(Entree& s )
30{
31 return s ;
32}
33
34/*! @brief KEL_(0,fa7),KEL_(1,fa7) are the local indices of the 2 faces surrounding the sub-facet of local index fa7.
35 *
36 * The local index of the sub-facet fa7 is that of the vertex it belongs to.
37 *
38 */
39Hexa_VEF::Hexa_VEF()
40{
41 int tmp[4][12]=
42 {
43 {1, 3, 4, 0, 5, 5, 5, 5, 0, 1, 3, 4},
44 {2, 2, 2, 2, 1, 3, 4, 0, 1, 3, 4, 0},
45 {0, 1, 3, 2, 4, 5, 7, 6, 4, 5, 7, 6},
46 {1, 3, 2, 0, 5, 7, 6, 4, 0, 1, 3, 2}
47 };
48 KEL_.resize(4,12);
49 for (int i=0; i<4; i++)
50 for (int j=0; j<12; j++)
51 KEL_(i,j)=tmp[i][j];
52}
53
54/*! @brief Fills the face_normales array in the Domaine_VEF.
55 *
56 */
57
58void Hexa_VEF::creer_face_normales(DoubleTab& Face_normales,
59 const IntTab& Face_sommets,
60 const IntTab& Face_voisins,
61 const IntTab& elem_faces,
62 const Domaine& domaine_geom) const
63{
64 const DoubleTab& les_coords = domaine_geom.coord_sommets();
65 int nb_face_tot = Face_normales.dimension_tot(0);
66 for (int num_Face=0; num_Face<nb_face_tot; num_Face++)
67 {
68 double x1, y1, z1, x2, y2, z2;
69 double nx, ny, nz;
70 int f0, no4;
71 int elem1;
72
73 // we have the 4 vertices of the face
74 int n0 = Face_sommets(num_Face, 0);
75 int n1 = Face_sommets(num_Face, 1);
76 int n2 = Face_sommets(num_Face, 2);
77 int n3 = Face_sommets(num_Face, 3);
78
79 // split the face into two triangles
80 // initialize Face_normales(num_face,dimension) to 0
81 // since we will sum the two triangle face normals
82 Face_normales(num_Face, 0) = 0;
83 Face_normales(num_Face, 1) = 0;
84 Face_normales(num_Face, 2) = 0;
85
86 // take vertices S1, S2 and S4
87
88 x1 = les_coords(n0, 0) - les_coords(n1, 0);
89 y1 = les_coords(n0, 1) - les_coords(n1, 1);
90 z1 = les_coords(n0, 2) - les_coords(n1, 2);
91
92 x2 = les_coords(n3, 0) - les_coords(n1, 0);
93 y2 = les_coords(n3, 1) - les_coords(n1, 1);
94 z2 = les_coords(n3, 2) - les_coords(n1, 2);
95
96 nx = (y1 * z2 - y2 * z1) / 2;
97 ny = (-x1 * z2 + x2 * z1) / 2;
98 nz = (x1 * y2 - x2 * y1) / 2;
99
100 // Orientation of the normal from elem1 to elem2:
101 // find the vertex of elem1 that does not lie on the face
102
103 elem1 = Face_voisins(num_Face, 0);
104
105 if ((f0 = elem_faces(elem1, 0)) == num_Face)
106 f0 = elem_faces(elem1, 1);
107
108 if ((no4 = Face_sommets(f0, 0)) != n0 && no4 != n1
109 && no4 != n2 && no4 != n3) { /* Do nothing */}
110 else if ((no4 = Face_sommets(f0, 1)) != n0 && no4 != n1
111 && no4 != n2 && no4 != n3) { /* Do nothing */}
112 else if ((no4 = Face_sommets(f0, 2)) != n0 && no4 != n1
113 && no4 != n2 && no4 != n3) { /* Do nothing */}
114 else
115 no4 = Face_sommets(f0, 3);
116
117 x1 = les_coords(no4, 0) - les_coords(n0, 0);
118 y1 = les_coords(no4, 1) - les_coords(n0, 1);
119 z1 = les_coords(no4, 2) - les_coords(n0, 2);
120
121 if ((nx * x1 + ny * y1 + nz * z1) > 0)
122 {
123 Face_normales(num_Face, 0) += -nx;
124 Face_normales(num_Face, 1) += -ny;
125 Face_normales(num_Face, 2) += -nz;
126 }
127 else
128 {
129 Face_normales(num_Face, 0) += nx;
130 Face_normales(num_Face, 1) += ny;
131 Face_normales(num_Face, 2) += nz;
132 }
133
134 // take vertices S1, S4 and S3
135
136 x1 = les_coords(n0, 0) - les_coords(n2, 0);
137 y1 = les_coords(n0, 1) - les_coords(n2, 1);
138 z1 = les_coords(n0, 2) - les_coords(n2, 2);
139
140 x2 = les_coords(n3, 0) - les_coords(n2, 0);
141 y2 = les_coords(n3, 1) - les_coords(n2, 1);
142 z2 = les_coords(n3, 2) - les_coords(n2, 2);
143
144 nx = (y1 * z2 - y2 * z1) / 2;
145 ny = (-x1 * z2 + x2 * z1) / 2;
146 nz = (x1 * y2 - x2 * y1) / 2;
147
148 // Orientation of the normal from elem1 to elem2:
149 // find the vertex of elem1 that does not lie on the face
150
151 elem1 = Face_voisins(num_Face, 0);
152
153 if ((f0 = elem_faces(elem1, 0)) == num_Face)
154 f0 = elem_faces(elem1, 1);
155
156 if ((no4 = Face_sommets(f0, 0)) != n0 && no4 != n1
157 && no4 != n2 && no4 != n3) { /* Do nothing */}
158 else if ((no4 = Face_sommets(f0, 1)) != n0 && no4 != n1
159 && no4 != n2 && no4 != n3) { /* Do nothing */}
160 else if ((no4 = Face_sommets(f0, 2)) != n0 && no4 != n1
161 && no4 != n2 && no4 != n3) { /* Do nothing */}
162 else
163 no4 = Face_sommets(f0, 3);
164
165 x1 = les_coords(no4, 0) - les_coords(n0, 0);
166 y1 = les_coords(no4, 1) - les_coords(n0, 1);
167 z1 = les_coords(no4, 2) - les_coords(n0, 2);
168
169 if ((nx * x1 + ny * y1 + nz * z1) > 0)
170 {
171 Face_normales(num_Face, 0) += -nx;
172 Face_normales(num_Face, 1) += -ny;
173 Face_normales(num_Face, 2) += -nz;
174 }
175 else
176 {
177 Face_normales(num_Face, 0) += nx;
178 Face_normales(num_Face, 1) += ny;
179 Face_normales(num_Face, 2) += nz;
180 }
181 }
182}
183
184
185/*! @brief Computes the facet normals for standard elements.
186 *
187 */
188void Hexa_VEF::creer_facette_normales(const Domaine_VEF& dom_VEF,
189 const IntVect& rang_elem_non_std) const
190{
191 const Domaine& domaine_geom = dom_VEF.domaine();
192 auto& facette_normales = const_cast<Domaine_VEF&>(dom_VEF).facette_normales();
193 const DoubleTab& les_coords = domaine_geom.coord_sommets();
194 const IntTab& les_Polys = domaine_geom.les_elems();
195 int nb_elem_tot = domaine_geom.nb_elem_tot();
196 int i, fa7;
197 int num_som[8];
198 double x[8][3];
199 double x_[4][3];
200 double xg[3];
201 double xj0[3];
202 double xjnum1[3];
203 double xjnum2[3];
204 double u[3];
205 double v[3];
206 double psc;
207 double pv[3];
208
209 if (facette_normales.dimension(0) != nb_elem_tot)
210 facette_normales.resize(nb_elem_tot,nb_facette(),3);
211
212 for(i=0; i<nb_elem_tot; i++)
213 {
214 if (rang_elem_non_std(i)==-1)
215 {
216 num_som[0]=les_Polys(i,0);
217 num_som[1]=les_Polys(i,1);
218 num_som[2]=les_Polys(i,2);
219 num_som[3]=les_Polys(i,3);
220 num_som[4]=les_Polys(i,4);
221 num_som[5]=les_Polys(i,5);
222 num_som[6]=les_Polys(i,6);
223 num_som[7]=les_Polys(i,7);
224
225 x[0][0]=les_coords(num_som[0],0);
226 x[0][1]=les_coords(num_som[0],1);
227 x[0][2]=les_coords(num_som[0],2);
228 x[1][0]=les_coords(num_som[1],0);
229 x[1][1]=les_coords(num_som[1],1);
230 x[1][2]=les_coords(num_som[1],2);
231 x[2][0]=les_coords(num_som[2],0);
232 x[2][1]=les_coords(num_som[2],1);
233 x[2][2]=les_coords(num_som[2],2);
234 x[3][0]=les_coords(num_som[3],0);
235 x[3][1]=les_coords(num_som[3],1);
236 x[3][2]=les_coords(num_som[3],2);
237 x[4][0]=les_coords(num_som[4],0);
238 x[4][1]=les_coords(num_som[4],1);
239 x[4][2]=les_coords(num_som[4],2);
240 x[5][0]=les_coords(num_som[5],0);
241 x[5][1]=les_coords(num_som[5],1);
242 x[5][2]=les_coords(num_som[5],2);
243 x[6][0]=les_coords(num_som[6],0);
244 x[6][1]=les_coords(num_som[6],1);
245 x[6][2]=les_coords(num_som[6],2);
246 x[7][0]=les_coords(num_som[7],0);
247 x[7][1]=les_coords(num_som[7],1);
248 x[7][2]=les_coords(num_som[7],2);
249
250 xg[0]=(x[0][0]+x[1][0]+x[2][0]+x[3][0]+x[4][0]+x[5][0]+x[6][0]+x[7][0])*0.125;
251 xg[1]=(x[0][1]+x[1][1]+x[2][1]+x[3][1]+x[4][1]+x[5][1]+x[6][1]+x[7][1])*0.125;
252 xg[2]=(x[0][2]+x[1][2]+x[2][2]+x[3][2]+x[4][2]+x[5][2]+x[6][2]+x[7][2])*0.125;
253
254
255 for (fa7=0; fa7<nb_facette(); fa7++)
256 {
257 // the sub-facet fa7 of a standard element has vertices:
258 // fa7 and G with coordinates xg:
259 u[0]=x[KEL_(2,fa7)][0]-xg[0];
260 u[1]=x[KEL_(2,fa7)][1]-xg[1];
261 u[2]=x[KEL_(2,fa7)][2]-xg[2];
262 v[0]=x[KEL_(3,fa7)][0]-xg[0];
263 v[1]=x[KEL_(3,fa7)][1]-xg[1];
264 v[2]=x[KEL_(3,fa7)][2]-xg[2];
265 prodvect(u,v,pv);
266
267 // Normal orientation:
268
269 // retrieve the face index, then the vertex indices, and compute the face center of gravity
270
271 int num1 = KEL_(0,fa7);
272 int num2 = KEL_(1,fa7);
273
274 switch (num1)
275 {
276 case 0:
277 {
278 num_som[0]=les_Polys(i,0);
279 num_som[1]=les_Polys(i,2);
280 num_som[2]=les_Polys(i,4);
281 num_som[3]=les_Polys(i,6);
282 }
283 break;
284 case 1:
285 {
286 num_som[0]=les_Polys(i,0);
287 num_som[1]=les_Polys(i,1);
288 num_som[2]=les_Polys(i,4);
289 num_som[3]=les_Polys(i,5);
290 }
291 break;
292 case 2:
293 {
294 num_som[0]=les_Polys(i,0);
295 num_som[1]=les_Polys(i,1);
296 num_som[2]=les_Polys(i,2);
297 num_som[3]=les_Polys(i,3);
298 }
299 break;
300 case 3:
301 {
302 num_som[0]=les_Polys(i,1);
303 num_som[1]=les_Polys(i,3);
304 num_som[2]=les_Polys(i,5);
305 num_som[3]=les_Polys(i,7);
306 }
307 break;
308 case 4:
309 {
310 num_som[0]=les_Polys(i,2);
311 num_som[1]=les_Polys(i,3);
312 num_som[2]=les_Polys(i,6);
313 num_som[3]=les_Polys(i,7);
314 }
315 break;
316 case 5:
317 {
318 num_som[0]=les_Polys(i,4);
319 num_som[1]=les_Polys(i,5);
320 num_som[2]=les_Polys(i,6);
321 num_som[3]=les_Polys(i,7);
322 }
323 break;
324 }
325
326 x_[0][0]=les_coords(num_som[0],0);
327 x_[0][1]=les_coords(num_som[0],1);
328 x_[0][2]=les_coords(num_som[0],2);
329 x_[1][0]=les_coords(num_som[1],0);
330 x_[1][1]=les_coords(num_som[1],1);
331 x_[1][2]=les_coords(num_som[1],2);
332 x_[2][0]=les_coords(num_som[2],0);
333 x_[2][1]=les_coords(num_som[2],1);
334 x_[2][2]=les_coords(num_som[2],2);
335 x_[3][0]=les_coords(num_som[3],0);
336 x_[3][1]=les_coords(num_som[3],1);
337 x_[3][2]=les_coords(num_som[3],2);
338
339 xjnum1[0]=(x_[0][0]+x_[1][0]+x_[2][0]+x_[3][0])*0.25;
340 xjnum1[1]=(x_[0][1]+x_[1][1]+x_[2][1]+x_[3][1])*0.25;
341 xjnum1[2]=(x_[0][2]+x_[1][2]+x_[2][2]+x_[3][2])*0.25;
342
343 switch (num2)
344 {
345 case 0:
346 {
347 num_som[0]=les_Polys(i,0);
348 num_som[1]=les_Polys(i,2);
349 num_som[2]=les_Polys(i,4);
350 num_som[3]=les_Polys(i,6);
351 }
352 break;
353 case 1:
354 {
355 num_som[0]=les_Polys(i,0);
356 num_som[1]=les_Polys(i,1);
357 num_som[2]=les_Polys(i,4);
358 num_som[3]=les_Polys(i,5);
359 }
360 break;
361 case 2:
362 {
363 num_som[0]=les_Polys(i,0);
364 num_som[1]=les_Polys(i,1);
365 num_som[2]=les_Polys(i,2);
366 num_som[3]=les_Polys(i,3);
367 }
368 break;
369 case 3:
370 {
371 num_som[0]=les_Polys(i,1);
372 num_som[1]=les_Polys(i,3);
373 num_som[2]=les_Polys(i,5);
374 num_som[3]=les_Polys(i,7);
375 }
376 break;
377 case 4:
378 {
379 num_som[0]=les_Polys(i,2);
380 num_som[1]=les_Polys(i,3);
381 num_som[2]=les_Polys(i,6);
382 num_som[3]=les_Polys(i,7);
383 }
384 break;
385 case 5:
386 {
387 num_som[0]=les_Polys(i,4);
388 num_som[1]=les_Polys(i,5);
389 num_som[2]=les_Polys(i,6);
390 num_som[3]=les_Polys(i,7);
391 }
392 break;
393 }
394
395 x_[0][0]=les_coords(num_som[0],0);
396 x_[0][1]=les_coords(num_som[0],1);
397 x_[0][2]=les_coords(num_som[0],2);
398 x_[1][0]=les_coords(num_som[1],0);
399 x_[1][1]=les_coords(num_som[1],1);
400 x_[1][2]=les_coords(num_som[1],2);
401 x_[2][0]=les_coords(num_som[2],0);
402 x_[2][1]=les_coords(num_som[2],1);
403 x_[2][2]=les_coords(num_som[2],2);
404 x_[3][0]=les_coords(num_som[3],0);
405 x_[3][1]=les_coords(num_som[3],1);
406 x_[3][2]=les_coords(num_som[3],2);
407
408 xjnum2[0]=(x_[0][0]+x_[1][0]+x_[2][0]+x_[3][0])*0.25;
409 xjnum2[1]=(x_[0][1]+x_[1][1]+x_[2][1]+x_[3][1])*0.25;
410 xjnum2[2]=(x_[0][2]+x_[1][2]+x_[2][2]+x_[3][2])*0.25;
411
412 xj0[0]= xjnum1[0]-xjnum2[0];
413 xj0[1]= xjnum1[1]-xjnum2[1];
414 xj0[2]= xjnum1[2]-xjnum2[2];
415
416 psc=xj0[0]*pv[0] + xj0[1]*pv[1] + xj0[2]*pv[2] ;
417
418 if (psc < 0)
419 {
420 facette_normales(i,fa7,0) = pv[0]/2;
421 facette_normales(i,fa7,1) = pv[1]/2;
422 facette_normales(i,fa7,2) = pv[2]/2;
423 }
424 else
425 {
426 facette_normales(i,fa7,0)= -pv[0]/2;
427 facette_normales(i,fa7,1)= -pv[1]/2;
428 facette_normales(i,fa7,2)= -pv[2]/2;
429 }
430 }
431 }
432 }
433}
434
435void Hexa_VEF::creer_normales_facettes_Cl(DoubleTab& normales_facettes_Cl,
436 int fa7,
437 int num_elem,const DoubleTab& x,
438 const DoubleVect& xg, const
439 Domaine& domaine_geom) const
440{
441 double u[3];
442 double v[3];
443 double xj0[3];
444 double xjnum1[3];
445 double xjnum2[3];
446
447 int num_som[8];
448 double x_[4][3];
449 double pv[3];
450 double psc;
451
452 const DoubleTab& les_coords = domaine_geom.coord_sommets();
453 const IntTab& les_Polys = domaine_geom.les_elems();
454
455 int i2 = KEL_(2,fa7);
456 int i3 = KEL_(3,fa7);
457
458 u[0]=x(i2,0)-xg[0];
459 u[1]=x(i2,1)-xg[1];
460 u[2]=x(i2,2)-xg[2];
461
462 v[0]=x(i3,0)-xg[0];
463 v[1]=x(i3,1)-xg[1];
464 v[2]=x(i3,2)-xg[2];
465
466 prodvect(u,v,pv);
467
468 // Normal orientation:
469
470 // retrieve the face index, then the vertex indices, and compute the face center of gravity
471
472 int num1 = KEL_(0,fa7);
473 int num2 = KEL_(1,fa7);
474
475 switch (num1)
476 {
477 case 0:
478 {
479 num_som[0]=les_Polys(num_elem,0);
480 num_som[1]=les_Polys(num_elem,2);
481 num_som[2]=les_Polys(num_elem,4);
482 num_som[3]=les_Polys(num_elem,6);
483 }
484 break;
485 case 1:
486 {
487 num_som[0]=les_Polys(num_elem,0);
488 num_som[1]=les_Polys(num_elem,1);
489 num_som[2]=les_Polys(num_elem,4);
490 num_som[3]=les_Polys(num_elem,5);
491 }
492 break;
493 case 2:
494 {
495 num_som[0]=les_Polys(num_elem,0);
496 num_som[1]=les_Polys(num_elem,1);
497 num_som[2]=les_Polys(num_elem,2);
498 num_som[3]=les_Polys(num_elem,3);
499 }
500 break;
501 case 3:
502 {
503 num_som[0]=les_Polys(num_elem,1);
504 num_som[1]=les_Polys(num_elem,3);
505 num_som[2]=les_Polys(num_elem,5);
506 num_som[3]=les_Polys(num_elem,7);
507 }
508 break;
509 case 4:
510 {
511 num_som[0]=les_Polys(num_elem,2);
512 num_som[1]=les_Polys(num_elem,3);
513 num_som[2]=les_Polys(num_elem,6);
514 num_som[3]=les_Polys(num_elem,7);
515 }
516 break;
517 case 5:
518 {
519 num_som[0]=les_Polys(num_elem,4);
520 num_som[1]=les_Polys(num_elem,5);
521 num_som[2]=les_Polys(num_elem,6);
522 num_som[3]=les_Polys(num_elem,7);
523 }
524 break;
525 default:
526 {
527 exit();
528 num_som[0]=-1;
529 num_som[1]=-1;
530 num_som[2]=-1;
531 num_som[3]=-1;
532 }
533 break;
534 }
535
536 x_[0][0]=les_coords(num_som[0],0);
537 x_[0][1]=les_coords(num_som[0],1);
538 x_[0][2]=les_coords(num_som[0],2);
539 x_[1][0]=les_coords(num_som[1],0);
540 x_[1][1]=les_coords(num_som[1],1);
541 x_[1][2]=les_coords(num_som[1],2);
542 x_[2][0]=les_coords(num_som[2],0);
543 x_[2][1]=les_coords(num_som[2],1);
544 x_[2][2]=les_coords(num_som[2],2);
545 x_[3][0]=les_coords(num_som[3],0);
546 x_[3][1]=les_coords(num_som[3],1);
547 x_[3][2]=les_coords(num_som[3],2);
548
549 xjnum1[0]=(x_[0][0]+x_[1][0]+x_[2][0]+x_[3][0])*0.25;
550 xjnum1[1]=(x_[0][1]+x_[1][1]+x_[2][1]+x_[3][1])*0.25;
551 xjnum1[2]=(x_[0][2]+x_[1][2]+x_[2][2]+x_[3][2])*0.25;
552
553
554 switch (num2)
555 {
556 case 0:
557 {
558 num_som[0]=les_Polys(num_elem,0);
559 num_som[1]=les_Polys(num_elem,2);
560 num_som[2]=les_Polys(num_elem,4);
561 num_som[3]=les_Polys(num_elem,6);
562 }
563 break;
564 case 1:
565 {
566 num_som[0]=les_Polys(num_elem,0);
567 num_som[1]=les_Polys(num_elem,1);
568 num_som[2]=les_Polys(num_elem,4);
569 num_som[3]=les_Polys(num_elem,5);
570 }
571 break;
572 case 2:
573 {
574 num_som[0]=les_Polys(num_elem,0);
575 num_som[1]=les_Polys(num_elem,1);
576 num_som[2]=les_Polys(num_elem,2);
577 num_som[3]=les_Polys(num_elem,3);
578 }
579 break;
580 case 3:
581 {
582 num_som[0]=les_Polys(num_elem,1);
583 num_som[1]=les_Polys(num_elem,3);
584 num_som[2]=les_Polys(num_elem,5);
585 num_som[3]=les_Polys(num_elem,7);
586 }
587 break;
588 case 4:
589 {
590 num_som[0]=les_Polys(num_elem,2);
591 num_som[1]=les_Polys(num_elem,3);
592 num_som[2]=les_Polys(num_elem,6);
593 num_som[3]=les_Polys(num_elem,7);
594 }
595 break;
596 case 5:
597 {
598 num_som[0]=les_Polys(num_elem,4);
599 num_som[1]=les_Polys(num_elem,5);
600 num_som[2]=les_Polys(num_elem,6);
601 num_som[3]=les_Polys(num_elem,7);
602 }
603 break;
604 }
605
606
607 x_[0][0]=les_coords(num_som[0],0);
608 x_[0][1]=les_coords(num_som[0],1);
609 x_[0][2]=les_coords(num_som[0],2);
610 x_[1][0]=les_coords(num_som[1],0);
611 x_[1][1]=les_coords(num_som[1],1);
612 x_[1][2]=les_coords(num_som[1],2);
613 x_[2][0]=les_coords(num_som[2],0);
614 x_[2][1]=les_coords(num_som[2],1);
615 x_[2][2]=les_coords(num_som[2],2);
616 x_[3][0]=les_coords(num_som[3],0);
617 x_[3][1]=les_coords(num_som[3],1);
618 x_[3][2]=les_coords(num_som[3],2);
619
620 xjnum2[0]=(x_[0][0]+x_[1][0]+x_[2][0]+x_[3][0])*0.25;
621 xjnum2[1]=(x_[0][1]+x_[1][1]+x_[2][1]+x_[3][1])*0.25;
622 xjnum2[2]=(x_[0][2]+x_[1][2]+x_[2][2]+x_[3][2])*0.25;
623
624 xj0[0]= xjnum1[0]-xjnum2[0];
625 xj0[1]= xjnum1[1]-xjnum2[1];
626 xj0[2]= xjnum1[2]-xjnum2[2];
627
628 psc=xj0[0]*pv[0] + xj0[1]*pv[1] + xj0[2]*pv[2] ;
629
630 if (psc < 0)
631 {
632 normales_facettes_Cl(num_elem,fa7,0) = -pv[0]/2;
633 normales_facettes_Cl(num_elem,fa7,1) = -pv[1]/2;
634 normales_facettes_Cl(num_elem,fa7,2) = -pv[2]/2;
635 }
636 else
637 {
638 normales_facettes_Cl(num_elem,fa7,0) = pv[0]/2;
639 normales_facettes_Cl(num_elem,fa7,1) = pv[1]/2;
640 normales_facettes_Cl(num_elem,fa7,2) = pv[2]/2;
641 }
642}
643
644/*! @brief Modifies the interlaced volumes for face j of a non-standard element.
645 *
646 */
647void Hexa_VEF::modif_volumes_entrelaces(int j,int elem,
648 const Domaine_VEF& le_dom_VEF,
649 DoubleVect& volumes_entrelaces_Cl,
650 int type_cl) const
651{
652 Cerr << "Hexa_VEF::modif_volumes_entrelaces() does nothing for now" << finl;
653
654}
655
656/*! @brief Modifies the interlaced volumes for the joint face j of a non-standard element.
657 *
658 */
659void Hexa_VEF::modif_volumes_entrelaces_faces_joints(int j,int elem,
660 const Domaine_VEF& le_dom_VEF,
661 DoubleVect& volumes_entrelaces_Cl,
662 int type_cl) const
663{
664
665
666}
667
668/*! @brief
669 *
670 */
671void Hexa_VEF::calcul_vc(const ArrOfInt& Face,ArrOfDouble& vc,
672 const ArrOfDouble& vs,const DoubleTab& vsom,
673 const Champ_Inc_base& vitesse,int type_cl, const DoubleVect& porosite_face) const
674{
675
676 vc[0] = vs[0]/6;
677 vc[1] = vs[1]/6;
678 vc[2] = vs[2]/6;
679}
680
681/*! @brief Computes the coordinates xg of the center of a non-standard element.
682 *
683 * Also computes idirichlet = number of Dirichlet faces of the element.
684 * If idirichlet=2, n1 is the index of the vertex coinciding with G.
685 *
686 */
687void Hexa_VEF::calcul_xg(DoubleVect& xg, const DoubleTab& x,
688 const int type_elem_Cl,int& idirichlet,int& n1,int& ,int& ) const
689{
690 int dim=xg.size();
691 for (int j=0; j<dim; j++)
692 xg[j]=(x(0,j)+x(1,j)+x(2,j)+x(3,j)+x(4,j)+x(5,j)+x(6,j)+x(7,j))*0.125;
693 idirichlet=0;
694}
695
696
697/*! @brief Modifies normales_facettes_Cl when idirichlet=3.
698 *
699 * idirichlet = number of Dirichlet faces of the element.
700 * If idirichlet=3, n1 is the index of the vertex coinciding with G.
701 *
702 */
703void Hexa_VEF::modif_normales_facettes_Cl(DoubleTab& normales_facettes_Cl,
704 int fa7,int num_elem,
705 int idirichlet,int n1,int ,int ) const
706{
707 switch (idirichlet)
708 {
709 case 0:
710 break;
711 case 1:
712 break;
713 case 2:
714 break;
715 case 3:
716 {
717 // normales_facettes_Cl(num_elem,n1,0) = 0;
718 // normales_facettes_Cl(num_elem,n1,1) = 0;
719 break;
720 }
721 }
722}
Champ_Inc_base
Class Champ_Inc_base.
Definition Champ_Inc_base.h:53
Domaine_32_64::nb_elem_tot
int_t nb_elem_tot() const
Definition Domaine.h:132
Domaine_32_64::les_elems
IntTab_t & les_elems()
Definition Domaine.h:129
Domaine_32_64::coord_sommets
const DoubleTab_t & coord_sommets() const
Definition Domaine.h:112
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:53
Domaine_dis_base::domaine
const Domaine & domaine() const
Definition Domaine_dis_base.h:46
Elem_VEF_base
Definition Elem_VEF_base.h:27
Elem_VEF_base::KEL_
IntTab KEL_
Definition Elem_VEF_base.h:41
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Hexa_VEF
Definition Hexa_VEF.h:23
Hexa_VEF::Hexa_VEF
Hexa_VEF()
KEL_(0,fa7),KEL_(1,fa7) are the local indices of the 2 faces surrounding the sub-facet of local index...
Definition Hexa_VEF.cpp:39
Hexa_VEF::modif_normales_facettes_Cl
void modif_normales_facettes_Cl(DoubleTab &, int, int, int, int, int, int) const override
Modifies normales_facettes_Cl when idirichlet=3.
Definition Hexa_VEF.cpp:703
Hexa_VEF::nb_facette
int nb_facette() const override
Definition Hexa_VEF.h:29
Hexa_VEF::creer_normales_facettes_Cl
void creer_normales_facettes_Cl(DoubleTab &, int, int, const DoubleTab &, const DoubleVect &, const Domaine &) const override
Definition Hexa_VEF.cpp:435
Hexa_VEF::creer_facette_normales
void creer_facette_normales(const Domaine_VEF &, const IntVect &) const override
Computes the facet normals for standard elements.
Definition Hexa_VEF.cpp:188
Hexa_VEF::calcul_vc
void calcul_vc(const ArrOfInt &, ArrOfDouble &, const ArrOfDouble &, const DoubleTab &, const Champ_Inc_base &, int, const DoubleVect &) const override
Definition Hexa_VEF.cpp:671
Hexa_VEF::creer_face_normales
void creer_face_normales(DoubleTab &, const IntTab &, const IntTab &, const IntTab &, const Domaine &) const override
Fills the face_normales array in the Domaine_VEF.
Definition Hexa_VEF.cpp:58
Hexa_VEF::calcul_xg
void calcul_xg(DoubleVect &, const DoubleTab &, const int, int &, int &, int &, int &) const override
Computes the coordinates xg of the center of a non-standard element.
Definition Hexa_VEF.cpp:687
Hexa_VEF::modif_volumes_entrelaces_faces_joints
void modif_volumes_entrelaces_faces_joints(int, int, const Domaine_VEF &, DoubleVect &, int) const override
Modifies the interlaced volumes for the joint face j of a non-standard element.
Definition Hexa_VEF.cpp:659
Hexa_VEF::modif_volumes_entrelaces
void modif_volumes_entrelaces(int, int, const Domaine_VEF &, DoubleVect &, int) const override
Modifies the interlaced volumes for face j of a non-standard element.
Definition Hexa_VEF.cpp:647
Objet_U::que_suis_je
const Nom & que_suis_je() const
Returns the string identifying the class.
Definition Objet_U.cpp:104
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::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
TRUSTTab::dimension_tot
_SIZE_ dimension_tot(int) const override
Definition TRUSTTab.tpp:160
TRUSTVect::size
_SIZE_ size() const
Definition TRUSTVect.tpp:45