TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Op_Conv_kschemas_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 <Op_Conv_kschemas_VEF.h>
17#include <Periodique.h>
18#include <Neumann_sortie_libre.h>
19
20Implemente_base(Op_Conv_kschemas_VEF,"Op_Conv_kschemas_VEF_P1NC",Op_Conv_VEF_base);
21
22
23Sortie& Op_Conv_kschemas_VEF::printOn(Sortie& s ) const
24{
25 return s << que_suis_je() ;
26}
27
28Entree& Op_Conv_kschemas_VEF::readOn(Entree& s )
29{
30 return s ;
31}
32
33void Op_Conv_kschemas_VEF::associer(const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_cl_dis, const Champ_Inc_base& ch)
34{
35 // CCa le 28/05/99 Le schema Kquick ne marche pas en paralle !!
36 if (Process::is_parallel())
37 {
38 Cerr << "ATTENTION le kquick ne marche pas en parallele !!!" << finl;
39 exit();
40 }
41
42 Op_Conv_VEF_base::associer(domaine_dis, domaine_cl_dis, ch);
43}
44
45
46//////////////////////////////////////////////////////////////
47// Fonctions pour les kschemas.
48////////////////////////////////////////////////////////////////
49
50// convkschemas : fonction utilitaire pour la convection
51
52void convkschemas(const double K, const int ncomp, int dimension, const int poly ,
53 const int poly1, const int poly2,const int jel0,
54 const int jel1,const double psc ,const DoubleTab& tab1 ,
55 DoubleVect& fluent, DoubleVect& flux,
56 const DoubleVect& rx0, const DoubleTab& gradient_elem )
57{
58
59 int comp,amont,i,elem1,elem2;
60 double CF,UTC,deltat0,deltat1,deltat;
61 DoubleVect rx(dimension);
62 deltat0 = 0.;
63 deltat1 = 0.;
64
65 ////////////////////////////////////////////////////////////////////////
66 // Test sur les bords
67 ////////////////////////////////////////////////////////////////////////
68 if ((poly1==-1) || (poly2==-1))
69 {
70 if (psc >= 0)
71 {
72 amont = jel0;
73 fluent[jel1] += psc;
74 }
75 else
76 {
77 amont = jel1;
78 fluent[jel0] -= psc;
79 }
80
81 for (comp=0; comp<ncomp; comp++)
82 flux(comp) = tab1(amont,comp);
83
84 }
85 else
86 {
87 if (psc >= 0)
88 {
89 amont = jel0;
90 rx = rx0;
91 elem1 = poly;
92 elem2 = poly1;
93 fluent(jel1) += psc;
94 }
95 else
96 {
97 amont = jel1;
98 rx = rx0;
99 rx *= -1.;
100 elem1 = poly;
101 elem2 = poly2;
102 fluent(jel0) -= psc;
103 }
104
105 for (comp=0; comp<ncomp; comp++)
106 {
107 deltat0 = deltat1 = 0.0;
108 flux(comp) = tab1(amont,comp);
109 //Cerr << " flux(" << comp << ") ie phiamont= " << flux(comp) << finl;
110
111 for (i=0; i<dimension; i++)
112 {
113 deltat0 += gradient_elem(elem1,comp,i)*rx(i);
114 deltat1 += gradient_elem(elem2,comp,i)*rx(i);
115 }
116
117 if (K == 0.5)
118 {
119 deltat = deltat0 + deltat1;
120
121 if (std::fabs(deltat) <= 1.e-5)
122 {
123 CF = 0.125;
124 }
125 else
126 {
127 UTC = deltat1 / deltat;
128
129 if ( (UTC <= -1.) || (UTC >= 1.5) ) CF = 0.125;
130 else if ((UTC > -1.) && (UTC <= 0.)) CF = 0.5 + 0.375*UTC;
131 else if ((UTC > 0.) && (UTC <= 0.25)) CF = 0.5 - 0.625*sqrt(UTC);
132 else if ((UTC > 0.25) && (UTC < 1.5 )) CF = 0.25* std::fabs(UTC - 1.);
133 else
134 {
135 CF=0.;
136 Process::exit();
137 }
138 }
139 // Calcul du flux
140 flux(comp) += (0.5 - CF)*deltat0 + CF*deltat1 ;
141 //Cerr << " flux(" << comp << ")= " << flux(comp) << finl;
142 }
143 else
144 {
145 // Calcul du flux
146 flux(comp) += 0.25*((1.+K)*deltat0 + (1.-K)*deltat1) ;
147 }
148 }
149
150 }
151}
152
153//////////////////////////////////////////////////////////////////////////
154// Procedure AJOUTER
155/////////////////////////////////////////////////////////////////////////
156
157DoubleTab& Op_Conv_kschemas_VEF::ajouter(const DoubleTab& transporte,
158 DoubleTab& resu) const
159{
160 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
161 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
162 const DoubleVect& porosite_face = equation().milieu().porosite_face();
163 const Champ_Inc_base& la_vitesse=vitesse_.valeur();
164 const IntTab& elem_faces = domaine_VEF.elem_faces();
165 const DoubleTab& face_normales = domaine_VEF.face_normales();
166 const auto& facette_normales = domaine_VEF.facette_normales();
167 const Domaine& domaine = domaine_VEF.domaine();
168 const int nb_faces = domaine_VEF.nb_faces();
169 const int nfa7 = domaine_VEF.type_elem().nb_facette();
170 const int nb_elem = domaine_VEF.nb_elem();
171 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
172 const IntVect& rang_elem_non_std = domaine_VEF.rang_elem_non_std();
173 const IntTab& face_voisins = domaine_VEF.face_voisins();
174 const DoubleVect& volumes = domaine_VEF.volumes();
175 const DoubleTab& xv = domaine_VEF.xv();
176 const DoubleTab& xg = domaine_VEF.xp();
177 const DoubleTab& coord = domaine.coord_sommets();
178 int premiere_face_int = domaine_VEF.premiere_face_int();
179 const IntTab& les_Polys = domaine.les_elems();
180
181 const DoubleTab& normales_facettes_Cl = domaine_Cl_VEF.normales_facettes_Cl();
182
183 int nfac = domaine.nb_faces_elem();
184 int nsom = domaine.nb_som_elem();
185 int nb_som_facette = domaine.type_elem()->nb_som_face();
186
187 // Pour le traitement de la convection on distingue les polyedres
188 // standard qui ne "voient" pas les conditions aux limites et les
189 // polyedres non standard qui ont au moins une face sur le bord.
190 // Un polyedre standard a n facettes sur lesquelles on applique le
191 // schema de convection.
192 // Pour un polyedre non standard qui porte des conditions aux limites
193 // de Dirichlet, une partie des facettes sont portees par les faces.
194 // En bref pour un polyedre le traitement de la convection depend
195 // du type (triangle, tetraedre ...) et du nombre de faces de Dirichlet.
196
197 double psc;
198 int poly,poly1,poly2,face_adj,fa7,i,j,n_bord;
199 int num_face, rang ,itypcl;
200 int num10,num20,num3,num_som;
201
202 // MODIF SB su 10/09/03
203 // Pour les 3 elements suivants, il y a autant de sommets que de face
204 // constituant l'element geometrique
205 // PB avec les hexa, 8 sommets et 6 faces, donc l'utilisation du tableau
206 // face[i] ne fonctionne plus
207 // la methode retenue pour eviter de calculer la vitesse aux sommets sans
208 // les fonctions de forme n'est donc pas utilisable,
209 // pour l'hexa on n'a pas acces a la face.
210 // il existe le tableau Face=>sommets mais pas l'inverse.
211 // trop couteux et pour le moment on n'etend pas les porosites aux hexa
212
213 int istetra=0;
214 const Elem_VEF_base& type_elemvef= domaine_VEF.type_elem();
215 Nom nom_elem=type_elemvef.que_suis_je();
216 if ((nom_elem=="Tetra_VEF")||(nom_elem=="Tri_VEF")) istetra=1;
217
218 const int ncomp_ch_transporte= transporte.line_size();
219 int fac,elem1,elem2,comp0;
220 int nb_faces_ = domaine_VEF.nb_faces();
221 IntVect face(nfac);
222
223 DoubleVect flux(ncomp_ch_transporte);
224 DoubleVect fluxsom(ncomp_ch_transporte);
225 DoubleVect fluxg(ncomp_ch_transporte);
226
227 // Traitement particulier pour les faces de periodicite
228 int nb_faces_perio = 0;
229 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
230 {
231 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
232 if (sub_type(Periodique,la_cl.valeur()))
233 {
234 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
235 int num1 = le_bord.num_premiere_face();
236 int num2 = num1 + le_bord.nb_faces();
237 for (num_face=num1; num_face<num2; num_face++)
238 nb_faces_perio++;
239 }
240 }
241
242 DoubleTab tab(nb_faces_perio,ncomp_ch_transporte);
243
244 nb_faces_perio=0;
245 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
246 {
247 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
248 if (sub_type(Periodique,la_cl.valeur()))
249 {
250 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
251 int num1 = le_bord.num_premiere_face();
252 int num2 = num1 + le_bord.nb_faces();
253 for (num_face=num1; num_face<num2; num_face++)
254 {
255 for (int comp=0; comp<ncomp_ch_transporte; comp++)
256 tab(nb_faces_perio,comp) = resu(num_face,comp);
257 nb_faces_perio++;
258 }
259 }
260 }
261
262
263 ///////////////////////////////////////////////////////////////////////////////////////////////
264 // <
265 // calcul des gradients; < [ Ujp*np/vol(j) ]
266 // j
267 ////////////////////////////////////////////////////////////////////////////////////////////////
268 DoubleTab gradient_elem(0, ncomp_ch_transporte, dimension);
269 domaine_VEF.domaine().creer_tableau_elements(gradient_elem);
270 // Boucle sur les faces
271
272
273 for (fac=0; fac< premiere_face_int; fac++)
274 {
275 elem1=face_voisins(fac,0);
276 if(ncomp_ch_transporte==1)
277 for (i=0; i<dimension; i++)
278 {
279 gradient_elem(elem1, 0, i) +=
280 face_normales(fac,i)*transporte(fac);
281 }
282 else
283 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
284 for (i=0; i<dimension; i++)
285 gradient_elem(elem1, comp0, i) +=
286 face_normales(fac,i)*transporte(fac,comp0);
287 // dUcomp/dXi
288 } // fin du for faces
289
290 for (; fac<nb_faces_; fac++)
291 {
292 elem1=face_voisins(fac,0);
293 elem2=face_voisins(fac,1);
294 if(ncomp_ch_transporte==1)
295 for (i=0; i<dimension; i++)
296 {
297 gradient_elem(elem1, 0, i) +=
298 face_normales(fac,i)*transporte(fac);
299 gradient_elem(elem2, 0, i) -=
300 face_normales(fac,i)*transporte(fac);
301 }
302 else
303 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
304 for (i=0; i<dimension; i++)
305 {
306 gradient_elem(elem1, comp0, i) +=
307 face_normales(fac,i)*transporte(fac,comp0);
308 gradient_elem(elem2, comp0, i) -=
309 face_normales(fac,i)*transporte(fac,comp0);
310 }
311 // dUcomp/dXi
312 } // fin du for faces
313
314 for (int elem=0; elem<nb_elem; elem++)
315 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
316 for (i=0; i<dimension; i++)
317 gradient_elem(elem,comp0,i) /= volumes(elem);
318
319 gradient_elem.echange_espace_virtuel();
320
321 ////////////////////////////////////////////////////////////////////////////////////
322 // On a les gradient_elem par elements
323 ////////////////////////////////////////////////////////////////////////////////////
324
325 DoubleVect vs(dimension);
326 DoubleVect vc(dimension);
327 DoubleTab vsom(nsom,dimension);
328 DoubleVect cc(dimension);
329 double xm;
330
331 // On remet a zero le tableau qui sert pour
332 // le calcul du pas de temps de stabilite
333 fluent_ = 0;
334
335 // Les polyedres non standard sont ranges en 2 groupes dans le Domaine_VEF:
336 // - polyedres bords et joints
337 // - polyedres bords et non joints
338 // On traite les polyedres en suivant l'ordre dans lequel ils figurent
339 // dans le domaine
340
341 //////////////////////////////////////////////////////////////////////////////////////
342 // boucle sur les polys
343 //////////////////////////////////////////////////////////////////////////////////////
344 const IntTab& KEL=domaine_VEF.type_elem().KEL();
345 for (poly=0; poly<nb_elem; poly++)
346 {
347 rang = rang_elem_non_std(poly);
348 if (rang==-1)
349 itypcl=0;
350 else
351 itypcl=domaine_Cl_VEF.type_elem_Cl(rang);
352
353 // calcul des numeros des faces du polyedre
354 for (face_adj=0; face_adj<nfac; face_adj++)
355 face(face_adj)= elem_faces(poly,face_adj);
356
357 int scom;
358 DoubleVect rx0(dimension);
359
360 // calcul de la vitesse aux sommets des polyedres
361 for (j=0; j<dimension; j++)
362 {
363 vs(j) = la_vitesse.valeurs()(face(0),j)*porosite_face(face(0));
364 for (i=1; i<nfac; i++)
365 vs(j)+= la_vitesse.valeurs()(face(i),j)*porosite_face(face(i));
366 }
367
368 //int ncomp;
369 if (istetra==1)
370 {
371 for (j=0; j<nsom; j++)
372 {
373 for (int ncomp=0; ncomp<Objet_U::dimension; ncomp++)
374 vsom(j,ncomp) =vs[ncomp] - Objet_U::dimension*la_vitesse.valeurs()(face[j],ncomp)*porosite_face(face[j]);
375 }
376 }
377 else
378 {
379 // pour que cela soit valide avec les hexa
380 // On va utliser les fonctions de forme implementees dans la classe Champs_P1_impl ou Champs_Q1_impl
381 //int ncomp;
382 for (j=0; j<nsom; j++)
383 {
384 num_som = domaine.sommet_elem(poly,j);
385 for (int ncomp=0; ncomp<dimension; ncomp++)
386 {
387 vsom(j,ncomp) = la_vitesse.valeur_a_sommet_compo(num_som,poly,ncomp);
388 }
389 }
390 }
391 // calcul de la vitesse au centre de gravite
392 domaine_VEF.type_elem().calcul_vc(face,vc,vs,vsom,vitesse(),itypcl,porosite_face);
393
394 // Boucle sur les facettes du polyedre non standard:
395 for (fa7=0; fa7<nfa7; fa7++)
396 {
397 //Cerr << "la facette etudiee est " << fa7 << finl;
398 // fa7 separe num1 et num2. num3 est la troisieme face (2D).
399
400 num10 = face(KEL(0,fa7));
401 num20 = face(KEL(1,fa7));
402 num3 = face(KEL(2,fa7));
403
404 // Determination des elements voisins aux faces num1 et num2
405
406 poly1 = face_voisins(num10,0);
407 if (poly1==poly)
408 poly1 = face_voisins(num10,1);
409
410 poly2 = face_voisins(num20,0);
411 if (poly2==poly)
412 poly2 = face_voisins(num20,1);
413
414 scom = les_Polys(poly,KEL(2,fa7));
415
416 // calcul des rx0, distance entre les milieux des 'num i'
417
418 for (i=0; i<dimension; i++)
419 rx0(i) = xv(num20,i)-xv(num10,i);
420
421 // normales aux facettes
422
423 if (rang==-1)
424 for (i=0; i<dimension; i++)
425 cc[i] = facette_normales(poly, fa7, i);
426 else
427 for (i=0; i<dimension; i++)
428 cc[i] = normales_facettes_Cl(rang,fa7,i);
429
430 /////////////////////////////////////////////////////////////////////////
431 // On traite le point pour lequel vitesse = 0.5(vitsommet + vitmilieu)
432 /////////////////////////////////////////////////////////////////////////
433
434 for (i=0; i<nb_som_facette-1; i++)
435 {
436 //////////////////////////////////////////////////////////////////////////
437 //Determination de PhiIJ au milieu entre le sommet et le milieu de num3
438 /////////////////////////////////////////////////////////////////////////
439
440 psc = 0;
441 for (j=0; j<dimension; j++)
442 psc+=((vsom(KEL(i+2,fa7),j) + la_vitesse.valeurs()(num3,j) * porosite_face(num3)))*cc[j];
443 psc *=0.5;
444 convkschemas(K,ncomp_ch_transporte,dimension,poly,poly1,poly2,num10,num20,psc,transporte,
445 fluent_,flux,rx0,gradient_elem);
446
447 ////////////////////////////////////////////////////////////////////////////////////////////////////////
448 //Limiteur pour le gradient. Calcul effectue en meme temps que le calcul du flux.
449 // gradient(K0) = teta*gradient(K0)+ (1-teta)gradient(K1ou K2)
450 ////////////////////////////////////////////////////////////////////////////////////////////////////////
451
452 double teta = 0.5;
453
454 /////////////////////////////////////////////////////////////////////////
455 // On traite les sommets qui sont aussi des sommets du polyedre
456 /////////////////////////////////////////////////////////////////////////
457
458 // XXX XXX XXX : Attention : we can not factorize more... the code is not the same
459 if (ncomp_ch_transporte == 1)
460 {
461 for (j=0; j<dimension; j++)
462 {
463 xm = 0.5 *(coord(scom,j)+xv(num3,j));
464 if (psc >= 0)
465 {
466 if (poly1==-1)
467 fluxsom(0) += gradient_elem(poly,0,j)*(coord(scom,j)-xm);
468 else
469 fluxsom(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly1,0,j))*(coord(scom,j)-xm);
470 }
471 else
472 {
473 if (poly2==-1)
474 fluxsom(0) += gradient_elem(poly,0,j)*(coord(scom,j)-xm);
475 else
476 fluxsom(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly2,0,j))*(coord(scom,j)-xm);
477 }
478 }
479 fluxsom(0) += flux(0);
480 fluxsom(0) *= psc;
481 }
482 else
483 {
484 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
485 {
486 fluxsom(comp0) = flux(comp0);
487 for (j=0; j<dimension; j++)
488 {
489 xm = 0.5 *(coord(scom,j)+xv(num3,j));
490 if (psc >= 0)
491 {
492 if (poly1==-1)
493 fluxsom(comp0) += gradient_elem(poly,comp0,j)*(coord(scom,j)-xm);
494 else
495 fluxsom(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly1,comp0,j))*(coord(scom,j)-xm);
496 }
497 else
498 {
499 if (poly2==-1)
500 fluxsom(comp0) += gradient_elem(poly,comp0,j)*(coord(scom,j)-xm);
501 else
502 fluxsom(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly2,comp0,j))*(coord(scom,j)-xm);
503 }
504 }
505 fluxsom(comp0) *= psc;
506 }
507 }
508
509 ////////////////////////////////////////////////////////////////////////////
510 // on traite le centre de gravite
511 ////////////////////////////////////////////////////////////////////////////
512
513 // XXX XXX XXX : Attention : we can not factorize more... the code is not the same
514 if (ncomp_ch_transporte == 1)
515 {
516 for (j=0; j<dimension; j++)
517 {
518 xm = 0.5 *(coord(scom,j)+xv(num3,j));
519 if (psc >= 0)
520 {
521 if (poly1==-1)
522 fluxg(0) += gradient_elem(poly,0,j)*(xg(poly,j)-xm);
523 else
524 fluxg(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly1,0,j))*(xg(poly,j)-xm);
525 }
526
527 else
528 {
529 if (poly2==-1)
530 fluxg(0) += gradient_elem(poly,0,j)*(xg(poly,j)-xm);
531 else
532 fluxg(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly2,0,j))*(xg(poly,j)-xm);
533 }
534 }
535
536 fluxg(0) += flux(0);
537 fluxg(0) *= psc;
538 }
539 else
540 {
541 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
542 {
543 fluxg(comp0) = flux(comp0);
544 for (j=0; j<dimension; j++)
545 {
546 xm = 0.5 *(coord(scom,j)+xv(num3,j));
547 if (psc >= 0)
548 {
549 if (poly1==-1)
550 fluxg(comp0) += gradient_elem(poly,comp0,j)*(xg(poly,j)-xm);
551 else
552 fluxg(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly1,comp0,j))*(xg(poly,j)-xm);
553 }
554
555 else
556 {
557 if (poly2==-1)
558 fluxg(comp0) += gradient_elem(poly,comp0,j)*(xg(poly,j)-xm);
559 else
560 fluxg(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly2,comp0,j))*(xg(poly,j)-xm);
561 }
562 }
563 fluxg(comp0) *= psc;
564 }
565 }
566 //////////////////////////////////////////////////////////////////////////////
567 // Integration de u.n.flux
568 /////////////////////////////////////////////////////////////////////////////
569 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
570 {
571 resu(num10,comp0) -= ( 0.5*(fluxsom(comp0)+fluxg(comp0)) );
572 resu(num20,comp0) += ( 0.5*(fluxsom(comp0)+fluxg(comp0)) );
573 }
574 }
575 }
576 } // fin de la boucle
577
578 // Traitement des elements joints d'epaisseur 1
579 for (poly=0; poly<nb_elem_tot; poly++)
580 {
581 // On regarde si une face du polyedre est une face joint
582 for (face_adj=0; face_adj<nfac; face_adj++)
583 if(face_adj<nb_faces) break;
584 if(face_adj<nfac)
585 {
586 rang = rang_elem_non_std(poly);
587 if (rang==-1)
588 itypcl=0;
589 else
590 itypcl=domaine_Cl_VEF.type_elem_Cl(rang);
591
592 // calcul des numeros des faces du polyedre
593 for (face_adj=0; face_adj<nfac; face_adj++)
594 {
595 face(face_adj)= elem_faces(poly,face_adj);
596 //Cerr << "les faces de l'elements sont : " << face(face_adj) << finl;
597 }
598
599 int scom;
600 DoubleVect rx0(dimension);
601
602 // calcul de la vitesse aux sommets des polyedres
603 for (j=0; j<dimension; j++)
604 {
605 vs(j) = la_vitesse.valeurs()(face(0),j)*porosite_face(face(0));
606 for (i=1; i<nfac; i++)
607 vs(j)+= la_vitesse.valeurs()(face(i),j)*porosite_face(face(j));
608 }
609 int ncomp;
610 for (j=0; j<nsom; j++)
611 {
612 num_som = domaine.sommet_elem(poly,j);
613 for (ncomp=0; ncomp<dimension; ncomp++)
614 vsom(j,ncomp) = la_vitesse.valeur_a_sommet_compo(num_som,poly,ncomp);
615 }
616 // calcul de la vitesse au centre de gravite
617
618 domaine_VEF.type_elem().calcul_vc(face,vc,vs,vsom,vitesse(),itypcl,porosite_face);
619
620
621 // Boucle sur les facettes du polyedre non standard:
622
623 for (fa7=0; fa7<nfa7; fa7++)
624 {
625 //Cerr << "la facette etudiee est " << fa7 << finl;
626 // fa7 separe num1 et num2. num3 est la troisieme face (2D).
627
628 num10 = face(KEL(0,fa7));
629 num20 = face(KEL(1,fa7));
630 num3 = face(KEL(2,fa7));
631
632 // Determination des elements voisins aux faces num1 et num2
633
634 poly1 = face_voisins(num10,0);
635 if (poly1==poly)
636 {
637 poly1 = face_voisins(num10,1);
638 }
639
640 poly2 = face_voisins(num20,0);
641 if (poly2==poly)
642 {
643 poly2 = face_voisins(num20,1);
644 }
645
646 scom = les_Polys(poly,KEL(2,fa7));
647
648 // calcul des rx0, distance entre les milieux des 'num i'
649
650 for (i=0; i<dimension; i++)
651 rx0(i) = xv(num20,i)-xv(num10,i);
652
653 // normales aux facettes
654
655 if (rang==-1)
656 for (i=0; i<dimension; i++)
657 cc[i] = facette_normales(poly, fa7, i);
658 else
659 for (i=0; i<dimension; i++)
660 cc[i] = normales_facettes_Cl(rang,fa7,i);
661
662 /////////////////////////////////////////////////////////////////////////
663 // On traite le point pour lequel vitesse = 0.5(vitsommet + vitmilieu)
664 /////////////////////////////////////////////////////////////////////////
665
666 for (i=0; i<nb_som_facette-1; i++)
667 {
668 //////////////////////////////////////////////////////////////////////////
669 //Determination de PhiIJ au milieu entre le sommet et le milieu de num3
670 /////////////////////////////////////////////////////////////////////////
671
672 psc = 0;
673 for (j=0; j<dimension; j++)
674 psc+=((vsom(KEL(i+2,fa7),j) + la_vitesse.valeurs()(num3,j) * porosite_face(num3)))*cc[j];
675 psc *=0.5;
676 convkschemas(K,ncomp_ch_transporte,dimension,poly,poly1,poly2,num10,num20,psc,transporte,
677 fluent_,flux,rx0,gradient_elem);
678
679
680 ////////////////////////////////////////////////////////////////////////////////////////////////////////
681 //Limiteur pour le gradient. Calcul effectue en meme temps que le calcul du flux.
682 // gradient(K0) = teta*gradient(K0)+ (1-teta)gradient(K1ou K2)
683 ////////////////////////////////////////////////////////////////////////////////////////////////////////
684
685 double teta = 0.5;
686
687 /////////////////////////////////////////////////////////////////////////
688 // On traite les sommets qui sont aussi des sommets du polyedre
689 /////////////////////////////////////////////////////////////////////////
690
691 // XXX XXX XXX : Attention : we can not factorize more... the code is not the same
692 if (ncomp_ch_transporte == 1)
693 {
694 for (j=0; j<dimension; j++)
695 {
696 xm = 0.5 *(coord(scom,j)+xv(num3,j));
697 if (psc >= 0)
698 {
699 if (poly1==-1)
700 fluxsom(0) += gradient_elem(poly,0,j)*(coord(scom,j)-xm);
701 else
702 fluxsom(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly1,0,j))*(coord(scom,j)-xm);
703 }
704 else
705 {
706 if (poly2==-1)
707 fluxsom(0) += gradient_elem(poly,0,j)*(coord(scom,j)-xm);
708 else
709 fluxsom(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly2,0,j))*(coord(scom,j)-xm);
710 }
711 }
712 fluxsom(0) += flux(0);
713 fluxsom(0) *= psc;
714 }
715 else
716 {
717 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
718 {
719 fluxsom(comp0) = flux(comp0);
720 for (j=0; j<dimension; j++)
721 {
722 xm = 0.5 *(coord(scom,j)+xv(num3,j));
723 if (psc >= 0)
724 {
725 if (poly1==-1)
726 fluxsom(comp0) += gradient_elem(poly,comp0,j)*(coord(scom,j)-xm);
727 else
728 fluxsom(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly1,comp0,j))*(coord(scom,j)-xm);
729 }
730 else
731 {
732 if (poly2==-1)
733 fluxsom(comp0) += gradient_elem(poly,comp0,j)*(coord(scom,j)-xm);
734 else
735 fluxsom(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly2,comp0,j))*(coord(scom,j)-xm);
736 }
737 }
738 fluxsom(comp0) *= psc;
739 }
740 }
741
742 ////////////////////////////////////////////////////////////////////////////
743 // on traite le centre de gravite
744 ////////////////////////////////////////////////////////////////////////////
745
746 // XXX XXX XXX : Attention : we can not factorize more... the code is not the same
747 if (ncomp_ch_transporte == 1)
748 {
749 for (j=0; j<dimension; j++)
750 {
751 xm = 0.5 *(coord(scom,j)+xv(num3,j));
752 if (psc >= 0)
753 {
754 if (poly1==-1)
755 fluxg(0) += gradient_elem(poly,0,j)*(xg(poly,j)-xm);
756 else
757 fluxg(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly1,0,j))*(xg(poly,j)-xm);
758 }
759 else
760 {
761 if (poly2==-1)
762 fluxg(0) += gradient_elem(poly,0,j)*(xg(poly,j)-xm);
763 else
764 fluxg(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly2,0,j))*(xg(poly,j)-xm);
765 }
766 }
767 fluxg(0) += flux(0);
768 fluxg(0) *= psc;
769 }
770 else
771 {
772 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
773 {
774 fluxg(comp0) = flux(comp0);
775 for (j=0; j<dimension; j++)
776 {
777 xm = 0.5 *(coord(scom,j)+xv(num3,j));
778 if (psc >= 0)
779 {
780 if (poly1==-1)
781 fluxg(comp0) += gradient_elem(poly,comp0,j)*(xg(poly,j)-xm);
782 else
783 fluxg(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly1,comp0,j))*(xg(poly,j)-xm);
784 }
785 else
786 {
787 if (poly2==-1)
788 fluxg(comp0) += gradient_elem(poly,comp0,j)*(xg(poly,j)-xm);
789 else
790 fluxg(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly2,comp0,j))*(xg(poly,j)-xm);
791 }
792 }
793 fluxg(comp0) *= psc;
794 }
795 }
796
797 //////////////////////////////////////////////////////////////////////////////
798 // Integration de u.n.flux
799 /////////////////////////////////////////////////////////////////////////////
800 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
801 {
802 resu(num10,comp0) -= ( 0.5*(fluxsom(comp0)+fluxg(comp0)) );
803 resu(num20,comp0) += ( 0.5*(fluxsom(comp0)+fluxg(comp0)) );
804 }
805 }
806 }
807 }
808 } // fin de la boucle
809 int voisine;
810 nb_faces_perio = 0;
811 double diff1,diff2;
812
813 // Dimensionnement du tableau des flux convectifs au bord du domaine
814 // de calcul
815 DoubleTab& flux_b = flux_bords_;
816 flux_b.resize(domaine_VEF.nb_faces_bord(),ncomp_ch_transporte);
817 flux_b = 0.;
818
819 // Boucle sur les bords pour traiter les conditions aux limites
820 // il y a prise en compte d'un terme de convection pour les
821 // conditions aux limites de Neumann_sortie_libre seulement
822
823 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
824 {
825 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
826
827 if (sub_type(Neumann_sortie_libre,la_cl.valeur()))
828 {
829 const Neumann_sortie_libre& la_sortie_libre = ref_cast(Neumann_sortie_libre,la_cl.valeur());
830 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
831 int num1 = le_bord.num_premiere_face();
832 int num2 = num1 + le_bord.nb_faces();
833 for (num_face=num1; num_face<num2; num_face++)
834 {
835 psc =0;
836 for (i=0; i<dimension; i++)
837 psc += la_vitesse.valeurs()(num_face,i)*face_normales(num_face,i)*porosite_face(num_face);
838 if (psc>0)
839 for (i=0; i<ncomp_ch_transporte; i++)
840 {
841 resu(num_face,i) -= psc*transporte(num_face,i);
842 flux_b(num_face,i) -= psc*transporte(num_face,i);
843 }
844 else
845 {
846 for (i=0; i<ncomp_ch_transporte; i++)
847 {
848 resu(num_face,i) -= psc*la_sortie_libre.val_ext(num_face-num1,i);
849 flux_b(num_face,i) -= psc*la_sortie_libre.val_ext(num_face-num1,i);
850 }
851 fluent_(num_face) -= psc;
852 }
853 }
854 }
855 else if (sub_type(Periodique,la_cl.valeur()))
856 {
857 const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
858 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
859 int num1 = le_bord.num_premiere_face();
860 int num2 = num1 + le_bord.nb_faces();
861 IntVect fait(le_bord.nb_faces());
862 fait = 0;
863 for (num_face=num1; num_face<num2; num_face++)
864 {
865 if (fait[num_face-num1] == 0)
866 {
867 voisine = la_cl_perio.face_associee(num_face-num1) + num1;
868 for (int comp=0; comp<ncomp_ch_transporte; comp++)
869 {
870 diff1 = resu(num_face,comp)-tab(nb_faces_perio,comp);
871 diff2 = resu(voisine,comp)-tab(nb_faces_perio+voisine-num_face,comp);
872 resu(voisine,comp) += diff1;
873 resu(num_face,comp) += diff2;
874 flux_b(voisine,comp) += diff1;
875 flux_b(num_face,comp) += diff2;
876 }
877 fait[num_face-num1]= 1;
878 fait[voisine-num1] = 1;
879 }
880 nb_faces_perio++;
881 }
882 }
883 }
884 modifier_flux(*this);
885 return resu;
886}
Champ_Inc_base
Classe Champ_Inc_base.
Definition Champ_Inc_base.h:53
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
Definition Champ_Inc_base.h:77
Champ_base::valeur_a_sommet_compo
virtual double valeur_a_sommet_compo(int, int, int) const
renvoi la compo eme corrdonne des valeurs a l'element le_poly au sommet sommet
Definition Champ_base.cpp:321
Cond_lim
classe Cond_lim Classe generique servant a representer n'importe quelle classe
Definition Cond_lim.h:31
Domaine_32_64::creer_tableau_elements
virtual void creer_tableau_elements(Array_base &, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT) const
creation d'un tableau parallele de valeurs aux elements.
Definition Domaine.cpp:851
Domaine_Cl_VEF
Definition Domaine_Cl_VEF.h:40
Domaine_Cl_VEF::type_elem_Cl
int type_elem_Cl(int i) const
Definition Domaine_Cl_VEF.h:62
Domaine_Cl_VEF::normales_facettes_Cl
DoubleTab & normales_facettes_Cl()
Definition Domaine_Cl_VEF.h:56
Domaine_Cl_dis_base
classe Domaine_Cl_dis_base Les objets Domaine_Cl_dis_base representent les conditions aux limites
Definition Domaine_Cl_dis_base.h:37
Domaine_Cl_dis_base::les_conditions_limites
const Cond_lim & les_conditions_limites(int) const
Renvoie la i-ieme condition aux limites.
Definition Domaine_Cl_dis_base.cpp:386
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:54
Domaine_VEF::rang_elem_non_std
IntVect & rang_elem_non_std()
Definition Domaine_VEF.h:86
Domaine_VEF::type_elem
const Elem_VEF_base & type_elem() const
Definition Domaine_VEF.h:75
Domaine_VEF::facette_normales
auto & facette_normales()
Definition Domaine_VEF.h:84
Domaine_VF::nb_faces
int nb_faces() const
renvoie le nombre global de faces.
Definition Domaine_VF.h:471
Domaine_VF::face_normales
virtual double face_normales(int face, int comp) const
Definition Domaine_VF.h:47
Domaine_VF::xv
double xv(int num_face, int k) const
Definition Domaine_VF.h:76
Domaine_VF::volumes
double volumes(int i) const
Definition Domaine_VF.h:113
Domaine_VF::elem_faces
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
Domaine_VF::xp
double xp(int num_elem, int k) const
Definition Domaine_VF.h:77
Domaine_VF::premiere_face_int
int premiere_face_int() const
une face est interne ssi elle separe deux elements.
Definition Domaine_VF.h:463
Domaine_VF::face_voisins
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
Definition Domaine_VF.h:418
Domaine_VF::nb_faces_bord
int nb_faces_bord() const
renvoie le nombre de faces sur lesquelles sont appliquees les conditions limites :
Definition Domaine_VF.h:513
Domaine_dis_base
classe Domaine_dis_base Cette classe est la base de la hierarchie des domaines discretisees.
Definition Domaine_dis_base.h:39
Domaine_dis_base::nb_elem_tot
int nb_elem_tot() const
Definition Domaine_dis_base.h:50
Domaine_dis_base::nb_elem
int nb_elem() const
Definition Domaine_dis_base.h:49
Domaine_dis_base::nb_front_Cl
int nb_front_Cl() const
Definition Domaine_dis_base.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::calcul_vc
virtual void calcul_vc(const ArrOfInt &, ArrOfDouble &, const ArrOfDouble &, const DoubleTab &, const Champ_Inc_base &, int, const DoubleVect &) const =0
Elem_VEF_base::KEL
const IntTab & KEL() const
Definition Elem_VEF_base.h:31
Elem_VEF_base::nb_facette
virtual int nb_facette() const =0
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Equation_base::milieu
virtual const Milieu_base & milieu() const =0
Front_VF
class Front_VF
Definition Front_VF.h:36
Front_VF::nb_faces
int nb_faces() const
Definition Front_VF.h:53
Front_VF::num_premiere_face
int num_premiere_face() const
Definition Front_VF.h:63
Milieu_base::porosite_face
DoubleVect & porosite_face()
Definition Milieu_base.h:62
MorEqn::equation
const Equation_base & equation() const
Renvoie la reference sur l'equation pointe par MorEqn::mon_equation.
Definition MorEqn.h:62
Neumann_sortie_libre
classe Neumann_sortie_libre Cette classe represente une frontiere ouverte sans vitesse imposee
Definition Neumann_sortie_libre.h:34
Neumann_sortie_libre::val_ext
double val_ext(int i) const override
Renvoie la valeur de la i-eme composante du champ impose a l'exterieur de la frontiere.
Definition Neumann_sortie_libre.cpp:120
Nom
class Nom Une chaine de caractere pour nommer les objets de TRUST
Definition Nom.h:31
Objet_U::dimension
static int dimension
Definition Objet_U.h:99
Objet_U::que_suis_je
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
Definition Objet_U.cpp:104
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::printOn
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
Op_Conv_VEF_base
class Op_Conv_VEF_base
Definition Op_Conv_VEF_base.h:37
Op_Conv_VEF_base::fluent_
DoubleVect fluent_
Definition Op_Conv_VEF_base.h:61
Op_Conv_VEF_base::associer
void associer(const Domaine_dis_base &, const Domaine_Cl_dis_base &, const Champ_Inc_base &) override
Definition Op_Conv_VEF_base.cpp:231
Op_Conv_VEF_base::vitesse
const Champ_Inc_base & vitesse() const
Definition Op_Conv_VEF_base.h:41
Op_Conv_kschemas_VEF
class Op_Conv_kschemas_VEF
Definition Op_Conv_kschemas_VEF.h:34
Op_Conv_kschemas_VEF::associer
void associer(const Domaine_dis_base &, const Domaine_Cl_dis_base &, const Champ_Inc_base &) override
Definition Op_Conv_kschemas_VEF.cpp:33
Op_Conv_kschemas_VEF::K
double K
Definition Op_Conv_kschemas_VEF.h:44
Op_Conv_kschemas_VEF::ajouter
DoubleTab & ajouter(const DoubleTab &, DoubleTab &) const override
Definition Op_Conv_kschemas_VEF.cpp:157
Op_VEF_Face::modifier_flux
void modifier_flux(const Operateur_base &) const
Definition Op_VEF_Face.cpp:338
Operateur_base::flux_bords_
DoubleTab flux_bords_
Definition Operateur_base.h:140
Periodique
classe Periodique Cette classe represente une condition aux limites periodique.
Definition Periodique.h:31
Periodique::face_associee
int face_associee(int i) const
Definition Periodique.h:35
Process::is_parallel
static bool is_parallel()
Definition Process.cpp:110
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
TRUSTTab::resize
void resize(_SIZE_ n, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTTab.tpp:469
TRUSTVect::line_size
int line_size() const
Definition TRUSTVect.tpp:67
TRUSTVect::echange_espace_virtuel
virtual void echange_espace_virtuel(IsExchangeBlocking exchange_type=IsExchangeBlocking::DefaultBlocking, const std::string kernel_name="noname")
Definition TRUSTVect.tpp:282