TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Terme_Source_inc_th_VDF_Face.cpp
1/****************************************************************************
2* Copyright (c) 2024, 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 <Terme_Source_inc_th_VDF_Face.h>
17#include <math.h>
18#include <Champ_Uniforme.h>
19#include <Domaine_VDF.h>
20#include <Domaine_Cl_VDF.h>
21#include <Neumann_sortie_libre.h>
22#include <Dirichlet.h>
23#include <Dirichlet_homogene.h>
24#include <Symetrie.h>
25#include <Periodique.h>
26#include <Pb_Fluide_base.h>
27#include <Convection_Diffusion_Temperature.h>
28#include <Navier_Stokes_std.h>
29#include <Fluide_base.h>
30#include <TRUSTTrav.h>
31
32Implemente_instanciable(Terme_Source_inc_th_VDF_Face,"inc_th_VDF_P0_VDF",Source_base);
33
34//// printOn
35//
36
37Sortie& Terme_Source_inc_th_VDF_Face::printOn(Sortie& s ) const
38{
39 return s << que_suis_je() ;
40}
41
42
43//// readOn
44//
45
46Entree& Terme_Source_inc_th_VDF_Face::readOn(Entree& s )
47{
48 return Terme_Source_inc_th::lire_donnees(s);
49}
50
51void Terme_Source_inc_th_VDF_Face::associer_pb(const Probleme_base& pb)
52{
53 if (sub_type(Pb_Fluide_base,pb) && (pb.nombre_d_equations()>1))
54 {
55 const Equation_base& eqn = pb.equation(1);
56 if (sub_type(Convection_Diffusion_Temperature,eqn))
57 {
58 const Navier_Stokes_std& eqn_hydr = ref_cast(Navier_Stokes_std,pb.equation(0));
59 const Convection_Diffusion_Temperature& eqn_thermo = ref_cast(Convection_Diffusion_Temperature,eqn);
60 Terme_Source_inc_th::associer_eqn(eqn_hydr);
61 Terme_Source_inc_th::associer_eqn(eqn_thermo);
62 }
63 }
64 else
65 {
66 Cerr << "Error for the method Terme_Source_inc_th_VDF_Face::associer_pb" << finl;
67 Cerr << "The source term "<<que_suis_je()<<" cannot be activated for the problem "<<pb.que_suis_je()<< finl;
68 exit();
69 }
70 const Fluide_base& fluide=eq_thermique().fluide();
71 beta_t_=fluide.beta_t();
72}
73
74
75void Terme_Source_inc_th_VDF_Face::associer_domaines(const Domaine_dis_base& domaine_dis,
76 const Domaine_Cl_dis_base& domaine_Cl_dis)
77{
78 le_dom_VDF = ref_cast(Domaine_VDF, domaine_dis);
79 le_dom_Cl_VDF = ref_cast(Domaine_Cl_VDF, domaine_Cl_dis);
80}
81
82// void Terme_Source_inc_th_VDF_Face::mettre_a_jour(double temps)
83// //void Turbulence_hyd_sous_maille_inc_VDF::mettre_a_jour(double )
84// {
85// int nb_elem_tot = le_dom_VDF->domaine().nb_elem_tot();
86// DoubleTrav cell_cent_vel(nb_elem_tot,dimension);
87// calculer_cell_cent_vel(cell_cent_vel);
88// }
89
90void Terme_Source_inc_th_VDF_Face::ajouter_blocs(matrices_t matrices, DoubleTab& secmem, const tabs_t& semi_impl) const
91{
92 const Domaine_VDF& domaine_VDF = le_dom_VDF.valeur();
93 const Domaine_Cl_VDF& domaine_Cl_VDF = le_dom_Cl_VDF.valeur();
94 const IntTab& face_voisins = domaine_VDF.face_voisins();
95 const IntVect& orientation = domaine_VDF.orientation();
96 //const DoubleTab& xv = domaine_VDF.xv();
97 const DoubleVect& porosite_surf = equation().milieu().porosite_face();
98 const DoubleVect& volumes_entrelaces = domaine_VDF.volumes_entrelaces();
99 const Domaine& domaine = le_dom_VDF->domaine();
100 int nb_elem = domaine.nb_elem();
101 const Champ_Uniforme& ch_beta = ref_cast(Champ_Uniforme,beta());
102 const DoubleTab& vitesse = eq_hydraulique().inconnue().valeurs();
103 const DoubleVect& temperature = eq_thermique().inconnue().valeurs();
104 DoubleVect cell_cent_temp(temperature);
105 int ndeb,nfin,ncomp,num_face;
106 // double debit_e=0.,dt,vol,Surf;
107 //double vol,Surf;
108 double vol;
109 DoubleTrav s(3);
110 s = 0.;
111 //int nb_elem_tot = le_dom_VDF->domaine().nb_elem_tot();
112 // double dQij=0;
113 // double si;
114
115 // static int deb=0;
116 // static double source=0.;
117 // static double debnm1;
118 // static double debit_ref=0.;
119 // Cerr << " debut de inc" << finl;
120 // double tps = eq_hydraulique().inconnue().temps();
121 // const DoubleVect& vit = eq_hydraulique().inconnue().valeurs();
122
123 // calculer_cell_cent_vel
124 // -------------------------------------------------------------------
125
126
127 // Exemple XTOF
128 // Cerr << " debut de parallel " << finl;
129 DoubleTab cell_cent_vel(0, dimension);
130 le_dom_VDF->domaine().creer_tableau_elements(cell_cent_vel);
131 // DoubleTrav toto(cell_cent_vel);
132
133 DoubleTab temp1(cell_cent_vel);
134 DoubleTab temp2(cell_cent_vel);
135 DoubleTab filt_vel(cell_cent_vel);
136 DoubleTab inc_dv(cell_cent_vel);
137
138 // DoubleTrav cell_cent_vel(nb_elem_tot,dimension);
139
140 // -------------------------------------------------------------------
141 // Cerr << " vitesse " << vitesse << finl;
142 // Cerr << " temperature " << cell_cent_temp << finl;
143 DoubleVect temp1_t(cell_cent_temp);
144 DoubleVect temp2_t(cell_cent_temp);
145 DoubleVect filt_temp(cell_cent_temp);
146 DoubleVect inc_dt(cell_cent_temp);
147 DoubleVect dQij_j(cell_cent_temp);
148
149 // const DoubleTab& vitesse = mon_equation->inconnue().valeurs();
150 // const Domaine_VDF& domaine_VDF = le_dom_VDF.valeur();
151 // int nb_elem_tot = domaine_VDF.domaine().nb_elem_tot();
152 const IntTab& elem_faces = domaine_VDF.elem_faces();
153 int element_number;
154 int num0=-1,num1,num2=-1,num3,num4=-1,num5=-1;
155 // Cerr << " vitesse " << vitesse << finl;
156 // This is to calculate the cell centered velocity values
157
158 DoubleTab uij_filt(cell_cent_vel);
159 // DoubleTab uij_filt(nb_elem,dimension,dimension);
160 for (element_number=0 ; element_number<nb_elem ; element_number ++)
161 {
162 // Cerr << me() << " " << element_number << finl;
163 cell_cent_temp(element_number)*=ch_beta.valeurs()(0,0);
164 temp1_t(element_number)=0.;
165 temp2_t(element_number)=0.;
166 filt_temp(element_number)=0.;
167 inc_dt(element_number)=0.;
168 dQij_j(element_number)=0.;
169 num0 = elem_faces(element_number,0);
170 num1 = elem_faces(element_number,1);
171 num2 = elem_faces(element_number,2);
172 num3 = elem_faces(element_number,3);
173 // uij_filt(element_number,0,0)=uij_filt_1(element_number,0);
174 // uij_filt(element_number,1,0)=uij_filt_1(element_number,1);
175 // uij_filt(element_number,0,1)=uij_filt_2(element_number,0);
176 // uij_filt(element_number,1,1)=uij_filt_2(element_number,1);
177 if (dimension == 3)
178 {
179 // uij_filt(element_number,0,2)=uij_filt_3(element_number,0);
180 // uij_filt(element_number,1,2)=uij_filt_3(element_number,1);
181 // uij_filt(element_number,2,2)=uij_filt_3(element_number,2);
182 // uij_filt(element_number,2,0)=uij_filt_1(element_number,2);
183 // uij_filt(element_number,2,1)=uij_filt_2(element_number,2);
184 num4 = elem_faces(element_number,4);
185 num5 = elem_faces(element_number,5);
186 inc_dv(element_number,0)=0.;
187 inc_dv(element_number,1)=0.;
188 inc_dv(element_number,2)=0.;
189 }
190
191 if (dimension == 2)
192 {
193 cell_cent_vel(element_number,0)=0.5*(vitesse[num0]+vitesse[num2]);
194 cell_cent_vel(element_number,1)=0.5*(vitesse[num1]+vitesse[num3]);
195 inc_dv(element_number,0)=0.;
196 inc_dv(element_number,1)=0.;
197 }
198 else
199 {
200 cell_cent_vel(element_number,0)=0.5*(vitesse[num0]+vitesse[num3]);
201 cell_cent_vel(element_number,1)=0.5*(vitesse[num1]+vitesse[num4]);
202 cell_cent_vel(element_number,2)=0.5*(vitesse[num2]+vitesse[num5]);
203 }
204 }
205 cell_cent_vel.echange_espace_virtuel();
206
207 // Cerr << " calculer_filter_field" << finl;
208 //
209 // calculer_filter_field(cell_cent_vel,filt_vel);
210 //
211 // calculate filter vector, unfilter tensor and
212 // increment fields in the same loop to speed up calculation
213 //
214 // -------------------------------------------------------------------
215
216 // DoubleTab toto(cell_cent_vel);
217
218 DoubleTab Qij(cell_cent_vel);
219 DoubleTab Lij(cell_cent_vel);
220 DoubleTab ttemp1(cell_cent_vel);
221 DoubleTab ttemp2(cell_cent_vel);
222 // DoubleTab Qij(uij_filt);
223 // DoubleTab Lij(uij_filt);
224 // DoubleTab ttemp1(uij_filt);
225 // DoubleTab ttemp2(uij_filt);
226 // DoubleTrav filt_vel(nb_elem_tot,dimension);
227 // DoubleTrav inc_dv(nb_elem_tot,dimension);
228 // inc_dv.resize(nb_elem_tot,dimension);
229 // DoubleTrav uij_filt(nb_elem_tot,dimension,dimension);
230
231 // -------------------------------------------------------------------
232 // const Domaine_VDF& domaine_VDF = le_dom_VDF.valeur();
233 // const IntTab& face_voisins = domaine_VDF.face_voisins();
234 // int nb_elem_tot = domaine_VDF.domaine().nb_elem_tot();
235 // const IntTab& elem_faces = domaine_VDF.elem_faces();
236 // int element_number;
237 // int num0,num1,num2,num3,num4,num5;
238 int f0,f1,f2,f3,f4,f5;
239 int i;
240 int num0_1,num1_1,num2_1,num3_1,num4_1,num5_1;
241 int f0_1,f1_1,f2_1,f3_1,f4_1,f5_1;
242
243 // DoubleTrav temp1(nb_elem_tot,dimension);
244 // DoubleTrav temp2(nb_elem_tot,dimension);
245
246 // This is to calculate the filtered velocity field
247 // Cerr << me() << " boucle 1 " << finl;
248 if (dimension ==2)
249 {
250 for (element_number=0 ; element_number<nb_elem ; element_number ++)
251 {
252 // for (i=0 ; i<dimension ; i++)
253 // for (j=0 ; j<dimension ; j++)
254 // uij_filt(element_number,i,j)=cell_cent_vel(element_number,i)*
255 // cell_cent_vel(element_number,j);
256 for (i=0 ; i<dimension ; i++)
257 {
258 uij_filt(element_number,i)=cell_cent_vel(element_number,i)*
259 cell_cent_temp(element_number);
260 }
261 f0 = elem_faces(element_number,0);
262 num0 = face_voisins(f0,0);
263 f2 = elem_faces(element_number,2);
264 num2 = face_voisins(f2,1);
265 if ( (num0 == -1) || (num2 == -1) )
266 {
267 for (i=0 ; i<dimension ; i ++)
268 {
269 temp1(element_number,i)=cell_cent_vel(element_number,i);
270 inc_dv(element_number,i)=0.;
271 }
272 temp1_t(element_number)=cell_cent_temp(element_number);
273 inc_dt(element_number)=0.;
274 }
275 else
276 {
277 {
278 for (i=0 ; i<dimension ; i ++)
279 {
280 temp1(element_number,i)=0.25*(cell_cent_vel(num0,i)+
281 2.0*cell_cent_vel(element_number,i)+
282 cell_cent_vel(num2,i));
283 }
284 temp1_t(element_number)=0.25*(cell_cent_temp(num0)+
285 2.0*cell_cent_temp(element_number)+
286 cell_cent_temp(num2));
287 }
288 f0_1 = elem_faces(num0,0);
289 num0_1 = face_voisins(f0_1,0);
290 f2_1 = elem_faces(num2,2);
291 num2_1 = face_voisins(f2_1,1);
292 if ( (num0_1 == -1) || (num2_1 == -1) )
293 {
294 for (i=0 ; i<dimension ; i ++)
295 {
296 inc_dv(element_number,i)+=(cell_cent_vel(num0,i))-(cell_cent_vel(num2,i));
297 }
298 inc_dt(element_number)+=(cell_cent_temp(num0))-(cell_cent_temp(num2));
299 }
300 else
301 {
302 for (i=0 ; i<dimension ; i ++)
303 {
304 inc_dv(element_number,i)+=(cell_cent_vel(num0_1,i))-(cell_cent_vel(num2_1,i));
305 }
306 inc_dt(element_number)+=(cell_cent_temp(num0_1))-(cell_cent_temp(num2_1));
307 }
308 }
309 }
310 temp1.echange_espace_virtuel();
311 temp1_t.echange_espace_virtuel();
312
313 for (element_number=0 ; element_number<nb_elem ; element_number ++)
314 {
315 f1 = elem_faces(element_number,1);
316 num1 = face_voisins(f1,0);
317 f3 = elem_faces(element_number,3);
318 num3 = face_voisins(f3,1);
319 if ( (num1 == -1) || (num3 == -1) )
320 {
321 for (i=0 ; i<dimension ; i ++)
322 {
323 filt_vel(element_number,i)=temp1(element_number,i);
324 inc_dv(element_number,i)=0.;
325 }
326 filt_temp(element_number)=temp1_t(element_number);
327 inc_dt(element_number)=0.;
328 }
329 else
330 {
331 for (i=0 ; i<dimension ; i ++)
332 {
333 filt_vel(element_number,i)=0.25*(temp1(num1,i)+
334 2.0*temp1(element_number,i)+
335 temp1(num3,i));
336 }
337 filt_temp(element_number)=0.25*(temp1_t(num1)+
338 2.0*temp1_t(element_number)+
339 temp1_t(num3));
340 f1_1 = elem_faces(num0,1);
341 num1_1 = face_voisins(f1_1,0);
342 f3_1 = elem_faces(num3,3);
343 num3_1 = face_voisins(f3_1,1);
344 if ( (num1_1 == -1) || (num3_1 == -1) )
345 {
346 for (i=0 ; i<dimension ; i ++)
347 inc_dv(element_number,i)+=(cell_cent_vel(num1,i))-(cell_cent_vel(num3,i));
348 inc_dt(element_number)+=(cell_cent_temp(num1))-(cell_cent_temp(num3));
349 }
350 else
351 {
352 for (i=0 ; i<dimension ; i ++)
353 {
354 inc_dv(element_number,i)+=(cell_cent_vel(num1_1,i))-(cell_cent_vel(num3_1,i));
355 }
356 inc_dt(element_number)+=(cell_cent_temp(num1_1))-(cell_cent_temp(num3_1));
357 }
358 }
359 }
360 }
361 else // dimension == 3
362 {
363 for (element_number=0 ; element_number<nb_elem ; element_number ++)
364 {
365 // for (i=0 ; i<dimension ; i++)
366 // for (j=0 ; j<dimension ; j++)
367 // uij_filt(element_number,i,j)=cell_cent_vel(element_number,i)*
368 // cell_cent_vel(element_number,j);
369 // Cerr << " elem " << element_number << finl;
370 for (i=0 ; i<dimension ; i++)
371 {
372 uij_filt(element_number,i)=cell_cent_vel(element_number,i)*
373 cell_cent_temp(element_number);
374 }
375 f0 = elem_faces(element_number,0);
376 num0 = face_voisins(f0,0);
377 f3 = elem_faces(element_number,3);
378 num3 = face_voisins(f3,1);
379 if ( (num0 == -1) || (num3 == -1) )
380 {
381 for (i=0 ; i<dimension ; i ++)
382 {
383 temp1(element_number,i)=cell_cent_vel(element_number,i);
384 inc_dv(element_number,i)+=0.;
385 }
386 temp1_t(element_number)=cell_cent_temp(element_number);
387 inc_dt(element_number)+=0.;
388 }
389 else
390 {
391 for (i=0 ; i<dimension ; i ++)
392 {
393 temp1(element_number,i)=0.25*(cell_cent_vel(num0,i)+
394 2.0*cell_cent_vel(element_number,i)+
395 cell_cent_vel(num3,i));
396 }
397 temp1_t(element_number)=0.25*(cell_cent_temp(num0)+
398 2.0*cell_cent_temp(element_number)+
399 cell_cent_temp(num3));
400 // Cerr << " not num0 et num3 " << num0 << " " << num3 << finl;
401 f0_1 = elem_faces(num0,0);
402 num0_1 = face_voisins(f0_1,0);
403 f3_1 = elem_faces(num3,3);
404 num3_1 = face_voisins(f3_1,1);
405 if ( (num0_1 == -1) || (num3_1 == -1) )
406 {
407 for (i=0 ; i<dimension ; i ++)
408 inc_dv(element_number,i)+=(cell_cent_vel(num0,i))-(cell_cent_vel(num3,i));
409 inc_dt(element_number)+=(cell_cent_temp(num0))-(cell_cent_temp(num3));
410 }
411 else
412 {
413 // Cerr << " not num0_1 et num3_1 " << num0_1 << " " << num3_1 << finl;
414 for (i=0 ; i<dimension ; i ++)
415 {
416 inc_dv(element_number,i)+=(cell_cent_vel(num0_1,i))-(cell_cent_vel(num3_1,i));
417 }
418 inc_dt(element_number)+=(cell_cent_temp(num0_1))-(cell_cent_temp(num3_1));
419 }
420 }
421 }
422
423 temp1.echange_espace_virtuel();
424 temp1_t.echange_espace_virtuel();
425 Cerr << " part 1 " << finl;
426 for (element_number=0 ; element_number<nb_elem ; element_number ++)
427 {
428 f1 = elem_faces(element_number,1);
429 num1 = face_voisins(f1,0);
430 f4 = elem_faces(element_number,4);
431 num4 = face_voisins(f4,1);
432 if ( (num1 == -1) || (num4 == -1) )
433 {
434 for (i=0 ; i<dimension ; i ++)
435 {
436 temp2(element_number,i)=temp1(element_number,i);
437 inc_dv(element_number,i)+=0.;
438 }
439 temp2_t(element_number)=temp1_t(element_number);
440 inc_dt(element_number)+=0.;
441 }
442 else
443 {
444 for (i=0 ; i<dimension ; i ++)
445 {
446 temp2(element_number,i)=0.25*(temp1(num1,i)+
447 2.0*temp1(element_number,i)+
448 temp1(num4,i));
449 }
450 temp2_t(element_number)=0.25*(temp1_t(num1)+
451 2.0*temp1_t(element_number)+
452 temp1_t(num4));
453 // Cerr << " not num1 et num4 " << num1 << " " << num4 << finl;
454 f1_1 = elem_faces(num1,1);
455 num1_1 = face_voisins(f1_1,0);
456 f4_1 = elem_faces(num4,4);
457 num4_1 = face_voisins(f4_1,1);
458 if ( (num1_1 == -1) || (num4_1 == -1) )
459 {
460 for (i=0 ; i<dimension ; i ++)
461 inc_dv(element_number,i)+=(cell_cent_vel(num1,i))-(cell_cent_vel(num4,i));
462 inc_dt(element_number)+=(cell_cent_temp(num1))-(cell_cent_temp(num4));
463 }
464 else
465 {
466 // Cerr << " not num1_1 et num4_1 " << num1_1 << " " << num4_1 << finl;
467 for (i=0 ; i<dimension ; i ++)
468 {
469 inc_dv(element_number,i)+=(cell_cent_vel(num1_1,i))-(cell_cent_vel(num4_1,i));
470 }
471 inc_dt(element_number)+=(cell_cent_temp(num1_1))-(cell_cent_temp(num4_1));
472 }
473 }
474 }
475
476 temp2.echange_espace_virtuel();
477 temp2_t.echange_espace_virtuel();
478 for (element_number=0 ; element_number<nb_elem ; element_number ++)
479 {
480 f2 = elem_faces(element_number,2);
481 num2 = face_voisins(f2,0);
482 f5 = elem_faces(element_number,5);
483 num5 = face_voisins(f5,1);
484 if ( (num2 == -1) || (num5 == -1) )
485 {
486 for (i=0 ; i<dimension ; i ++)
487 {
488 filt_vel(element_number,i)=temp2(element_number,i);
489 inc_dv(element_number,i)+=0.;
490 }
491 filt_temp(element_number)=temp2_t(element_number);
492 inc_dt(element_number)+=0.;
493 }
494 else
495 {
496 for (i=0 ; i<dimension ; i ++)
497 {
498 filt_vel(element_number,i)=0.25*(temp2(num2,i)+
499 2.0*temp2(element_number,i)+
500 temp2(num5,i));
501 }
502 filt_temp(element_number)=0.25*(temp2_t(num2)+
503 2.0*temp2_t(element_number)+
504 temp2_t(num5));
505 // Cerr << " not num2 et num5 " << num2 << " " << num5 << finl;
506 f2_1 = elem_faces(num2,2);
507 num2_1 = face_voisins(f2_1,0);
508 f5_1 = elem_faces(num5,5);
509 num5_1 = face_voisins(f5_1,1);
510 if ( (num2_1 == -1) || (num5_1 == -1) )
511 {
512 for (i=0 ; i<dimension ; i ++)
513 inc_dv(element_number,i)+=(cell_cent_vel(num2,i))-(cell_cent_vel(num5,i));
514 inc_dt(element_number)+=(cell_cent_temp(num2))-(cell_cent_temp(num5));
515 }
516 else
517 {
518 // Cerr << " not num2_1 et num5_1 " << num2_1 << " " << num5_1 << finl;
519 for (i=0 ; i<dimension ; i ++)
520 {
521 inc_dv(element_number,i)+=(cell_cent_vel(num2_1,i))-(cell_cent_vel(num5_1,i));
522 }
523 inc_dt(element_number)+=(cell_cent_temp(num2_1))-(cell_cent_temp(num5_1));
524 }
525 }
526 }
527 }
528 filt_vel.echange_espace_virtuel();
529 inc_dv.echange_espace_virtuel();
530 filt_temp.echange_espace_virtuel();
531 inc_dt.echange_espace_virtuel();
532 uij_filt.echange_espace_virtuel();
533 Cerr << " first loop " << finl;
534 //
535 // calculer_Qij(inc_dv,Qij);
536 // -------------------------------------------------------------------
537
538 // DoubleTrav Qij(nb_elem_tot,dimension,dimension);
539
540 // -------------------------------------------------------------------
541 // const Domaine_VDF& domaine_VDF = le_dom_VDF.valeur();
542 // int nb_elem_tot = domaine_VDF.domaine().nb_elem_tot();
543 // int element_number;
544
545 // This is to calculate the Qij tensor
546
547 // for (element_number=0 ; element_number<nb_elem_tot ; element_number ++)
548 // {
549 // for (int i=0 ; i<dimension ; i++)
550 // for (int j=0 ; j<dimension ; j++)
551 // Qij(element_number,i,j)=inc_dv(element_number,i)*inc_dv(element_number,j);
552 // }
553
554 // Cerr << " calculer_Lij " << finl;
555 //
556 // calculer_Lij(cell_cent_vel,filt_vel,Lij);
557 //
558 // calculate Qij and uij_filt in the same loop to speed up calculation
559 //
560 // -------------------------------------------------------------------
561
562 // DoubleTrav Lij(nb_elem_tot,dimension,dimension);
563
564 // -------------------------------------------------------------------
565 // const Domaine_VDF& domaine_VDF = le_dom_VDF.valeur();
566 // int nb_elem_tot = domaine_VDF.domaine().nb_elem_tot();
567 // int element_number;
568
569 // DoubleTrav uij_filt(nb_elem_tot,dimension,dimension);
570 // DoubleTrav ttemp1(nb_elem_tot,dimension,dimension);
571 // DoubleTrav ttemp2(nb_elem_tot,dimension,dimension);
572
573 // This is to calculate the Lij term for the C coefficient
574
575 // for (element_number=0 ; element_number<nb_elem_tot ; element_number ++)
576 // {
577 // for (int i=0 ; i<dimension ; i++)
578 // for (int j=0 ; j<dimension ; j++)
579 // {
580 // uij_filt(element_number,i,j)=cell_cent_vel(element_number,i)*
581 // cell_cent_vel(element_number,j);
582 // }
583 // }
584
585 // calculer_filter_tensor(uij_filt);
586 // Cerr << " calculer_uij_filt " << finl;
587
588 // This is to calculate the filtered tensor field
589 // int j;
590 // Cerr << me() << " boucle 2 " << finl;
591 if (dimension ==2)
592 {
593 for (element_number=0 ; element_number<nb_elem ; element_number ++)
594 {
595 f0 = elem_faces(element_number,0);
596 num0 = face_voisins(f0,0);
597 f2 = elem_faces(element_number,2);
598 num2 = face_voisins(f2,1);
599 if ( (num0 == -1) || (num2 == -1) )
600 // for (i=0 ; i<dimension ; i ++)
601 // for (j=0 ; j<dimension ; j ++)
602 // ttemp1(element_number,i,j)=uij_filt(element_number,i,j);
603 for (i=0 ; i<dimension ; i++)
604 {
605 ttemp1(element_number,i)=uij_filt(element_number,i);
606 }
607 else
608 {
609 // for (i=0 ; i<dimension ; i ++)
610 // for (j=0 ; j<dimension ; j ++)
611 // {
612 // ttemp1(element_number,i,j)=0.25*(uij_filt(num0,i,j)+
613 // 2.0*uij_filt(element_number,i,j)+
614 // uij_filt(num2,i,j));
615 // }
616 for (i=0 ; i<dimension ; i++)
617 {
618 ttemp1(element_number,i)=0.25*(uij_filt(num0,i)+
619 2.0*uij_filt(element_number,i)+
620 uij_filt(num2,i));
621 }
622 }
623 }
624
625 ttemp1.echange_espace_virtuel();
626 for (element_number=0 ; element_number<nb_elem ; element_number ++)
627 {
628 f1 = elem_faces(element_number,1);
629 num1 = face_voisins(f1,0);
630 f3 = elem_faces(element_number,3);
631 num3 = face_voisins(f3,1);
632 if ( (num1 == -1) || (num3 == -1) )
633 // for (i=0 ; i<dimension ; i ++)
634 // for (j=0 ; j<dimension ; j ++)
635 // uij_filt(element_number,i,j)=ttemp1(element_number,i,j);
636 for (i=0 ; i<dimension ; i++)
637 {
638 uij_filt(element_number,i)=ttemp1(element_number,i);
639 }
640 else
641 {
642 // for (i=0 ; i<dimension ; i ++)
643 // for (j=0 ; j<dimension ; j ++)
644 // {
645 // uij_filt(element_number,i,j)=0.25*(ttemp1(num1,i,j)+
646 // 2.0*ttemp1(element_number,i,j)+
647 // ttemp1(num3,i,j));
648 // }
649 for (i=0 ; i<dimension ; i++)
650 {
651 uij_filt(element_number,i)=0.25*(ttemp1(num0,i)+
652 2.0*ttemp1(element_number,i)+
653 ttemp1(num2,i));
654 }
655 }
656 }
657 }
658 else
659 {
660 for (element_number=0 ; element_number<nb_elem ; element_number ++)
661 {
662 f0 = elem_faces(element_number,0);
663 num0 = face_voisins(f0,0);
664 f3 = elem_faces(element_number,3);
665 num3 = face_voisins(f3,1);
666 if ( (num0 == -1) || (num3 == -1) )
667 // for (i=0 ; i<dimension ; i ++)
668 // for (j=0 ; j<dimension ; j ++)
669 // ttemp1(element_number,i,j)=uij_filt(element_number,i,j);
670 for (i=0 ; i<dimension ; i++)
671 {
672 ttemp1(element_number,i)=uij_filt(element_number,i);
673 }
674 else
675 {
676 // for (i=0 ; i<dimension ; i ++)
677 // for (j=0 ; j<dimension ; j ++)
678 // {
679 // ttemp1(element_number,i,j)=0.25*(uij_filt(num0,i,j)+
680 // 2.0*uij_filt(element_number,i,j)+
681 // uij_filt(num3,i,j));
682 // }
683 for (i=0 ; i<dimension ; i++)
684 {
685 ttemp1(element_number,i)=0.25*(uij_filt(num0,i)+
686 2.0*uij_filt(element_number,i)+
687 uij_filt(num2,i));
688 }
689 }
690 }
691
692 ttemp1.echange_espace_virtuel();
693 for (element_number=0 ; element_number<nb_elem ; element_number ++)
694 {
695 f1 = elem_faces(element_number,1);
696 num1 = face_voisins(f1,0);
697 f4 = elem_faces(element_number,4);
698 num4 = face_voisins(f4,1);
699 if ( (num1 == -1) || (num4 == -1) )
700 // for (i=0 ; i<dimension ; i ++)
701 // for (j=0 ; j<dimension ; j ++)
702 // ttemp2(element_number,i,j)=ttemp1(element_number,i,j);
703 for (i=0 ; i<dimension ; i++)
704 {
705 ttemp2(element_number,i)=ttemp1(element_number,i);
706 }
707 else
708 {
709 // for (i=0 ; i<dimension ; i ++)
710 // for (j=0 ; j<dimension ; j ++)
711 // {
712 // ttemp2(element_number,i,j)=0.25*(ttemp1(num1,i,j)+
713 // 2.0*ttemp1(element_number,i,j)+
714 // ttemp1(num4,i,j));
715 // }
716 }
717 }
718
719 ttemp2.echange_espace_virtuel();
720 for (element_number=0 ; element_number<nb_elem ; element_number ++)
721 {
722 f2 = elem_faces(element_number,2);
723 num2 = face_voisins(f2,0);
724 f5 = elem_faces(element_number,5);
725 num5 = face_voisins(f5,1);
726 if ( (num2 == -1) || (num5 == -1) )
727 // for (i=0 ; i<dimension ; i ++)
728 // for (j=0 ; j<dimension ; j ++)
729 // uij_filt(element_number,i,j)=ttemp2(element_number,i,j);
730 for (i=0 ; i<dimension ; i++)
731 {
732 uij_filt(element_number,i)=ttemp2(element_number,i);
733 }
734 else
735 {
736 // for (i=0 ; i<dimension ; i ++)
737 // for (j=0 ; j<dimension ; j ++)
738 // {
739 // uij_filt(element_number,i,j)=0.25*(ttemp2(num2,i,j)+
740 // 2.0*ttemp2(element_number,i,j)+
741 // ttemp2(num5,i,j));
742 // }
743 for (i=0 ; i<dimension ; i++)
744 {
745 uij_filt(element_number,i)=0.25*(ttemp2(num0,i)+
746 2.0*ttemp2(element_number,i)+
747 ttemp2(num2,i));
748 }
749 }
750 }
751 }
752 // finish calculating Lij
753
754 uij_filt.echange_espace_virtuel();
755 // Cerr << me() << " boucle 3 " << finl;
756 for (element_number=0 ; element_number<nb_elem ; element_number ++)
757 {
758 // for (int i=0 ; i<dimension ; i++)
759 // for (int j=0 ; j<dimension ; j++)
760 // {
761 // Qij(element_number,i,j)=inc_dv(element_number,i)*inc_dv(element_number,j);
762 // Lij(element_number,i,j)=uij_filt(element_number,i,j)-
763 // filt_vel(element_number,i)*filt_vel(element_number,j);
764 // }
765 for (int ii=0 ; ii<dimension ; ii++)
766 {
767 Qij(element_number,ii)=inc_dv(element_number,ii)*inc_dt(element_number);
768 Lij(element_number,ii)=uij_filt(element_number,ii)-
769 filt_vel(element_number,ii)*filt_temp(element_number);
770 }
771 }
772
773 Lij.echange_espace_virtuel();
774 Qij.echange_espace_virtuel();
775 //
776 // calculer_f(Lij,Qij);
777 // -------------------------------------------------------------------
778
779 DoubleVect f(nb_elem);
780
781 // -------------------------------------------------------------------
782 DoubleVect temp_1(nb_elem);
783 DoubleVect temp_2(nb_elem);
784
785 // Evaluate the dynamic model coefficient f
786 int elem;
787 if (dimension == 2)
788 for (elem=0; elem<nb_elem; elem++)
789 {
790 // temp_1(elem)=0.;
791 // temp_2(elem)=0.;
792 // for (int k=0 ; k<dimension ; k++)
793 // for (int l=0 ; l<dimension ; l++)
794 // {
795 // temp_1(elem)+=Lij(elem,k,k);
796 // temp_2(elem)+=4.*Qij(elem,l,l);
797 // }
798
799 temp_1(elem)=Lij(elem,0)+Lij(elem,1);
800 temp_2(elem)=4.*Qij(elem,0)+4.*Qij(elem,1);
801 if(std::fabs(temp_2(elem)) < 1.e-12)
802 f(elem)=0.;
803 else
804 f(elem)=temp_1(elem)/temp_2(elem);
805 }
806 else
807 for (elem=0; elem<nb_elem; elem++)
808 {
809 // temp_1(elem)=0.;
810 // temp_2(elem)=0.;
811
812 temp_1(elem)=Lij(elem,0)+Lij(elem,1)+Lij(elem,2);
813 temp_2(elem)=4.*Qij(elem,0)+4.*Qij(elem,1)+4.*Qij(elem,2);
814 if(std::fabs(temp_2(elem)) < 1.e-12)
815 f(elem)=0.;
816 else
817 f(elem)=temp_1(elem)/temp_2(elem);
818 }
819
820 f.echange_espace_virtuel();
821 // Cerr << " calculer_dQij_j" << f << " " << Qij << finl;
822 //
823 // calculer_dQij_j(f,Qij,dQij_j);
824 // -------------------------------------------------------------------
825
826 // DoubleTrav dQij_j(nb_elem_tot,dimension);
827 // dQij_j.resize(nb_elem_tot,dimension);
828
829 // -------------------------------------------------------------------
830 //const IntTab& Qdm = domaine_VDF.Qdm();
831 // const IntVect& orientation = domaine_VDF.orientation();
832
833 //
834 // Calculate the source term dQij_j = d/dxj ( f Qij )
835 //
836 // Cerr << me() << " boucle 4 " << nb_elem << finl;
837 if (dimension ==2)
838 {
839 for (element_number=0 ; element_number<nb_elem ; element_number ++)
840 {
841 f0 = elem_faces(element_number,0);
842 num0 = face_voisins(f0,0);
843 if (num0 == -1)
844 num0 = element_number;
845 f1 = elem_faces(element_number,1);
846 num1 = face_voisins(f1,0);
847 if (num1 == -1)
848 num1 = element_number;
849 f2 = elem_faces(element_number,2);
850 num2 = face_voisins(f2,1);
851 if (num2 == -1)
852 num2 = element_number;
853 f3 = elem_faces(element_number,3);
854 num3 = face_voisins(f3,1);
855 if (num3 == -1)
856 num3 = element_number;
857 dQij_j(element_number)=0.5*((f(num2)*Qij(num2,0)-f(num0)*Qij(num0,0))/
858 domaine_VDF.dim_elem(element_number,0));
859 // dQij_j(element_number,0)+=0.5*((f(num3)*Qij(num3,0)-f(num1)*Qij(num1,0))/
860 // domaine_VDF.dim_elem(element_number,1));
861 dQij_j(element_number)+=0.5*((f(num2)*Qij(num2,1)-f(num0)*Qij(num0,1))/
862 domaine_VDF.dim_elem(element_number,0));
863 // dQij_j(element_number,1)+=0.5*((f(num3)*Qij(num3,1)-f(num1)*Qij(num1,1))/
864 // domaine_VDF.dim_elem(element_number,1));
865 // dQij_j(element_number,0)=0.5*((f(num2)*Qij_1(num2,0)-f(num0)*Qij_1(num0,0))/
866 // domaine_VDF.dim_elem(element_number,0));
867 // dQij_j(element_number,0)+=0.5*((f(num3)*Qij_2(num3,0)-f(num1)*Qij_2(num1,0))/
868 // domaine_VDF.dim_elem(element_number,1));
869 // dQij_j(element_number,1)=0.5*((f(num2)*Qij_1(num2,1)-f(num0)*Qij_1(num0,1))/
870 // domaine_VDF.dim_elem(element_number,0));
871 // dQij_j(element_number,1)+=0.5*((f(num3)*Qij_2(num3,1)-f(num1)*Qij_2(num1,1))/
872 // domaine_VDF.dim_elem(element_number,1));
873 }
874 }
875 else
876 {
877 for (element_number=0 ; element_number<nb_elem ; element_number ++)
878 {
879 f0 = elem_faces(element_number,0);
880 num0 = face_voisins(f0,0);
881 if (num0 == -1 || num0 >= nb_elem)
882 {
883 // Cerr << me() << " bord num0 " << finl;
884 num0 = element_number;
885 }
886 f1 = elem_faces(element_number,1);
887 num1 = face_voisins(f1,0);
888 if (num1 == -1 || num0 >= nb_elem)
889 {
890 // Cerr << me() << " bord num1 " << finl;
891 num1 = element_number;
892 }
893 f2 = elem_faces(element_number,2);
894 num2 = face_voisins(f2,0);
895 if (num2 == -1 || num2 >= nb_elem)
896 {
897 // Cerr << me() << " bord num2 " << finl;
898 num2 = element_number;
899 }
900 f3 = elem_faces(element_number,3);
901 num3 = face_voisins(f3,1);
902 if (num3 == -1 || num3 >= nb_elem)
903 {
904 // Cerr << me() << " bord num3 " << finl;
905 num3 = element_number;
906 }
907 f4 = elem_faces(element_number,4);
908 num4 = face_voisins(f4,1);
909 if (num4 == -1 || num4 >= nb_elem)
910 {
911 // Cerr << me() << " bord num4 " << finl;
912 num4 = element_number;
913 }
914 f5 = elem_faces(element_number,5);
915 num5 = face_voisins(f5,1);
916 if (num5 == -1 || num5 >= nb_elem)
917 {
918 // Cerr << me() << " bord num5 " << finl;
919 num5 = element_number;
920 }
921 dQij_j(element_number)=0.5*((f(num3)*Qij(num3,0)-f(num0)*Qij(num0,0))/
922 domaine_VDF.dim_elem(element_number,0));
923 // dQij_j(element_number,0)+=0.5*((f(num4)*Qij(num4,0)-f(num1)*Qij(num1,0))/
924 // domaine_VDF.dim_elem(element_number,1));
925 // dQij_j(element_number,0)+=0.5*((f(num5)*Qij(num5,0)-f(num2)*Qij(num2,0))/
926 // domaine_VDF.dim_elem(element_number,2));
927 // dQij_j(element_number,0)=0.5*((f(num3)*Qij_1(num3,0)-f(num0)*Qij_1(num0,0))/
928 // domaine_VDF.dim_elem(element_number,0));
929 // dQij_j(element_number,0)+=0.5*((f(num4)*Qij_2(num4,0)-f(num1)*Qij_2(num1,0))/
930 // domaine_VDF.dim_elem(element_number,1));
931 // dQij_j(element_number,0)+=0.5*((f(num5)*Qij_3(num5,0)-f(num2)*Qij_3(num2,0))/
932 // domaine_VDF.dim_elem(element_number,2));
933
934
935 // dQij_j(element_number,1)=0.5*((f(num3)*Qij(num3,1)-f(num0)*Qij(num0,1))/
936 // domaine_VDF.dim_elem(element_number,0));
937 dQij_j(element_number)+=0.5*((f(num4)*Qij(num4,1)-f(num1)*Qij(num1,1))/
938 domaine_VDF.dim_elem(element_number,1));
939 // dQij_j(element_number,1)+=0.5*((f(num5)*Qij(num5,1)-f(num2)*Qij(num2,1))/
940 // domaine_VDF.dim_elem(element_number,2));
941 // dQij_j(element_number,1)=0.5*((f(num3)*Qij_1(num3,1)-f(num0)*Qij_1(num0,1))/
942 // domaine_VDF.dim_elem(element_number,0));
943 // dQij_j(element_number,1)+=0.5*((f(num4)*Qij_2(num4,1)-f(num1)*Qij_2(num1,1))/
944 // domaine_VDF.dim_elem(element_number,1));
945 // dQij_j(element_number,1)+=0.5*((f(num5)*Qij_3(num5,1)-f(num2)*Qij_3(num2,1))/
946 // domaine_VDF.dim_elem(element_number,2));
947
948 // dQij_j(element_number,2)=0.5*((f(num3)*Qij(num3,2)-f(num0)*Qij(num0,2))/
949 // domaine_VDF.dim_elem(element_number,0));
950 // dQij_j(element_number,2)+=0.5*((f(num4)*Qij(num4,2)-f(num1)*Qij(num1,2))/
951 // domaine_VDF.dim_elem(element_number,1));
952 dQij_j(element_number)+=0.5*((f(num5)*Qij(num5,2)-f(num2)*Qij(num2,2))/
953 domaine_VDF.dim_elem(element_number,2));
954 // dQij_j(element_number,2)=0.5*((f(num3)*Qij_1(num3,2)-f(num0)*Qij_1(num0,2))/
955 // domaine_VDF.dim_elem(element_number,0));
956 // dQij_j(element_number,2)+=0.5*((f(num4)*Qij_2(num4,2)-f(num1)*Qij_2(num1,2))/
957 // domaine_VDF.dim_elem(element_number,1));
958 // dQij_j(element_number,2)+=0.5*((f(num5)*Qij_3(num5,2)-f(num2)*Qij_3(num2,2))/
959 // domaine_VDF.dim_elem(element_number,2));
960 }
961 }
962 dQij_j.echange_espace_virtuel();
963 double min_f=1.e6, max_f=-1.e6;
964 double min_q=1.e6, max_q=-1.e6;
965 double min_dq=1.e6, max_dq=-1.e6;
966 double min_dq1=1.e6, max_dq1=-1.e6, sum_t=0.;
967 for (element_number=0 ; element_number<nb_elem ; element_number ++)
968 {
969 if (inc_dv(element_number,0)<min_q)
970 min_q=inc_dv(element_number,0);
971 if (inc_dv(element_number,0)>max_q)
972 max_q=inc_dv(element_number,0);
973 if (inc_dt(element_number)<min_f)
974 min_f=inc_dt(element_number);
975 if (inc_dt(element_number)>max_f)
976 max_f=inc_dt(element_number);
977 if (inc_dv(element_number,1)<min_dq)
978 min_dq=inc_dv(element_number,1);
979 if (inc_dv(element_number,1)>max_dq)
980 max_dq=inc_dv(element_number,1);
981 if (inc_dv(element_number,2)<min_dq1)
982 min_dq1=inc_dv(element_number,2);
983 if (inc_dv(element_number,2)>max_dq1)
984 max_dq1=inc_dv(element_number,2);
985 sum_t+=inc_dt(element_number);
986 }
987 Cerr << me() << " min/max " << min_f << " " << max_f << " " << min_dq1 << " " << max_dq1 << finl;
988 Cerr << " min/max " << (sum_t/nb_elem) << " " << min_q << " " << max_q << " " << min_dq << " " << max_dq << finl;
989
990 double inc_s=0;
991 // Cerr << "Force de volume (gradient de pression moteur) = " << source << finl;
992
993 // Boucle sur les conditions limites pour traiter les faces de bord
994 // Cerr << me() << " cond lims " << finl;
995 for (int n_bord=0; n_bord<domaine_VDF.nb_front_Cl(); n_bord++)
996 {
997
998 // pour chaque Condition Limite on regarde son type
999 // Si face de Dirichlet ou de Symetrie on ne fait rien
1000 // Si face de Neumann on calcule la contribution au terme source
1001
1002 const Cond_lim& la_cl = domaine_Cl_VDF.les_conditions_limites(n_bord);
1003
1004 if (sub_type(Neumann_sortie_libre,la_cl.valeur()))
1005 {
1006
1007 //const Neumann_sortie_libre& la_cl_neumann = ref_cast(Neumann_sortie_libre,la_cl.valeur());
1008 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1009 ndeb = le_bord.num_premiere_face();
1010 nfin = ndeb + le_bord.nb_faces();
1011
1012 for (num_face=ndeb; num_face<nfin; num_face++)
1013 {
1014 vol = volumes_entrelaces(num_face)*porosite_surf(num_face);
1015 ncomp = orientation(num_face);
1016 num0 = face_voisins(num_face,0);
1017 if (num0 == -1)
1018 num0 = face_voisins(num_face,1);
1019 inc_s=dQij_j(num0)*vol;
1020 inc_s=0.;
1021 // resu(num_face)+= s(ncomp)*vol;
1022 secmem(num_face)+= inc_s;
1023 }
1024 }
1025
1026 else if (sub_type(Symetrie,la_cl.valeur()))
1027 { /* Do nothing */}
1028 else if ( (sub_type(Dirichlet,la_cl.valeur()))
1029 ||
1030 (sub_type(Dirichlet_homogene,la_cl.valeur()))
1031 )
1032 { /* Do nothing */}
1033 else if (sub_type(Periodique,la_cl.valeur()))
1034 {
1035 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1036 ndeb = le_bord.num_premiere_face();
1037 nfin = ndeb + le_bord.nb_faces();
1038
1039 for (num_face=ndeb; num_face<nfin; num_face++)
1040 {
1041 vol = volumes_entrelaces(num_face)*porosite_surf(num_face);
1042 ncomp = orientation(num_face);
1043 num0 = face_voisins(num_face,0);
1044 if (num0 == -1)
1045 num0 = face_voisins(num_face,1);
1046 inc_s=dQij_j(num0)*vol;
1047 inc_s=0.;
1048 // resu(num_face)+= s(ncomp)*vol;
1049 secmem(num_face)+= inc_s;
1050
1051 }
1052 }
1053 }
1054
1055 // Boucle sur les faces internes
1056 Cerr << me() << " faces internes thermique" << ch_beta.valeurs()(0,0) << finl;
1057
1058 double min_inc=1.e6, max_inc=-1.e6;
1059 int num_e1=-1, n_comp1=-1;
1060 int num_e2=-1, n_comp2=-1;
1061 ndeb = domaine_VDF.premiere_face_int();
1062 for (num_face =domaine_VDF.premiere_face_int(); num_face<domaine_VDF.nb_faces(); num_face++)
1063 {
1064
1065 vol = volumes_entrelaces(num_face)*porosite_surf(num_face);
1066 ncomp = orientation(num_face);
1067 num0 = face_voisins(num_face,0);
1068 num1 = face_voisins(num_face,1);
1069 inc_s=0.5*(dQij_j(num0)+dQij_j(num1))*vol;
1070 // resu(num_face)+= s(ncomp)*vol;
1071 if (inc_s<min_inc)
1072 {
1073 min_inc=inc_s;
1074 num_e1=num0;
1075 n_comp1=ncomp;
1076 }
1077 if (inc_s>max_inc)
1078 {
1079 max_inc=inc_s;
1080 num_e2=num0;
1081 n_comp2=ncomp;
1082 }
1083 }
1084 // inc_s=0.;
1085 secmem(num_face)+= inc_s;
1086 Cerr << min_inc << " " << num_e1 << " " << n_comp1 << finl;
1087 Cerr << max_inc << " " << num_e2 << " " << n_comp2 << finl;
1088 secmem.echange_espace_virtuel();
1089}
1090
1091DoubleTab& Terme_Source_inc_th_VDF_Face::calculer(DoubleTab& resu) const
1092{
1093 resu = 0;
1094 ajouter(resu);
1095 resu.echange_espace_virtuel();
1096 return resu;
1097}
1098
1099
Champ_Don_base::valeurs
DoubleTab & valeurs() override
Surcharge Champ_base::valeurs() Renvoie le tableau des valeurs.
Definition Champ_Don_base.h:49
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
Definition Champ_Inc_base.h:77
Champ_Uniforme
classe Champ_Uniforme Represente un champ constant dans l'espace et dans le temps.
Definition Champ_Uniforme.h:26
Cond_lim
classe Cond_lim Classe generique servant a representer n'importe quelle classe
Definition Cond_lim.h:31
Convection_Diffusion_Temperature_base::fluide
const Fluide_base & fluide() const
Definition Convection_Diffusion_Temperature_base.cpp:26
Convection_Diffusion_Temperature
classe Convection_Diffusion_Temperature Cas particulier de Convection_Diffusion_std
Definition Convection_Diffusion_Temperature.h:28
Convection_Diffusion_Temperature::inconnue
const Champ_Inc_base & inconnue() const override
Definition Convection_Diffusion_Temperature.h:34
Dirichlet_homogene
Classe Dirichlet_homogene Cette classe est la classe de base de la hierarchie des conditions aux limi...
Definition Dirichlet_homogene.h:31
Dirichlet
classe Dirichlet Cette classe est la classe de base de la hierarchie des conditions aux limites de ty...
Definition Dirichlet.h:31
Domaine_Cl_VDF
class Domaine_Cl_VDF
Definition Domaine_Cl_VDF.h:63
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_VDF
class Domaine_VDF
Definition Domaine_VDF.h:64
Domaine_VDF::dim_elem
double dim_elem(int, int) const
Definition Domaine_VDF.h:627
Domaine_VDF::orientation
int orientation(int) const override
inline DoubleVect& Domaine_VDF::porosite_face() {
Definition Domaine_VDF.h:297
Domaine_VF::nb_faces
int nb_faces() const
renvoie le nombre global de faces.
Definition Domaine_VF.h:471
Domaine_VF::volumes_entrelaces
DoubleVect & volumes_entrelaces()
Definition Domaine_VF.h:99
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::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_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_front_Cl
int nb_front_Cl() const
Definition Domaine_dis_base.h:53
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Equation_base
classe Equation_base Le role d'une equation est le calcul d'un ou plusieurs champs....
Definition Equation_base.h:76
Equation_base::milieu
virtual const Milieu_base & milieu() const =0
Fluide_base
classe Fluide_base Cette classe represente un d'un fluide incompressible ainsi que
Definition Fluide_base.h:38
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::beta_t
virtual const Champ_Don_base & beta_t() const
Renvoie beta_t du milieu.
Definition Milieu_base.cpp:887
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
Navier_Stokes_std
classe Navier_Stokes_std Cette classe porte les termes de l'equation de la dynamique
Definition Navier_Stokes_std.h:56
Navier_Stokes_std::inconnue
const Champ_Inc_base & inconnue() const override
Renvoie la vitesse (champ inconnue de l'equation) (version const).
Definition Navier_Stokes_std.cpp:599
Neumann_sortie_libre
classe Neumann_sortie_libre Cette classe represente une frontiere ouverte sans vitesse imposee
Definition Neumann_sortie_libre.h:34
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
Pb_Fluide_base
classe Pb_Fluide_base Cette classe a pour but de disposer d une classe amont pour
Definition Pb_Fluide_base.h:29
Periodique
classe Periodique Cette classe represente une condition aux limites periodique.
Definition Periodique.h:31
Probleme_base
classe Probleme_base C'est un Probleme_U qui n'est pas un couplage.
Definition Probleme_base.h:55
Probleme_base::nombre_d_equations
virtual int nombre_d_equations() const =0
Probleme_base::equation
virtual const Equation_base & equation(int) const =0
Process::me
static int me()
renvoie mon rang dans le groupe de communication courant.
Definition Process.cpp:125
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
Source_base
classe Source_base Un objet Source_base est un terme apparaissant au second membre d'une
Definition Source_base.h:42
Source_base::ajouter
virtual DoubleTab & ajouter(DoubleTab &) const
Definition Source_base.cpp:225
Symetrie
classe Symetrie Sur les faces de symetrie on a les proprietes suivantes:
Definition Symetrie.h:37
TRUSTVect::echange_espace_virtuel
virtual void echange_espace_virtuel(IsExchangeBlocking exchange_type=IsExchangeBlocking::DefaultBlocking, const std::string kernel_name="noname")
Definition TRUSTVect.tpp:282
Terme_Source_inc_th_VDF_Face
class Terme_Source_inc_th_VDF_Face Cette classe permet de conserver le debit dans une simulation
Definition Terme_Source_inc_th_VDF_Face.h:39
Terme_Source_inc_th_VDF_Face::beta
const Champ_Don_base & beta() const
Definition Terme_Source_inc_th_VDF_Face.h:73
Terme_Source_inc_th_VDF_Face::ajouter_blocs
void ajouter_blocs(matrices_t matrices, DoubleTab &secmem, const tabs_t &semi_impl) const override
Definition Terme_Source_inc_th_VDF_Face.cpp:90
Terme_Source_inc_th_VDF_Face::associer_pb
void associer_pb(const Probleme_base &) override
Definition Terme_Source_inc_th_VDF_Face.cpp:51
Terme_Source_inc_th_VDF_Face::associer_domaines
void associer_domaines(const Domaine_dis_base &, const Domaine_Cl_dis_base &) override
Definition Terme_Source_inc_th_VDF_Face.cpp:75
Terme_Source_inc_th_VDF_Face::calculer
DoubleTab & calculer(DoubleTab &) const override
Definition Terme_Source_inc_th_VDF_Face.cpp:1091
Terme_Source_inc_th::eq_hydraulique
const Navier_Stokes_std & eq_hydraulique() const
Definition Terme_Source_inc_th.h:46
Terme_Source_inc_th::lire_donnees
Entree & lire_donnees(Entree &)
Definition Terme_Source_inc_th.cpp:31
Terme_Source_inc_th::associer_eqn
void associer_eqn(const Navier_Stokes_std &)
Definition Terme_Source_inc_th.cpp:20
Terme_Source_inc_th::eq_thermique
const Convection_Diffusion_Temperature & eq_thermique() const
Definition Terme_Source_inc_th.h:51