TrioCFD 1.9.9_beta
TrioCFD documentation
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
23Sortie& Op_Conv_kschemas_centre_VEF::printOn(Sortie& s ) const
24{
25 return s << que_suis_je() ;
26}
27
28Entree& Op_Conv_kschemas_centre_VEF::readOn(Entree& s )
29{
30 return s ;
31}
32
33
34//////////////////////////////////////////////////////////////
35// Functions for kschemas_centre.
36////////////////////////////////////////////////////////////////
37
38// convkschemas_centre : utility function for 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 //**** Case where the upwind scheme is applied
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 // extrapolate the upwind face flux onto the line connecting the midpoints of the two 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 // **** Case where the centered scheme is applied
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 // Determine deltat phi0 and phi1
118 deltat0 += gradient_elem(poly1,0,i)*rx(i);
119 deltat1 -= gradient_elem(poly2,0,i)*rx(i);
120 }
121
122
123 // Flux computation
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 // Flux computation
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 // For the following 3 element types, the number of vertices equals the number of faces
198 // composing the geometric element.
199 // Problem with hexa: 8 vertices and 6 faces, so using the face[i] array no longer works
200 // the method retained to avoid computing the velocity at vertices without
201 // shape functions is therefore not applicable,
202 // for hexa elements the face-to-vertex mapping is not directly accessible.
203 // The Face=>vertices table exists but not the inverse.
204 // Too costly; for now porosity is not extended to hexa elements.
205
206 int istetra=0;
207 const Elem_VEF_base& type_elemvef= domaine_VEF.type_elem();
208 Nom nom_elem=type_elemvef.que_suis_je();
209 if ((nom_elem=="Tetra_VEF")||(nom_elem=="Tri_VEF"))
210 istetra=1;
211
212 // For convection processing, standard polyhedra (those not "seeing" boundary conditions)
213 // are distinguished from non-standard polyhedra that have at least one boundary face.
214 // A standard polyhedron has n facets on which the convection scheme is applied.
215 // For a non-standard polyhedron with Dirichlet boundary conditions,
216 // some of the facets are carried by the boundary faces.
217 // In short, for a polyhedron, the convection treatment depends
218 // on the type (triangle, tetrahedron, ...) and the number of Dirichlet faces.
219
220 double psc;
221 int poly,poly1,poly2,face_adj,fa7,i,j,n_bord;
222 int num_face, rang ,itypcl;
223 int num10,num20,num3,num_som;
224 int ncomp_ch_transporte;
225
226 if (transporte.nb_dim() == 1)
227 ncomp_ch_transporte=1;
228 else
229 ncomp_ch_transporte= transporte.dimension(1);
230
231 int fac,elem1,elem2,comp0;
232 int nb_faces_ = domaine_VEF.nb_faces();
233 IntVect face(nfac);
234
235 DoubleVect flux(ncomp_ch_transporte);
236 DoubleVect fluxsom(ncomp_ch_transporte);
237 DoubleVect fluxg(ncomp_ch_transporte);
238
239
240 // Special treatment for periodic faces
241 int nb_faces_perio = 0;
242 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
243 {
244 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
245 if (sub_type(Periodique,la_cl.valeur()))
246 {
247 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
248 int num1 = le_bord.num_premiere_face();
249 int num2 = num1 + le_bord.nb_faces();
250 for (num_face=num1; num_face<num2; num_face++)
251 nb_faces_perio++;
252 }
253 }
254
255 DoubleTab tab;
256 if (ncomp_ch_transporte == 1)
257 tab.resize(nb_faces_perio);
258 else
259 tab.resize(nb_faces_perio,ncomp_ch_transporte);
260
261 nb_faces_perio=0;
262 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
263 {
264 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
265 if (sub_type(Periodique,la_cl.valeur()))
266 {
267 // const Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
268 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
269 int num1 = le_bord.num_premiere_face();
270 int num2 = num1 + le_bord.nb_faces();
271 for (num_face=num1; num_face<num2; num_face++)
272 {
273 if (ncomp_ch_transporte == 1)
274 tab(nb_faces_perio) = resu(num_face);
275 else
276 for (int comp=0; comp<ncomp_ch_transporte; comp++)
277 tab(nb_faces_perio,comp) = resu(num_face,comp);
278 nb_faces_perio++;
279 }
280 }
281 }
282
283
284 ///////////////////////////////////////////////////////////////////////////////////////////////
285 // <
286 // gradient computation; < [ Ujp*np/vol(j) ]
287 // j
288 ////////////////////////////////////////////////////////////////////////////////////////////////
289 // Create the virtual space for gradient_elem
290 DoubleTab gradient_elem(0, ncomp_ch_transporte, dimension);
291 // (du/dx du/dy dv/dx dv/dy) for a polyhedron
292 domaine_VEF.domaine().creer_tableau_elements(gradient_elem);
293 // Loop over faces
294
295
296 for (fac=0; fac< premiere_face_int; fac++)
297 {
298 elem1=face_voisins(fac,0);
299 if(ncomp_ch_transporte==1)
300 for (i=0; i<dimension; i++)
301 {
302 gradient_elem(elem1, 0, i) +=
303 face_normales(fac,i)*transporte(fac);
304 }
305 else
306 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
307 for (i=0; i<dimension; i++)
308 gradient_elem(elem1, comp0, i) +=
309 face_normales(fac,i)*transporte(fac,comp0);
310 // dUcomp/dXi
311 } // end of for faces
312
313 for (; fac<nb_faces_; fac++)
314 {
315 elem1=face_voisins(fac,0);
316 elem2=face_voisins(fac,1);
317 if(ncomp_ch_transporte==1)
318 for (i=0; i<dimension; i++)
319 {
320 gradient_elem(elem1, 0, i) +=
321 face_normales(fac,i)*transporte(fac);
322 gradient_elem(elem2, 0, i) -=
323 face_normales(fac,i)*transporte(fac);
324 }
325 else
326 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
327 for (i=0; i<dimension; i++)
328 {
329 gradient_elem(elem1, comp0, i) +=
330 face_normales(fac,i)*transporte(fac,comp0);
331 gradient_elem(elem2, comp0, i) -=
332 face_normales(fac,i)*transporte(fac,comp0);
333 }
334 // dUcomp/dXi
335 } // end of for faces
336
337 for (int elem=0; elem<nb_elem; elem++)
338 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
339 for (i=0; i<dimension; i++)
340 gradient_elem(elem,comp0,i) /= volumes(elem);
341
342 gradient_elem.echange_espace_virtuel();
343
344 ////////////////////////////////////////////////////////////////////////////////////
345 // gradient_elem is now available per element
346 ////////////////////////////////////////////////////////////////////////////////////
347 DoubleVect coord_som(dimension);
348 DoubleVect vs(dimension);
349 DoubleVect vc(dimension);
350 DoubleTab vsom(nsom,dimension);
351 DoubleVect cc(dimension);
352 double xm;
353
354 // Reset the array used for stability time step computation
355 fluent_ = 0;
356
357 // Non-standard polyhedra are grouped in 2 categories in Domaine_VEF:
358 // - boundary and joint polyhedra
359 // - boundary and non-joint polyhedra
360 // Process polyhedra following the order in which they appear
361 // in the domain
362
363 //////////////////////////////////////////////////////////////////////////////////////
364 // loop over polyhedra
365 //////////////////////////////////////////////////////////////////////////////////////
366
367 const IntTab& KEL=type_elemvef.KEL();
368 for (poly=0; poly<nb_elem_tot; poly++)
369 {
370 //Cerr << "poly = " << poly << finl;
371
372 rang = rang_elem_non_std(poly);
373 if (rang==-1)
374 itypcl=0;
375 else
376 itypcl=domaine_Cl_VEF.type_elem_Cl(rang);
377
378 // compute face indices of the polyhedron
379
380 for (face_adj=0; face_adj<nfac; face_adj++)
381 {
382 face(face_adj)= elem_faces(poly,face_adj);
383 //Cerr << "les faces de l'elements sont : " << face(face_adj) << finl;
384 }
385
386 int scom;
387 DoubleVect rx0(dimension);
388
389 // compute velocity at the vertices of the polyhedra
390
391 for (j=0; j<dimension; j++)
392 {
393 vs(j) = la_vitesse.valeurs()(face(0),j)*porosite_face(face(0));
394 for (i=1; i<nfac; i++)
395 vs(j)+= la_vitesse.valeurs()(face(i),j)*porosite_face(face(i));
396 }
397
398 //int ncomp;
399 if (istetra==1)
400 {
401 for (j=0; j<nsom; j++)
402 {
403 for (int ncomp=0; ncomp<Objet_U::dimension; ncomp++)
404 vsom(j,ncomp) =vs[ncomp] - Objet_U::dimension*la_vitesse.valeurs()(face[j],ncomp)*porosite_face(face[j]);
405 }
406 }
407 else
408 {
409 // to make this valid for hexahedral elements
410 // Use the shape functions implemented in the class Champs_P1_impl or Champs_Q1_impl
411 int ncomp;
412 for (j=0; j<nsom; j++)
413 {
414 num_som = domaine.sommet_elem(poly,j);
415 for (ncomp=0; ncomp<dimension; ncomp++)
416 {
417 vsom(j,ncomp) = la_vitesse.valeur_a_sommet_compo(num_som,poly,ncomp);
418 }
419 }
420 }
421
422 // compute velocity at the center of gravity
423
424 domaine_VEF.type_elem().calcul_vc(face,vc,vs,vsom,vitesse(),itypcl,porosite_face);
425
426
427 // Loop over the facets of the non-standard polyhedron:
428
429 for (fa7=0; fa7<nfa7; fa7++)
430 {
431 //Cerr << "the facet being processed is " << fa7 << finl;
432 // fa7 separates num1 and num2. num3 is the third face (2D).
433
434 num10 = face(KEL(0,fa7));
435 num20 = face(KEL(1,fa7));
436 num3 = face(KEL(2,fa7));
437
438 // Determine the neighbouring elements at faces num1 and num2
439
440 poly1 = face_voisins(num10,0);
441 if (poly1==poly)
442 {
443 poly1 = face_voisins(num10,1);
444 }
445
446 poly2 = face_voisins(num20,0);
447 if (poly2==poly)
448 {
449 poly2 = face_voisins(num20,1);
450 }
451
452 scom = les_Polys(poly,KEL(2,fa7));
453
454 // compute rx0, distance between centers of 'num i' faces
455
456 for (i=0; i<dimension; i++)
457 rx0(i) = xv(num20,i)-xv(num10,i);
458
459 // facet normals
460
461 if (rang==-1)
462 for (i=0; i<dimension; i++)
463 cc[i] = facette_normales(poly, fa7, i);
464 else
465 for (i=0; i<dimension; i++)
466 cc[i] = normales_facettes_Cl(rang,fa7,i);
467
468
469 /////////////////////////////////////////////////////////////////////////
470 // Process the point where velocity = 0.5(vitsommet + vitmilieu)
471 /////////////////////////////////////////////////////////////////////////
472
473
474 for (i=0; i<nb_som_facette-1; i++)
475 {
476
477 //////////////////////////////////////////////////////////////////////////
478 //Determine PhiIJ at the midpoint between the vertex and the midpoint of num3
479 /////////////////////////////////////////////////////////////////////////
480
481 psc = 0;
482 for (j=0; j<dimension; j++)
483 psc+=((vsom(KEL(i+2,fa7),j) + la_vitesse.valeurs()(num3,j) * porosite_face(num3)))*cc[j];
484 psc *=0.5;
485 for (j=0; j<dimension; j++)
486 coord_som(j)=coord(scom,j);
487
488 convkschemas_centre(K,ncomp_ch_transporte,dimension,poly,poly1,poly2,num10,num20,psc,transporte,
489 fluent_,flux,rx0,gradient_elem,coord_som);
490
491
492 ////////////////////////////////////////////////////////////////////////////////////////////////////////
493 //Gradient limiter. Computed simultaneously with the flux computation.
494 // gradient(K0) = teta*gradient(K0)+ (1-teta)gradient(K1 or K2)
495 ////////////////////////////////////////////////////////////////////////////////////////////////////////
496
497 double teta;
498 teta = 1.;
499
500 /////////////////////////////////////////////////////////////////////////
501 // Process vertices that are also vertices of the polyhedron
502 /////////////////////////////////////////////////////////////////////////
503
504 if (ncomp_ch_transporte == 1)
505 {
506 for (j=0; j<dimension; j++)
507 {
508 xm = 0.5 *(coord(scom,j)+xv(num3,j));
509 if (psc >= 0)
510 {
511 if (poly1==-1)
512 {
513 fluxsom(0) += gradient_elem(poly,0,j)*(coord(scom,j)-xm);
514 }
515 else
516 {
517 fluxsom(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly1,0,j))*(coord(scom,j)-xm);
518 }
519 }
520
521 else
522 {
523 if (poly2==-1)
524 {
525 fluxsom(0) += gradient_elem(poly,0,j)*(coord(scom,j)-xm);
526 }
527 else
528 {
529 fluxsom(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly2,0,j))*(coord(scom,j)-xm);
530 }
531 }
532 }
533 fluxsom(0) += flux(0);
534 fluxsom(0) *= psc;
535 }
536
537 else
538 {
539 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
540 {
541 fluxsom(comp0) = flux(comp0);
542 for (j=0; j<dimension; j++)
543 {
544 xm = 0.5 *(coord(scom,j)+xv(num3,j));
545 if (psc >= 0)
546 {
547 if (poly1==-1)
548 {
549 fluxsom(comp0) += gradient_elem(poly,comp0,j)*(coord(scom,j)-xm);
550 }
551 else
552 {
553 fluxsom(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly1,comp0,j))*(coord(scom,j)-xm);
554 }
555 }
556
557 else
558 {
559 if (poly2==-1)
560 {
561 fluxsom(comp0) += gradient_elem(poly,comp0,j)*(coord(scom,j)-xm);
562 }
563 else
564 {
565 fluxsom(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly2,comp0,j))*(coord(scom,j)-xm);
566 }
567
568 }
569
570 }
571 fluxsom(comp0) *= psc;
572 // Cerr << "fluxsom(" << comp << ") = " << fluxsom(comp) << finl;
573 }
574 }
575
576
577 ////////////////////////////////////////////////////////////////////////////
578 // Process the center of gravity
579 ////////////////////////////////////////////////////////////////////////////
580
581 if (ncomp_ch_transporte == 1)
582 {
583 for (j=0; j<dimension; j++)
584 {
585 xm = 0.5 *(coord(scom,j)+xv(num3,j));
586 if (psc >= 0)
587 {
588 if (poly1==-1)
589 {
590 fluxg(0) += gradient_elem(poly,0,j)*(xg(poly,j)-xm);
591 }
592 else
593 {
594 fluxg(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly1,0,j))*(xg(poly,j)-xm);
595 }
596 }
597
598 else
599 {
600 if (poly2==-1)
601 {
602 fluxg(0) += gradient_elem(poly,0,j)*(xg(poly,j)-xm);
603 }
604 else
605 {
606 fluxg(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly2,0,j))*(xg(poly,j)-xm);
607 }
608 }
609
610 }
611
612 fluxg(0) += flux(0);
613 fluxg(0) *= psc;
614 }
615 else
616 {
617 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
618 {
619 fluxg(comp0) = flux(comp0);
620 for (j=0; j<dimension; j++)
621 {
622 xm = 0.5 *(coord(scom,j)+xv(num3,j));
623 if (psc >= 0)
624 {
625 if (poly1==-1)
626 {
627 fluxg(comp0) += gradient_elem(poly,comp0,j)*(xg(poly,j)-xm);
628 }
629 else
630 {
631 fluxg(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly1,comp0,j))*(xg(poly,j)-xm);
632 }
633 }
634
635 else
636 {
637 if (poly2==-1)
638 {
639 fluxg(comp0) += gradient_elem(poly,comp0,j)*(xg(poly,j)-xm);
640 }
641 else
642 {
643 fluxg(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly2,comp0,j))*(xg(poly,j)-xm);
644 }
645 }
646
647 }
648
649 fluxg(comp0) *= psc;
650 // Cerr << "fluxg(" << comp << ") = " << fluxg(comp) << finl;
651 }
652 }
653
654 //////////////////////////////////////////////////////////////////////////////
655 // Integration of u.n.flux
656 /////////////////////////////////////////////////////////////////////////////
657
658 if (ncomp_ch_transporte == 1)
659 {
660 resu(num10) -=flux(0)*psc;
661 resu(num20) += flux(0)*psc;
662
663 }
664 else
665 {
666 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
667 {
668 resu(num10,comp0) -= flux(comp0)*psc;
669 resu(num20,comp0) += flux(comp0)*psc;
670 }
671 }
672 }
673 }
674 } // end of loop
675
676 // Processing of joint elements of thickness 1
677 for (poly=0; poly<nb_elem_tot; poly++)
678 {
679 // Check if a face of the polyhedron is a joint face
680 for (face_adj=0; face_adj<nfac; face_adj++)
681 if(face_adj<nb_faces) break;
682 if(face_adj<nfac)
683 {
684 rang = rang_elem_non_std(poly);
685 if (rang==-1)
686 itypcl=0;
687 else
688 itypcl=domaine_Cl_VEF.type_elem_Cl(rang);
689
690 // compute face indices of the polyhedron
691
692 for (face_adj=0; face_adj<nfac; face_adj++)
693 {
694 face(face_adj)= elem_faces(poly,face_adj);
695 //Cerr << "les faces de l'elements sont : " << face(face_adj) << finl;
696 }
697
698 int scom;
699 DoubleVect rx0(dimension);
700
701 // compute velocity at the vertices of the polyhedra
702 for (j=0; j<dimension; j++)
703 {
704 vs(j) = la_vitesse.valeurs()(face(0),j)*porosite_face(face(0));
705 for (i=1; i<nfac; i++)
706 vs(j)+= la_vitesse.valeurs()(face(i),j)*porosite_face(face(i));
707 }
708
709
710 if (istetra==1)
711 {
712 for (j=0; j<nsom; j++)
713 {
714 for (int ncomp=0; ncomp<Objet_U::dimension; ncomp++)
715 vsom(j,ncomp) =vs[ncomp] - Objet_U::dimension*la_vitesse.valeurs()(face[j],ncomp)*porosite_face(face[j]);
716 }
717 }
718 else
719 {
720 // to be valid for hexa elements
721 // use shape functions implemented in Champs_P1_impl or Champs_Q1_impl
722 int ncomp;
723 for (j=0; j<nsom; j++)
724 {
725 num_som = domaine.sommet_elem(poly,j);
726 for (ncomp=0; ncomp<dimension; ncomp++)
727 {
728 vsom(j,ncomp) = la_vitesse.valeur_a_sommet_compo(num_som,poly,ncomp);
729 }
730 }
731 }
732
733 // compute velocity at the center of gravity
734
735 domaine_VEF.type_elem().calcul_vc(face,vc,vs,vsom,vitesse(),itypcl,porosite_face);
736
737
738 // Loop over the facets of the non-standard polyhedron:
739
740 for (fa7=0; fa7<nfa7; fa7++)
741 {
742 //Cerr << "the facet being processed is " << fa7 << finl;
743 // fa7 separates num1 and num2. num3 is the third face (2D).
744
745 num10 = face(KEL(0,fa7));
746 num20 = face(KEL(1,fa7));
747 num3 = face(KEL(2,fa7));
748
749 // Determine the neighbouring elements at faces num1 and num2
750
751 poly1 = face_voisins(num10,0);
752 if (poly1==poly)
753 {
754 poly1 = face_voisins(num10,1);
755 }
756
757 poly2 = face_voisins(num20,0);
758 if (poly2==poly)
759 {
760 poly2 = face_voisins(num20,1);
761 }
762
763 scom = les_Polys(poly,KEL(2,fa7));
764
765 // compute rx0, distance between centers of 'num i' faces
766
767 for (i=0; i<dimension; i++)
768 rx0(i) = xv(num20,i)-xv(num10,i);
769
770 // facet normals
771
772 if (rang==-1)
773 for (i=0; i<dimension; i++)
774 cc[i] = facette_normales(poly, fa7, i);
775 else
776 for (i=0; i<dimension; i++)
777 cc[i] = normales_facettes_Cl(rang,fa7,i);
778
779
780 /////////////////////////////////////////////////////////////////////////
781 // Process the point where velocity = 0.5(vitsommet + vitmilieu)
782 /////////////////////////////////////////////////////////////////////////
783
784
785 for (i=0; i<nb_som_facette-1; i++)
786 {
787
788 //////////////////////////////////////////////////////////////////////////
789 //Determine PhiIJ at the midpoint between the vertex and the midpoint of num3
790 /////////////////////////////////////////////////////////////////////////
791
792 psc = 0;
793 for (j=0; j<dimension; j++)
794 psc+=((vsom(KEL(i+2,fa7),j) + la_vitesse.valeurs()(num3,j) * porosite_face(num3)))*cc[j];
795 psc *=0.5;
796 for (j=0; j<dimension; j++)
797 coord_som(j)=coord(scom,j);
798 convkschemas_centre(K,ncomp_ch_transporte,dimension,poly,poly1,poly2,num10,num20,psc,transporte,
799 fluent_,flux,rx0,gradient_elem,coord_som);
800
801
802 ////////////////////////////////////////////////////////////////////////////////////////////////////////
803 //Gradient limiter. Computed simultaneously with the flux computation.
804 // gradient(K0) = teta*gradient(K0)+ (1-teta)gradient(K1 or K2)
805 ////////////////////////////////////////////////////////////////////////////////////////////////////////
806
807 double teta;
808 teta = 1.;
809
810 /////////////////////////////////////////////////////////////////////////
811 // Process vertices that are also vertices of the polyhedron
812 /////////////////////////////////////////////////////////////////////////
813
814 if (ncomp_ch_transporte == 1)
815 {
816 for (j=0; j<dimension; j++)
817 {
818 xm = 0.5 *(coord(scom,j)+xv(num3,j));
819 if (psc >= 0)
820 {
821 if (poly1==-1)
822 {
823 fluxsom(0) += gradient_elem(poly,0,j)*(coord(scom,j)-xm);
824 }
825 else
826 {
827 fluxsom(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly1,0,j))*(coord(scom,j)-xm);
828 }
829 }
830
831 else
832 {
833 if (poly2==-1)
834 {
835 fluxsom(0) += gradient_elem(poly,0,j)*(coord(scom,j)-xm);
836 }
837 else
838 {
839 fluxsom(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly2,0,j))*(coord(scom,j)-xm);
840 }
841 }
842 }
843 fluxsom(0) += flux(0);
844 fluxsom(0) *= psc;
845 }
846
847 else
848 {
849 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
850 {
851 fluxsom(comp0) = flux(comp0);
852 for (j=0; j<dimension; j++)
853 {
854 xm = 0.5 *(coord(scom,j)+xv(num3,j));
855 if (psc >= 0)
856 {
857 if (poly1==-1)
858 {
859 fluxsom(comp0) += gradient_elem(poly,comp0,j)*(coord(scom,j)-xm);
860 }
861 else
862 {
863 fluxsom(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly1,comp0,j))*(coord(scom,j)-xm);
864 }
865 }
866
867 else
868 {
869 if (poly2==-1)
870 {
871 fluxsom(comp0) += gradient_elem(poly,comp0,j)*(coord(scom,j)-xm);
872 }
873 else
874 {
875 fluxsom(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly2,comp0,j))*(coord(scom,j)-xm);
876 }
877
878 }
879
880 }
881 fluxsom(comp0) *= psc;
882 // Cerr << "fluxsom(" << comp << ") = " << fluxsom(comp) << finl;
883 }
884 }
885
886
887 ////////////////////////////////////////////////////////////////////////////
888 // Process the center of gravity
889 ////////////////////////////////////////////////////////////////////////////
890
891 if (ncomp_ch_transporte == 1)
892 {
893 for (j=0; j<dimension; j++)
894 {
895 xm = 0.5 *(coord(scom,j)+xv(num3,j));
896 if (psc >= 0)
897 {
898 if (poly1==-1)
899 {
900 fluxg(0) += gradient_elem(poly,0,j)*(xg(poly,j)-xm);
901 }
902 else
903 {
904 fluxg(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly1,0,j))*(xg(poly,j)-xm);
905 }
906 }
907
908 else
909 {
910 if (poly2==-1)
911 {
912 fluxg(0) += gradient_elem(poly,0,j)*(xg(poly,j)-xm);
913 }
914 else
915 {
916 fluxg(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly2,0,j))*(xg(poly,j)-xm);
917 }
918 }
919
920 }
921
922 fluxg(0) += flux(0);
923 fluxg(0) *= psc;
924 }
925 else
926 {
927 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
928 {
929 fluxg(comp0) = flux(comp0);
930 for (j=0; j<dimension; j++)
931 {
932 xm = 0.5 *(coord(scom,j)+xv(num3,j));
933 if (psc >= 0)
934 {
935 if (poly1==-1)
936 {
937 fluxg(comp0) += gradient_elem(poly,comp0,j)*(xg(poly,j)-xm);
938 }
939 else
940 {
941 fluxg(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly1,comp0,j))*(xg(poly,j)-xm);
942 }
943 }
944
945 else
946 {
947 if (poly2==-1)
948 {
949 fluxg(comp0) += gradient_elem(poly,comp0,j)*(xg(poly,j)-xm);
950 }
951 else
952 {
953 fluxg(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly2,comp0,j))*(xg(poly,j)-xm);
954 }
955 }
956
957 }
958
959 fluxg(comp0) *= psc;
960 // Cerr << "fluxg(" << comp << ") = " << fluxg(comp) << finl;
961 }
962 }
963
964 //////////////////////////////////////////////////////////////////////////////
965 // Integration of u.n.flux
966 /////////////////////////////////////////////////////////////////////////////
967
968 if (ncomp_ch_transporte == 1)
969 {
970 resu(num10) -= flux(0)*psc;
971 resu(num20) += flux(0)*psc;
972
973 }
974 else
975 {
976 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
977 {
978 resu(num10,comp0) -= flux(comp0)*psc;
979 resu(num20,comp0) +=flux(comp0)*psc;
980 }
981 }
982 }
983 }
984 }
985 } // end of loop
986 int voisine;
987 nb_faces_perio = 0;
988 double diff1,diff2;
989
990 // Size the array of convective fluxes at the domain boundary
991 DoubleTab& flux_b = flux_bords_;
992 flux_b.resize(domaine_VEF.nb_faces_bord(),ncomp_ch_transporte);
993 flux_b = 0.;
994
995 // Loop over boundary conditions
996 // a convection term is taken into account for
997 // Neumann_sortie_libre boundary conditions only
998
999 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
1000 {
1001 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1002
1003 if (sub_type(Neumann_sortie_libre,la_cl.valeur()))
1004 {
1005 const Neumann_sortie_libre& la_sortie_libre = ref_cast(Neumann_sortie_libre,la_cl.valeur());
1006 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1007 int num1 = le_bord.num_premiere_face();
1008 int num2 = num1 + le_bord.nb_faces();
1009 for (num_face=num1; num_face<num2; num_face++)
1010 {
1011 psc =0;
1012 for (i=0; i<dimension; i++)
1013 psc += la_vitesse.valeurs()(num_face,i)*face_normales(num_face,i)*porosite_face(num_face);
1014 if (psc>0)
1015 if (ncomp_ch_transporte == 1)
1016 {
1017 resu(num_face) -= psc*transporte(num_face);
1018 flux_b(num_face,0) -= psc*transporte(num_face);
1019 }
1020 else
1021 for (i=0; i<ncomp_ch_transporte; i++)
1022 {
1023 resu(num_face,i) -= psc*transporte(num_face,i);
1024 flux_b(num_face,i) -= psc*transporte(num_face,i);
1025 }
1026 else
1027 {
1028 if (ncomp_ch_transporte == 1)
1029 {
1030 resu(num_face) -= psc*la_sortie_libre.val_ext(num_face-num1);
1031 flux_b(num_face,0) -= psc*la_sortie_libre.val_ext(num_face-num1);
1032 }
1033 else
1034 for (i=0; i<ncomp_ch_transporte; i++)
1035 {
1036 resu(num_face,i) -= psc*la_sortie_libre.val_ext(num_face-num1,i);
1037 flux_b(num_face,i) -= psc*la_sortie_libre.val_ext(num_face-num1,i);
1038 }
1039 fluent_(num_face) -= psc;
1040 }
1041 }
1042 }
1043 else if (sub_type(Periodique,la_cl.valeur()))
1044 {
1045 const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
1046 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1047 int num1 = le_bord.num_premiere_face();
1048 int num2 = num1 + le_bord.nb_faces();
1049 IntVect fait(le_bord.nb_faces());
1050 fait = 0;
1051 for (num_face=num1; num_face<num2; num_face++)
1052 {
1053 if (fait[num_face-num1] == 0)
1054 {
1055 voisine = la_cl_perio.face_associee(num_face-num1) + num1;
1056
1057 if (ncomp_ch_transporte == 1)
1058 {
1059 diff1 = resu(num_face)-tab(nb_faces_perio);
1060 diff2 = resu(voisine)-tab(nb_faces_perio+voisine-num_face);
1061 resu(voisine) += diff1;
1062 resu(num_face) += diff2;
1063 flux_b(voisine,0) += diff1;
1064 flux_b(num_face,0) += diff2;
1065 }
1066 else
1067 for (int comp=0; comp<ncomp_ch_transporte; comp++)
1068 {
1069 diff1 = resu(num_face,comp)-tab(nb_faces_perio,comp);
1070 diff2 = resu(voisine,comp)-tab(nb_faces_perio+voisine-num_face,comp);
1071 resu(voisine,comp) += diff1;
1072 resu(num_face,comp) += diff2;
1073 flux_b(num_face,comp) += diff1;
1074 flux_b(num_face,comp) += diff2;
1075 }
1076
1077 fait[num_face-num1]= 1;
1078 fait[voisine-num1] = 1;
1079 }
1080 nb_faces_perio++;
1081 }
1082 }
1083 }
1084 modifier_flux(*this);
1085 return resu;
1086
1087}
Champ_Inc_base
Class Champ_Inc_base.
Definition Champ_Inc_base.h:53
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Returns the array of field values at the current time.
Definition Champ_Inc_base.h:77
Champ_base::valeur_a_sommet_compo
virtual double valeur_a_sommet_compo(int, int, int) const
Returns the compo-th coordinate of the values at the element le_poly at the vertex sommet.
Definition Champ_base.cpp:321
Cond_lim
class Cond_lim Generic class used to represent any class
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
Creates a parallel array of values at elements.
Definition Domaine.cpp:850
Domaine_Cl_VEF
Definition Domaine_Cl_VEF.h:38
Domaine_Cl_VEF::type_elem_Cl
int type_elem_Cl(int i) const
Definition Domaine_Cl_VEF.h:60
Domaine_Cl_VEF::normales_facettes_Cl
DoubleTab & normales_facettes_Cl()
Definition Domaine_Cl_VEF.h:54
Domaine_Cl_dis_base::les_conditions_limites
const Cond_lim & les_conditions_limites(int) const
Returns the i-th boundary condition.
Definition Domaine_Cl_dis_base.cpp:386
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:53
Domaine_VEF::rang_elem_non_std
IntVect & rang_elem_non_std()
Definition Domaine_VEF.h:85
Domaine_VEF::type_elem
const Elem_VEF_base & type_elem() const
Definition Domaine_VEF.h:74
Domaine_VEF::facette_normales
auto & facette_normales()
Definition Domaine_VEF.h:83
Domaine_VF::nb_faces
int nb_faces() const
Returns the total number of 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
Returns the index of the i-th face of element num_elem; the face numbering convention is.
Definition Domaine_VF.h:542
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
A face is internal if and only if it separates two elements.
Definition Domaine_VF.h:463
Domaine_VF::face_voisins
int face_voisins(int num_face, int i) const
Returns the neighbouring element of num_face in direction i.
Definition Domaine_VF.h:418
Domaine_VF::nb_faces_bord
int nb_faces_bord() const
Returns the number of faces on which boundary conditions are applied:
Definition Domaine_VF.h:512
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
Returns the reference to the equation pointed to by MorEqn::mon_equation.
Definition MorEqn.h:62
Neumann_sortie_libre
Neumann_sortie_libre This class represents an open boundary without imposed velocity.
Definition Neumann_sortie_libre.h:34
Neumann_sortie_libre::val_ext
double val_ext(int i) const override
Returns the value of the i-th component of the field imposed on the exterior of the boundary.
Definition Neumann_sortie_libre.cpp:120
Nom
class Nom: a character string for naming TRUST objects.
Definition Nom.h:31
Objet_U::dimension
static int dimension
Definition Objet_U.h:94
Objet_U::que_suis_je
const Nom & que_suis_je() const
Returns the string identifying the class.
Definition Objet_U.cpp:104
Objet_U::readOn
virtual Entree & readOn(Entree &)
Reads an Objet_U from an input stream. Virtual method to override.
Definition Objet_U.cpp:289
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Writes the object to an output stream. Virtual method to override.
Definition Objet_U.cpp:278
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::vitesse
const Champ_Inc_base & vitesse() const
Definition Op_Conv_VEF_base.h:41
Op_Conv_kschemas_centre_VEF
class Op_Conv_kschemas_centre_VEF
Definition Op_Conv_kschemas_centre_VEF.h:33
Op_Conv_kschemas_centre_VEF::K
double K
Definition Op_Conv_kschemas_centre_VEF.h:42
Op_Conv_kschemas_centre_VEF::ajouter
DoubleTab & ajouter(const DoubleTab &, DoubleTab &) const override
Definition Op_Conv_kschemas_centre_VEF.cpp:164
Op_Conv_kschemas_centre_VEF::convkschemas_centre
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
Definition Op_Conv_kschemas_centre_VEF.cpp:40
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
class Periodique This class represents a periodic boundary condition.
Definition Periodique.h:31
Periodique::face_associee
int face_associee(int i) const
Definition Periodique.h:35
Sortie
Base class for output streams.
Definition Sortie.h:52
TRUSTTab::nb_dim
int nb_dim() const
Definition TRUSTTab.h:199
TRUSTTab::resize
void resize(_SIZE_ n, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTTab.tpp:469
TRUSTTab::dimension
_SIZE_ dimension(int d) const
Definition TRUSTTab.tpp:133
TRUSTVect::echange_espace_virtuel
virtual void echange_espace_virtuel(IsExchangeBlocking exchange_type=IsExchangeBlocking::DefaultBlocking, const std::string kernel_name="noname")
Definition TRUSTVect.tpp:282