TrioCFD 1.9.9_beta
TrioCFD documentation
Loading...
Searching...
No Matches
Tri_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 <Tri_VEF.h>
17#include <Domaine.h>
18#include <Domaine_VEF.h>
19#include <Champ_P1NC.h>
20
21Implemente_instanciable_sans_constructeur(Tri_VEF,"Tri_VEF",Elem_VEF_base);
22
23// printOn and readOn
24
25
27{
28 return s << que_suis_je() << finl;
29}
30
32{
33 return s ;
34}
35
37{
38 /*! @brief KEL_(0,fa7),KEL_(1,fa7) are the local indices of the 2 faces surrounding facet with local index fa7.
39 *
40 */
41 int tmp[3][3]=
42 {
43 {1, 2, 0},
44 {2, 0, 1},
45 {0, 1, 2}
46 };
47 KEL_.resize(3,3);
48 for (int i=0; i<3; i++)
49 for (int j=0; j<3; j++)
50 KEL_(i,j)=tmp[i][j];
51}
52
53/*! @brief Fills the face_normales array in Domaine_VEF.
54 *
55 */
56void Tri_VEF::creer_face_normales(DoubleTab& Face_normales,
57 const IntTab& Face_sommets,
58 const IntTab& Face_voisins,
59 const IntTab& elem_faces,
60 const Domaine& domaine_geom) const
61{
62 const DoubleTab& les_coords = domaine_geom.coord_sommets();
63 int nb_face_tot = Face_normales.dimension_tot(0);
64 for (int num_Face=0; num_Face<nb_face_tot; num_Face++)
65 {
66 double x1, y1;
67 double nx, ny;
68 int no3;
69 int f0;
70 int n0 = Face_sommets(num_Face, 0);
71 int n1 = Face_sommets(num_Face, 1);
72 x1 = les_coords(n0, 0) - les_coords(n1, 0);
73 y1 = les_coords(n0, 1) - les_coords(n1, 1);
74 nx = -y1;
75 ny = x1;
76
77 // Orient the normal from elem1 towards elem2
78 // for that, find the node of elem1 that is not on the Face
79 int elem1 = Face_voisins(num_Face, 0);
80 if ((f0 = elem_faces(elem1, 0)) == num_Face)
81 f0 = elem_faces(elem1, 1);
82 if ((no3 = Face_sommets(f0, 0)) != n0 && no3 != n1) { /* Do nothing */}
83 else
84 no3 = Face_sommets(f0, 1);
85
86 x1 = les_coords(no3, 0) - les_coords(n0, 0);
87 y1 = les_coords(no3, 1) - les_coords(n0, 1);
88
89 if ((nx * x1 + ny * y1) > 0)
90 {
91 Face_normales(num_Face, 0) = -nx;
92 Face_normales(num_Face, 1) = -ny;
93 }
94 else
95 {
96 Face_normales(num_Face, 0) = nx;
97 Face_normales(num_Face, 1) = ny;
98 }
99 }
100}
101/*! @brief Computes the facet normals for standard elements.
102 *
103 */
105 const IntVect& rang_elem_non_std) const
106{
107 const Domaine& domaine_geom = dom_VEF.domaine();
108 auto& facette_normales = const_cast<Domaine_VEF&>(dom_VEF).facette_normales();
109 const DoubleTab& les_coords = domaine_geom.coord_sommets();
110 const IntTab& les_Polys = domaine_geom.les_elems();
111 int nb_elem_tot = domaine_geom.nb_elem_tot();
112
113 int i, fa7;
114 int i0,i1;
115 int num_som[3];
116 double x[3][2];
117 double xg[2];
118 double xj0[2];
119 double u[2];
120 double v[2];
121 double psc;
122
123 if (facette_normales.dimension(0) != nb_elem_tot)
124 facette_normales.resize(nb_elem_tot,3,2);
125
126 for(i=0; i<nb_elem_tot; i++)
127 {
128 if (rang_elem_non_std(i)==-1)
129 {
130 num_som[0]=les_Polys(i,0);
131 num_som[1]=les_Polys(i,1);
132 num_som[2]=les_Polys(i,2);
133 x[0][0]=les_coords(num_som[0],0);
134 x[0][1]=les_coords(num_som[0],1);
135 x[1][0]=les_coords(num_som[1],0);
136 x[1][1]=les_coords(num_som[1],1);
137 x[2][0]=les_coords(num_som[2],0);
138 x[2][1]=les_coords(num_som[2],1);
139 xg[0]=(x[0][0]+x[1][0]+x[2][0])/3;
140 xg[1]=(x[0][1]+x[1][1]+x[2][1])/3;
141 for (fa7=0; fa7<3; fa7++)
142 {
143 // facet fa7 has vertices fa7 and "G" with coordinates xg
144 u[0]= x[fa7][0]-xg[0];
145 u[1]= x[fa7][1]-xg[1];
146 v[0]= -u[1];
147 v[1]= u[0];
148 i0 = KEL_(0,fa7);
149 i1 = KEL_(1,fa7);
150 // Orient normals:
151 xj0[0]= x[i0][0]-x[i1][0];
152 xj0[1]= x[i0][1]-x[i1][1];
153 psc=xj0[0]*v[0] + xj0[1]*v[1];
154 if (psc < 0)
155 {
156 facette_normales(i,fa7,0) = -v[0];
157 facette_normales(i,fa7,1) = -v[1];
158 }
159 else
160 {
161 facette_normales(i,fa7,0) = v[0];
162 facette_normales(i,fa7,1) = v[1];
163 }
164 }
165 }
166 }
167}
168
169/*! @brief Fills the normales_facettes_Cl array in Domaine_Cl_VEF for facet fa7 of element num_elem.
170 *
171 */
172void Tri_VEF::creer_normales_facettes_Cl(DoubleTab& normales_facettes_Cl,
173 int fa7,
174 int num_elem,const DoubleTab& x,
175 const DoubleVect& xg, const Domaine& domaine_geom) const
176{
177 double u[2];
178 double v[2];
179 double xj0[2];
180 double psc;
181
182 u[0]= x(fa7,0)-xg[0];
183 u[1]= x(fa7,1)-xg[1];
184 v[0]= -u[1];
185 v[1]= u[0];
186
187 int i0 = KEL_(0,fa7);
188 int i1 = KEL_(1,fa7);
189
190 // Normal orientation:
191 xj0[0]= x(i0,0)-x(i1,0);
192 xj0[1]= x(i0,1)-x(i1,1);
193
194 psc=xj0[0]*v[0] + xj0[1]*v[1];
195 if (psc < 0)
196 {
197 normales_facettes_Cl(num_elem,fa7,0) = -v[0];
198 normales_facettes_Cl(num_elem,fa7,1) = -v[1];
199 }
200 else
201 {
202 normales_facettes_Cl(num_elem,fa7,0) = v[0];
203 normales_facettes_Cl(num_elem,fa7,1) = v[1];
204 }
205
206}
207
208/*! @brief Modifies the interlaced volumes for face j of a non-standard element.
209 *
210 */
212 const Domaine_VEF& le_dom_VEF,
213 DoubleVect& volumes_entrelaces_Cl,
214 int type_cl) const
215{
216 double surf_mod;
217 const DoubleVect& volumes_entrelaces = le_dom_VEF.volumes_entrelaces();
218 const IntTab& elem_faces = le_dom_VEF.elem_faces();
219
220 switch(type_cl)
221 {
222
223 // no Dirichlet face: impossible
224 case 0:
225 {
226 Cerr << "Tri_VEF::modif_volumes_entrelaces() type 0 not possible!\n";
227 break;
228 }
229
230 case 1: // one Dirichlet face: Face 2
231 {
232 surf_mod = volumes_entrelaces[j]/2 ;
233 volumes_entrelaces_Cl[elem_faces(elem,0)] += surf_mod;
234 volumes_entrelaces_Cl[elem_faces(elem,1)] += surf_mod;
235 break;
236 }
237
238 case 2: // one Dirichlet face: Face 1
239 {
240 surf_mod = volumes_entrelaces[j]/2 ;
241 volumes_entrelaces_Cl[elem_faces(elem,0)] += surf_mod;
242 volumes_entrelaces_Cl[elem_faces(elem,2)] += surf_mod;
243 break;
244 }
245
246 case 4: // one Dirichlet face: Face 0
247 {
248 surf_mod = volumes_entrelaces[j]/2 ;
249 volumes_entrelaces_Cl[elem_faces(elem,1)] += surf_mod;
250 volumes_entrelaces_Cl[elem_faces(elem,2)] += surf_mod;
251 break;
252 }
253
254 case 6: // two Dirichlet faces: Faces 0,1
255 {
256 surf_mod = volumes_entrelaces[elem_faces(elem,0)]
257 + volumes_entrelaces[elem_faces(elem,1)];
258 volumes_entrelaces_Cl[elem_faces(elem,2)] += surf_mod;
259 break;
260 }
261
262 case 3: // two Dirichlet faces: Faces 1,2
263 {
264 surf_mod = volumes_entrelaces[elem_faces(elem,2)]
265 + volumes_entrelaces[elem_faces(elem,1)];
266 volumes_entrelaces_Cl[elem_faces(elem,0)] += surf_mod;
267 break;
268 }
269
270 case 5: // two Dirichlet faces: Faces 0,2
271 {
272 surf_mod = volumes_entrelaces[elem_faces(elem,0)]
273 + volumes_entrelaces[elem_faces(elem,2)];
274 volumes_entrelaces_Cl[elem_faces(elem,1)] += surf_mod;
275 break;
276 }
277
278 default :
279 {
280 Cerr << "\n unknown type in Tri_VEF::modif_volumes_entrelaces: " << type_cl ;
281 exit();
282 }
283
284 } // end switch
285
286}
287
288/*! @brief Modifies the interlaced volumes for joint face j of a non-standard element.
289 *
290 */
292 const Domaine_VEF& le_dom_VEF,
293 DoubleVect& volumes_entrelaces_Cl,
294 int type_cl) const
295{
296 double surf_mod;
297 const DoubleVect& volumes_entrelaces = le_dom_VEF.volumes_entrelaces();
298 const IntTab& elem_faces = le_dom_VEF.elem_faces();
299
300 int face;
301 int nb_faces_cl = volumes_entrelaces_Cl.size();
302 switch(type_cl)
303 {
304
305 // no Dirichlet face: impossible
306 case 0:
307 {
308 Cerr << "Tri_VEF::modif_volumes_entrelaces() type 0 not possible!\n";
309 break;
310 }
311
312 case 1: // one Dirichlet face: Face 2
313 {
314 surf_mod = volumes_entrelaces[j]/2 ;
315 face=elem_faces(elem,0);
316 if(face<nb_faces_cl)
317 volumes_entrelaces_Cl[face] += surf_mod;
318 face=elem_faces(elem,1);
319 if(face<nb_faces_cl)
320 volumes_entrelaces_Cl[face] += surf_mod;
321 break;
322 }
323
324 case 2: // one Dirichlet face: Face 1
325 {
326 surf_mod = volumes_entrelaces[j]/2 ;
327 face=elem_faces(elem,0);
328 if(face<nb_faces_cl)
329 volumes_entrelaces_Cl[face] += surf_mod;
330 face=elem_faces(elem,2);
331 if(face<nb_faces_cl)
332 volumes_entrelaces_Cl[face] += surf_mod;
333 break;
334 }
335
336 case 4: // one Dirichlet face: Face 0
337 {
338 surf_mod = volumes_entrelaces[j]/2 ;
339 face=elem_faces(elem,1);
340 if(face<nb_faces_cl)
341 volumes_entrelaces_Cl[face] += surf_mod;
342 face=elem_faces(elem,2);
343 if(face<nb_faces_cl)
344 volumes_entrelaces_Cl[face] += surf_mod;
345 break;
346 }
347
348 case 6: // two Dirichlet faces: Faces 0,1
349 {
350 surf_mod = volumes_entrelaces[elem_faces(elem,0)]
351 + volumes_entrelaces[elem_faces(elem,1)];
352 face=elem_faces(elem,2);
353 if(face<nb_faces_cl)
354 volumes_entrelaces_Cl[face] += surf_mod;
355 break;
356 }
357
358 case 3: // two Dirichlet faces: Faces 1,2
359 {
360 surf_mod = volumes_entrelaces[elem_faces(elem,2)]
361 + volumes_entrelaces[elem_faces(elem,1)];
362 face=elem_faces(elem,0);
363 if(face<nb_faces_cl)
364 volumes_entrelaces_Cl[face] += surf_mod;
365 break;
366 }
367
368 case 5: // two Dirichlet faces: Faces 0,2
369 {
370 surf_mod = volumes_entrelaces[elem_faces(elem,0)]
371 + volumes_entrelaces[elem_faces(elem,2)];
372 face=elem_faces(elem,1);
373 if(face<nb_faces_cl)
374 volumes_entrelaces_Cl[face] += surf_mod;
375 break;
376 }
377
378 default :
379 {
380 Cerr << "\n unknown type in Tri_VEF::modif_volumes_entrelaces: " << type_cl ;
381 exit();
382 }
383
384 } // end switch
385
386}
387
388/*! @brief
389 *
390 */
391void Tri_VEF::calcul_vc(const ArrOfInt& Face,ArrOfDouble& vc,
392 const ArrOfDouble& vs,const DoubleTab& vsom,
393 const Champ_Inc_base& vitesse,int type_cl, const DoubleVect& porosite_face) const
394{
395 switch(type_cl)
396 {
397 case 0: // triangle has no Dirichlet face
398 {
399 vc[0] = vs[0]/3;
400 vc[1] = vs[1]/3;
401 break;
402 }
403
404 case 1: // triangle has one Dirichlet face: Face 2
405 {
406 vc[0]= vitesse.valeurs()(Face[2],0)*porosite_face[Face[2]];
407 vc[1]= vitesse.valeurs()(Face[2],1)*porosite_face[Face[2]];
408 //vc[0]= vitesse.valeurs()(Face[2],0);
409 //vc[1]= vitesse.valeurs()(Face[2],1);
410 break;
411 }
412
413 case 2: // triangle has one Dirichlet face: Face 1
414 {
415 vc[0]= vitesse.valeurs()(Face[1],0)*porosite_face[Face[1]];
416 vc[1]= vitesse.valeurs()(Face[1],1)*porosite_face[Face[1]];
417 //vc[0]= vitesse.valeurs()(Face[1],0);
418 //vc[1]= vitesse.valeurs()(Face[1],1);
419 break;
420 }
421
422 case 4: // triangle has one Dirichlet face: Face 0
423 {
424 vc[0]= vitesse.valeurs()(Face[0],0)*porosite_face[Face[0]];
425 vc[1]= vitesse.valeurs()(Face[0],1)*porosite_face[Face[0]];
426 // vc[0]= vitesse.valeurs()(Face[0],0);
427 //vc[1]= vitesse.valeurs()(Face[0],1);
428 break;
429 }
430
431 case 3: // triangle has two Dirichlet faces: faces 1 and 2
432 {
433 vc[0]= vsom(0,0);
434 vc[1]= vsom(0,1);
435 break;
436 }
437
438 case 5: // triangle has two Dirichlet faces: faces 0 and 2
439 {
440 vc[0]= vsom(1,0);
441 vc[1]= vsom(1,1);
442 break;
443 }
444
445 case 6: // triangle has two Dirichlet faces: faces 0 and 1
446 {
447 vc[0]= vsom(2,0);
448 vc[1]= vsom(2,1);
449 break;
450 }
451
452 } // end switch
453
454}
455
456/*! @brief Computes the coordinates xg of the centre of a non-standard element.
457 *
458 * Also computes idirichlet = number of Dirichlet faces of the element.
459 * If idirichlet=2, n1 is the index of the node coinciding with G.
460 *
461 */
462void Tri_VEF::calcul_xg(DoubleVect& xg, const DoubleTab& x,
463 const int type_elem_Cl,int& idirichlet,int& n1,int& ,int& ) const
464{
465 int dim=xg.size();
466 switch(type_elem_Cl)
467 {
468
469 case 0: // triangle has no Dirichlet face: it has 3 facets
470 // point G is the barycentre of the triangle vertices
471 {
472 for (int j=0; j<dim; j++)
473 xg[j]=(x(0,j)+x(1,j)+x(2,j))/3;
474
475 idirichlet=0;
476 break;
477 }
478
479 case 1: // triangle has one Dirichlet face: Face 2
480 // point G is the barycentre of the vertices of Face 2
481
482 {
483 for (int j=0; j<dim; j++)
484 xg[j]=(x(0,j)+x(1,j))/2;
485
486 idirichlet=1;
487 break;
488 }
489
490 case 2: // triangle has one Dirichlet face: Face 1
491 // point G is the barycentre of the vertices of Face 1
492
493 {
494 for (int j=0; j<dim; j++)
495 xg[j]=(x(0,j)+x(2,j))/2;
496
497 idirichlet=1;
498 break;
499 }
500
501 case 4: // triangle has one Dirichlet face: Face 0
502 // point G is the barycentre of the vertices of Face 0
503
504 {
505 for (int j=0; j<dim; j++)
506 xg[j]=(x(1,j)+x(2,j))/2;
507
508 idirichlet=1;
509 break;
510 }
511
512 case 6 : // triangle has two Dirichlet faces: faces 0,1
513 // point G is vertex 2 of the triangle
514
515 {
516 for (int j=0; j<dim; j++)
517 xg[j]=x(2,j);
518
519 idirichlet=2;
520 n1 = 2;
521 break;
522
523 }
524
525 case 5 : // triangle has two Dirichlet faces: faces 0,2
526 // point G is vertex 1 of the triangle
527
528 {
529 for (int j=0; j<dim; j++)
530 xg[j]=x(1,j);
531
532 idirichlet=2;
533 n1 = 1;
534 break;
535
536 }
537
538 case 3 : // triangle has two Dirichlet faces: faces 1,2
539 // point G is vertex 0 of the triangle
540
541 {
542 for (int j=0; j<dim; j++)
543 xg[j]=x(0,j);
544
545 idirichlet=2;
546 n1 = 0;
547 break;
548
549 }
550
551 } // end switch
552
553}
554
555/*! @brief Modifies normales_facettes_Cl when idirichlet=2.
556 *
557 * idirichlet = number of Dirichlet faces of the element.
558 * If idirichlet=2, n1 is the index of the node coinciding with G.
559 *
560 */
561void Tri_VEF::modif_normales_facettes_Cl(DoubleTab& normales_facettes_Cl,
562 int fa7,int num_elem,
563 int idirichlet,int n1,int ,int ) const
564{
565 switch (idirichlet)
566 {
567
568 case 0:
569 break;
570
571 case 1:
572 break;
573
574 case 2:
575 {
576
577 // fa7=n1;
578 //normales_facettes_Cl(num_elem,fa7,0) = 0;
579 //normales_facettes_Cl(num_elem,fa7,1) = 0;
580 // the call seems unnecessary in 2D
581 normales_facettes_Cl(num_elem,n1,0) = 0;
582 normales_facettes_Cl(num_elem,n1,1) = 0;
583 break;
584 }
585 }
586}
Class Champ_Inc_base.
DoubleTab & valeurs() override
Returns the array of field values at the current time.
int_t nb_elem_tot() const
Definition Domaine.h:132
IntTab_t & les_elems()
Definition Domaine.h:129
const DoubleTab_t & coord_sommets() const
Definition Domaine.h:112
class Domaine_VEF
Definition Domaine_VEF.h:53
DoubleVect & volumes_entrelaces()
Definition Domaine_VF.h:99
int elem_faces(int i, int j) const
Returns the index of the i-th face of element num_elem; the face numbering convention is.
Definition Domaine_VF.h:542
const Domaine & domaine() const
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
const Nom & que_suis_je() const
Returns the string identifying the class.
Definition Objet_U.cpp:104
virtual Entree & readOn(Entree &)
Reads an Objet_U from an input stream. Virtual method to override.
Definition Objet_U.cpp:289
virtual Sortie & printOn(Sortie &) const
Writes the object to an output stream. Virtual method to override.
Definition Objet_U.cpp:278
static void exit(int exit_code=-1)
Exit routine for TRUST within a Kokkos region.
Definition Process.cpp:466
Base class for output streams.
Definition Sortie.h:52
_SIZE_ dimension_tot(int) const override
Definition TRUSTTab.tpp:160
_SIZE_ size() const
Definition TRUSTVect.tpp:45
void calcul_xg(DoubleVect &, const DoubleTab &, const int, int &, int &, int &, int &) const override
Computes the coordinates xg of the centre of a non-standard element.
Definition Tri_VEF.cpp:462
void creer_face_normales(DoubleTab &, const IntTab &, const IntTab &, const IntTab &, const Domaine &) const override
Fills the face_normales array in Domaine_VEF.
Definition Tri_VEF.cpp:56
void modif_normales_facettes_Cl(DoubleTab &, int, int, int, int, int, int) const override
Modifies normales_facettes_Cl when idirichlet=2.
Definition Tri_VEF.cpp:561
void creer_facette_normales(const Domaine_VEF &, const IntVect &) const override
Computes the facet normals for standard elements.
Definition Tri_VEF.cpp:104
void modif_volumes_entrelaces_faces_joints(int, int, const Domaine_VEF &, DoubleVect &, int) const override
Modifies the interlaced volumes for joint face j of a non-standard element.
Definition Tri_VEF.cpp:291
void calcul_vc(const ArrOfInt &, ArrOfDouble &, const ArrOfDouble &, const DoubleTab &, const Champ_Inc_base &, int, const DoubleVect &) const override
Definition Tri_VEF.cpp:391
Tri_VEF()
Definition Tri_VEF.cpp:36
void creer_normales_facettes_Cl(DoubleTab &, int, int, const DoubleTab &, const DoubleVect &, const Domaine &) const override
Fills the normales_facettes_Cl array in Domaine_Cl_VEF for facet fa7 of element num_elem.
Definition Tri_VEF.cpp:172
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 Tri_VEF.cpp:211