TrioCFD 1.9.8
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 et 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) sont les numeros locaux des 2 faces qui entourent la facette de numero local 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 remplit le tableau face_normales dans le 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 // Orientation de la normale de elem1 vers elem2
78 // pour cela recherche du sommet de elem1 qui n'est pas sur la 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 calcule les normales des facettes pour des elem standards
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 // la fa7 a pour sommets fa7 et "G" de coordonnees 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 // Orientation des normales :
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 remplit le tableau normales_facettes_Cl dans le Domaine_Cl_VEF pour la facette fa7 de l'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 // Orientation des normales :
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 modifie les volumes entrelaces pour la face j d'un elem non standard
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 // pas de Face de Dirichlet : impossible
224 case 0:
225 {
226 Cerr << "Tri_VEF::modif_volumes_entrelaces() type 0 impossible!\n";
227 break;
228 }
229
230 case 1: // une Face de Dirichlet : 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: // une Face de Dirichlet : 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: // une Face de Dirichlet : 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: // deux Faces de Dirichlet : 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: // deux Faces de Dirichlet : 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: // deux Faces de Dirichlet : 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 type inconnu Tri_VEF::modif_volumes_entrelaces: " << type_cl ;
281 exit();
282 }
283
284 } // fin du switch
285
286}
287
288/*! @brief modifie les volumes entrelaces pour la face joint j d'un elem non standard
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 // pas de Face de Dirichlet : impossible
306 case 0:
307 {
308 Cerr << "Tri_VEF::modif_volumes_entrelaces() type 0 impossible!\n";
309 break;
310 }
311
312 case 1: // une Face de Dirichlet : 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: // une Face de Dirichlet : 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: // une Face de Dirichlet : 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: // deux Faces de Dirichlet : 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: // deux Faces de Dirichlet : 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: // deux Faces de Dirichlet : 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 type inconnu Tri_VEF::modif_volumes_entrelaces: " << type_cl ;
381 exit();
382 }
383
384 } // fin du 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: // le triangle n'a pas de Face de Dirichlet
398 {
399 vc[0] = vs[0]/3;
400 vc[1] = vs[1]/3;
401 break;
402 }
403
404 case 1: // le triangle a une Face de Dirichlet :la 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: // le triangle a une Face de Dirichlet :la 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: // le triangle a une Face de Dirichlet :la 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: // le triangle a deux faces de Dirichlet :les faces 1 et 2
432 {
433 vc[0]= vsom(0,0);
434 vc[1]= vsom(0,1);
435 break;
436 }
437
438 case 5: // le triangle a deux faces de Dirichlet :les faces 0 et 2
439 {
440 vc[0]= vsom(1,0);
441 vc[1]= vsom(1,1);
442 break;
443 }
444
445 case 6: // le triangle a deux faces de Dirichlet :les faces 0 et 1
446 {
447 vc[0]= vsom(2,0);
448 vc[1]= vsom(2,1);
449 break;
450 }
451
452 } // fin du switch
453
454}
455
456/*! @brief calcule les coord xg du centre d'un element non standard calcule aussi idirichlet=nb de faces de Dirichlet de l'element
457 *
458 * si idirichlet=2, n1 est le numero du sommet confondu avec G
459 *
460 */
461void Tri_VEF::calcul_xg(DoubleVect& xg, const DoubleTab& x,
462 const int type_elem_Cl,int& idirichlet,int& n1,int& ,int& ) const
463{
464 int dim=xg.size();
465 switch(type_elem_Cl)
466 {
467
468 case 0: // le triangle n'a pas de Face de dirichlet: il a 3 Facettes
469 // le point G est le barycentre des sommets du triangle
470 {
471 for (int j=0; j<dim; j++)
472 xg[j]=(x(0,j)+x(1,j)+x(2,j))/3;
473
474 idirichlet=0;
475 break;
476 }
477
478 case 1: // le triangle a une Face de dirichlet: la Face 2
479 // le point G est le barycentre des sommets de la Face 2
480
481 {
482 for (int j=0; j<dim; j++)
483 xg[j]=(x(0,j)+x(1,j))/2;
484
485 idirichlet=1;
486 break;
487 }
488
489 case 2: // le triangle a une Face de dirichlet: la Face 1
490 // le point G est le barycentre des sommets de la Face 1
491
492 {
493 for (int j=0; j<dim; j++)
494 xg[j]=(x(0,j)+x(2,j))/2;
495
496 idirichlet=1;
497 break;
498 }
499
500 case 4: // le triangle a une Face de dirichlet: la Face 0
501 // le point G est le barycentre des sommets de la Face 0
502
503 {
504 for (int j=0; j<dim; j++)
505 xg[j]=(x(1,j)+x(2,j))/2;
506
507 idirichlet=1;
508 break;
509 }
510
511 case 6 : // le triangle a deux faces de Dirichlet : les faces 0,1
512 // le point G est le sommet 2 du triangle
513
514 {
515 for (int j=0; j<dim; j++)
516 xg[j]=x(2,j);
517
518 idirichlet=2;
519 n1 = 2;
520 break;
521
522 }
523
524 case 5 : // le triangle a deux faces de Dirichlet : les faces 0,2
525 // le point G est le sommet 1 du triangle
526
527 {
528 for (int j=0; j<dim; j++)
529 xg[j]=x(1,j);
530
531 idirichlet=2;
532 n1 = 1;
533 break;
534
535 }
536
537 case 3 : // le triangle a deux faces de Dirichlet : les faces 1,2
538 // le point G est le sommet 0 du triangle
539
540 {
541 for (int j=0; j<dim; j++)
542 xg[j]=x(0,j);
543
544 idirichlet=2;
545 n1 = 0;
546 break;
547
548 }
549
550 } // fin du switch
551
552}
553
554/*! @brief modifie normales_facettes_Cl quand idirichlet=2 idirichlet=nb de faces de Dirichlet de l'element
555 *
556 * si idirichlet=2, n1 est le numero du sommet confondu avec G
557 *
558 */
559void Tri_VEF::modif_normales_facettes_Cl(DoubleTab& normales_facettes_Cl,
560 int fa7,int num_elem,
561 int idirichlet,int n1,int ,int ) const
562{
563 switch (idirichlet)
564 {
565
566 case 0:
567 break;
568
569 case 1:
570 break;
571
572 case 2:
573 {
574
575 // fa7=n1;
576 //normales_facettes_Cl(num_elem,fa7,0) = 0;
577 //normales_facettes_Cl(num_elem,fa7,1) = 0;
578 // l'appel semble inutile en 2D
579 normales_facettes_Cl(num_elem,n1,0) = 0;
580 normales_facettes_Cl(num_elem,n1,1) = 0;
581 break;
582 }
583 }
584}
Classe Champ_Inc_base.
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
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:54
DoubleVect & volumes_entrelaces()
Definition Domaine_VF.h:99
int elem_faces(int i, int j) const
renvoie le numero de le ieme face de la maille num_elem la facon dont ces faces sont numerotees est
Definition Domaine_VF.h:543
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
renvoie la chaine identifiant la classe.
Definition Objet_U.cpp:104
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
Definition Objet_U.cpp:293
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
Classe de base des flux de sortie.
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
calcule les coord xg du centre d'un element non standard calcule aussi idirichlet=nb de faces de Diri...
Definition Tri_VEF.cpp:461
void creer_face_normales(DoubleTab &, const IntTab &, const IntTab &, const IntTab &, const Domaine &) const override
remplit le tableau face_normales dans le Domaine_VEF
Definition Tri_VEF.cpp:56
void modif_normales_facettes_Cl(DoubleTab &, int, int, int, int, int, int) const override
modifie normales_facettes_Cl quand idirichlet=2 idirichlet=nb de faces de Dirichlet de l'element
Definition Tri_VEF.cpp:559
void creer_facette_normales(const Domaine_VEF &, const IntVect &) const override
calcule les normales des facettes pour des elem standards
Definition Tri_VEF.cpp:104
void modif_volumes_entrelaces_faces_joints(int, int, const Domaine_VEF &, DoubleVect &, int) const override
modifie les volumes entrelaces pour la face joint j d'un elem non standard
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
remplit le tableau normales_facettes_Cl dans le Domaine_Cl_VEF pour la facette fa7 de l'element num_e...
Definition Tri_VEF.cpp:172
void modif_volumes_entrelaces(int, int, const Domaine_VEF &, DoubleVect &, int) const override
modifie les volumes entrelaces pour la face j d'un elem non standard
Definition Tri_VEF.cpp:211