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