TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Modele_Jones_Launder_VDF.cpp
1/****************************************************************************
2* Copyright (c) 2015 - 2016, 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 <Modele_Jones_Launder_VDF.h>
17#include <Domaine_VDF.h>
18#include <Domaine_Cl_VDF.h>
19#include <TRUSTTrav.h>
20#include <Dirichlet.h>
21#include <Dirichlet_homogene.h>
22#include <Dirichlet_paroi_defilante.h>
23#include <Echange_externe_impose.h>
24#include <Periodique.h>
25#include <Symetrie.h>
26#include <Neumann.h>
27#include <Neumann_homogene.h>
28#include <Champ_Face_VDF.h>
29#include <Champ_Uniforme.h>
30#include <Milieu_base.h>
31
32Implemente_instanciable(Modele_Jones_Launder_VDF,"Modele_Jones_Launder_VDF",Modele_Fonc_Bas_Reynolds_Base);
33// XD Jones_Launder modele_fonction_bas_reynolds_base Jones_Launder INHERITS_BRACE Model described in ' Jones, W. P. and
34// XD_CONT Launder, B. E. (1972), The prediction of laminarization with a two-equation model of turbulence, Int. J. of
35// XD_CONT Heat and Mass transfer, Vol. 15, pp. 301-314.'
36// printOn et readOn
37
38Sortie& Modele_Jones_Launder_VDF::printOn(Sortie& s ) const
39{
40 return s;
41}
42
43Entree& Modele_Jones_Launder_VDF::readOn(Entree& is )
44{
45 Motcle motlu, accolade_fermee="}", accolade_ouverte="{";
46 is >> motlu;
47 if (motlu==accolade_ouverte)
48 {
49 is >> motlu;
50 if (motlu != accolade_fermee)
51 {
52 Cerr << "Erreur a la lecture du Modele fonc bas reynolds Jones et Launder" << finl;
53 Cerr << "On attendait } a la place de " << motlu << finl;
54 Process::exit();
55 }
56 }
57 else
58 {
59 Cerr << "Erreur a la lecture du Modele fonc bas reynolds Jones et Launder" << finl;
60 Cerr << "On attendait { a la place de " << motlu << finl;
61 Process::exit();
62 }
63 return is;
64}
65
66Entree& Modele_Jones_Launder_VDF::lire(const Motcle& , Entree& is)
67{
68 return is;
69}
70///////////////////////////////////////////////////////////////
71// Implementation des fonctions de la classe
72///////////////////////////////////////////////////////////////
73
74void Modele_Jones_Launder_VDF::associer(const Domaine_dis_base& domaine_dis,
75 const Domaine_Cl_dis_base& domaine_Cl_dis)
76{
77 // const Domaine_VDF& le_dom = ref_cast(Domaine_VDF,domaine_dis);
78 // const Domaine_Cl_VDF& le_dom_Cl = ref_cast(Domaine_Cl_VDF,domaine_Cl_dis);
79}
80
81
82
83DoubleTab& Modele_Jones_Launder_VDF::Calcul_D(DoubleTab& D,const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis,
84 const DoubleTab& vitesse,const DoubleTab& K_eps_Bas_Re, const Champ_Don_base& ch_visco ) const
85{
86 double visco=-1;
87 const DoubleTab& tab_visco=ch_visco.valeurs();
88 int is_visco_const=sub_type(Champ_Uniforme,ch_visco);
89 if (is_visco_const)
90 visco=tab_visco(0,0);
91
92 const Domaine_VDF& le_dom = ref_cast(Domaine_VDF,domaine_dis);
93 const Domaine_Cl_VDF& le_dom_Cl = ref_cast(Domaine_Cl_VDF,domaine_Cl_dis);
94 D = 0;
95 // return D;
96 // const DoubleVect& volumes = le_dom.volumes();
97 const DoubleVect& porosite_surf = domaine_Cl_dis.equation().milieu().porosite_face();
98 const DoubleVect& volume_entrelaces = le_dom.volumes_entrelaces();
99 // int nb_elem = le_dom.nb_elem();
100 int nb_elem_tot = le_dom.nb_elem_tot();
101 //const Domaine& domaine=le_dom.domaine();
102
103 //int nb_faces_elem = domaine.nb_faces_elem();
104 //IntTrav numfa(nb_faces_elem);
105 double coef;
106 // const IntTab& elem_faces = le_dom.elem_faces();
107 const IntTab& face_voisins = le_dom.face_voisins();
108 int nb_faces = le_dom.nb_faces();
109
110 double gradk;
111 int num_face,poly1,poly2,ori, ndeb, nfin;
112
113 // Calcul de Gradient de racine de K.
114 if (mp_min_vect(K_eps_Bas_Re)<0)
115 {
116 Cerr << "Il y'a des valeurs negatives dans les valeurs de K" << finl;
117 Cerr << "dans Modele_Jones_Launder_VDF::Calcul_D" << finl;
118 Cerr << "On arrete le calcul." << finl;
119 Process::exit();
120 }
121 // Boucle sur les bords pour traiter les conditions aux limites
122 for (int n_bord=0; n_bord<le_dom.nb_front_Cl(); n_bord++)
123
124 {
125 const Cond_lim& la_cl = le_dom_Cl.les_conditions_limites(n_bord);
126
127 if ( sub_type(Dirichlet,la_cl.valeur())||
128 sub_type(Dirichlet_homogene,la_cl.valeur()) ||
129 sub_type(Dirichlet_paroi_defilante,la_cl.valeur()) ||
130 sub_type(Echange_externe_impose,la_cl.valeur())
131 )
132
133
134 {
135 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
136 ndeb = le_bord.num_premiere_face();
137 nfin = ndeb + le_bord.nb_faces();
138
139 for (num_face=ndeb; num_face<nfin; num_face++)
140 {
141 gradk = 0;
142 poly1 = face_voisins(num_face,0);
143 if (poly1 != -1)
144 {
145 // coef = 0.5;
146 coef = volume_entrelaces(num_face)*porosite_surf(num_face)*0.5;
147 gradk = ( - sqrt(K_eps_Bas_Re(poly1,0)))/le_dom.dist_norm_bord(num_face);
148 if (!is_visco_const)
149 visco=tab_visco[poly1];
150 D[poly1] += 2*visco*(gradk*gradk)*coef;
151 }
152 else
153 {
154 poly2 = face_voisins(num_face,1);
155 // coef = 0.5;
156 coef = volume_entrelaces(num_face)*porosite_surf(num_face)*0.5;
157 gradk = ((sqrt(K_eps_Bas_Re(poly2,0)) ))/le_dom.dist_norm_bord(num_face);
158 //
159 if (!is_visco_const)
160 visco=tab_visco[poly2];
161 D[poly2] += 2*visco*(gradk*gradk)*coef;
162 }
163 }
164 }
165 else if (sub_type(Periodique,la_cl.valeur()))
166 {
167 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
168 ndeb = le_bord.num_premiere_face();
169 nfin = ndeb + le_bord.nb_faces();
170
171 for (num_face=ndeb; num_face<nfin; num_face++)
172 {
173 gradk = 0;
174 poly1 = face_voisins(num_face,0);
175 poly2 = face_voisins(num_face,1);
176 ori = le_dom.orientation(num_face);
177 // coef = 0.5;
178
179 coef = volume_entrelaces(num_face)*porosite_surf(num_face);
180
181
182 gradk = (sqrt(K_eps_Bas_Re(poly2,0))-sqrt(K_eps_Bas_Re(poly1,0)))/le_dom.dist_elem_period(poly1,poly2,ori);
183 if (!is_visco_const)
184 visco=tab_visco[poly1];
185 D[poly1] += 2*visco*(gradk*gradk)*coef;
186 if (!is_visco_const)
187 visco=tab_visco[poly2];
188 D[poly2] += 2*visco*(gradk*gradk)*coef;
189 }
190
191 }
192 else if (sub_type(Symetrie,la_cl.valeur()))
193 ;
194 else if ( (sub_type(Neumann,la_cl.valeur()))
195 ||
196 (sub_type(Neumann_homogene,la_cl.valeur()))
197 )
198 {
199 // do nothing
200 ;
201 }
202 else
203 {
204 Cerr<<la_cl->que_suis_je()<< "not implemented in calculer_D"<<finl;
205 Process::exit();
206 }
207 }
208
209 // Traitement des faces internes
210 for (num_face=le_dom.premiere_face_int(); num_face<nb_faces; num_face++)
211 {
212 poly1 = face_voisins(num_face,0);
213 poly2 = face_voisins(num_face,1);
214 ori = le_dom.orientation(num_face);
215 // coef = 0.5;
216
217 coef = volume_entrelaces(num_face)*porosite_surf(num_face);
218
219 gradk = (sqrt(K_eps_Bas_Re(poly2,0))-sqrt(K_eps_Bas_Re(poly1,0)))/(le_dom.xp(poly2,ori)- le_dom.xp(poly1,ori));
220 // Cerr<<" ici "<< num_face<< " "<<gradk*gradk/K_eps_Bas_Re(poly2,0)<<" K "<<K_eps_Bas_Re(poly2,0)/K_eps_Bas_Re(0,0)<<finl;
221 if (num_face==-396)
222 Cerr << "K_eps_Bas_Re(poly2,0)=" << K_eps_Bas_Re(poly2,0)/K_eps_Bas_Re(0,0) << " K_eps_Bas_Re(poly1,0)=" << K_eps_Bas_Re(poly1,0)/K_eps_Bas_Re(0,0) << " test "<<sqrt( K_eps_Bas_Re(poly2,0))/sqrt( K_eps_Bas_Re(0,0))<<" "<<sqrt( K_eps_Bas_Re(poly1,0))/sqrt( K_eps_Bas_Re(0,0))<< " "<<(sqrt(K_eps_Bas_Re(poly2,0))-sqrt(K_eps_Bas_Re(poly1,0)))/sqrt( K_eps_Bas_Re(0,0))<<" "<< K_eps_Bas_Re(0,0)<<finl;
223 if (!is_visco_const)
224 visco=tab_visco[poly1];
225 D[poly1] += 2*visco*(gradk*gradk)*coef;
226 if (!is_visco_const)
227 visco=tab_visco[poly2];
228 D[poly2] += 2*visco*(gradk*gradk)*coef;
229 }
230 // GF on a nb_face_elem contributions par elem ....
231 // mais 2 contributions dans chaque direction
232
233 // provisoire
234 D/=2;
235 const DoubleVect& volumes= le_dom.volumes();
236 for (int i=0; i<nb_elem_tot; i++)
237 D(i)/=volumes(i);
238 //Cerr<<D.mp_min_vect()<<" DDDDDDDDDDDDDD "<<D.mp_max_vect()<<finl;
239 //D=0;
240 return D;
241 // D abord sans le 1/3 2/3
242}
243
244// void Modele_Jones_Launder_VDF::associer_domaines(const Domaine_dis_base& domaine_dis,
245// const Domaine_Cl_dis_base& domaine_Cl_dis)
246// {
247// le_dom_VDF = ref_cast(Domaine_VDF, domaine_dis);
248// le_dom_Cl_VDF = ref_cast(Domaine_Cl_VDF, domaine_Cl_dis);
249// }
250
251
252DoubleTab& Modele_Jones_Launder_VDF::Calcul_E(DoubleTab& E,const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis, const DoubleTab& vit,const DoubleTab& K_eps_Bas_Re,const Champ_Don_base& ch_visco, const DoubleTab& visco_turb ) const
253{
254 double visco=-1;
255 const DoubleTab& tab_visco=ch_visco.valeurs();
256 int is_visco_const=sub_type(Champ_Uniforme,ch_visco);
257 if (is_visco_const)
258 visco=tab_visco(0,0);
259 const Domaine_VDF& le_dom = ref_cast(Domaine_VDF,domaine_dis);
260 const Domaine_Cl_VDF& le_dom_Cl = ref_cast(Domaine_Cl_VDF,domaine_Cl_dis);
261 E = 0;
262 //return E;
263 // provisoire
264 if (1)
265 {
266 const IntTab& face_voisins = le_dom.face_voisins();
267 const IntTab& elem_faces = le_dom.elem_faces();
268
269 const Champ_Face_VDF& vitesse = ref_cast(Champ_Face_VDF,eq_hydraulique->inconnue());
270 int nb_elem_tot=le_dom.nb_elem_tot();
271 assert (vitesse.valeurs().line_size() == 1);
272 DoubleTab gij(0,dimension,dimension, vitesse.valeurs().line_size());
273 le_dom.domaine().creer_tableau_elements(gij);
274 // Rque methode non const Pourquoi ?
275
276 ref_cast_non_const(Champ_Face_VDF,vitesse).calcul_duidxj(vitesse.valeurs(),gij,ref_cast(Domaine_Cl_VDF,eq_hydraulique->domaine_Cl_dis()));
277 gij.echange_espace_virtuel();
278
279 DoubleTab der_seconde(dimension,dimension,dimension);
280 // der_seconde (i,j,k) d2 ui dxj dxk
281 for (int elem=0; elem<nb_elem_tot; elem++)
282 {
283 // calcul de la derivee seconde
284 ArrOfInt num(dimension*2);
285 for (int dir=0; dir<dimension; dir++)
286 for (int sens=0; sens<2; sens++)
287 {
288 int f=elem_faces(elem,dir+dimension*sens);
289 int num0=face_voisins(f,sens);
290 if (num0 == -1)
291 num0=elem;
292 num[dir+dimension*sens]=num0;
293 }
294 for (int i=0; i<dimension; i++)
295 for (int j=0; j<dimension; j++)
296 for (int k=0; k<dimension; k++)
297 {
298 der_seconde(i,j,k)=0.5*(gij(num[k+dimension],i,j,0)-gij(num[k],i,j,0))/le_dom.dim_elem(elem,k);
299 }
300
301 if (!is_visco_const)
302 visco=tab_visco[elem];
303 if (dimension==2)
304 {
305 double val2=der_seconde(0,1,1)*der_seconde(0,1,1); //d2u/dy2 ^2
306 double val3=der_seconde(1,0,0)*der_seconde(1,0,0); //d2v/dx2 ^2
307
308 E[elem]= 2*visco*visco_turb(elem)*(val2+val3);
309 }
310 else
311 {
312 double val2 = der_seconde(1,0,0)+der_seconde(2,0,0); //d2v_dx2(fac)+d2w_dx2(fac);
313 val2 *= val2;
314 double val3 = der_seconde(0,1,1)+der_seconde(2,1,1); //d2u_dy2(fac)+d2w_dy2(fac);
315 val3 *= val3;
316 double val4 = der_seconde(0,2,2)+der_seconde(1,2,2); //d2u_dz2(fac)+d2v_dz2(fac);
317 val4 *= val4;
318 E[elem] += 2*visco*visco_turb(elem)*(val2+val3+val4);
319 }
320 }
321 // E/=2.;
322 //Cerr<<E.mp_min_vect()<<" E "<<E.mp_max_vect()<<finl;
323 return E;
324
325 }
326
327
328 int nb_faces_tot = le_dom.nb_faces_tot();
329 const IntTab& face_voisins = le_dom.face_voisins();
330 int n_bord,poly1, poly2,fac=-1;
331
332 DoubleTab derivee_premiere(0, Objet_U::dimension, Objet_U::dimension);
333 le_dom.creer_tableau_faces(derivee_premiere);
334 DoubleTab derivee_seconde(derivee_premiere);
335
336 calcul_derivees_premieres_croisees(derivee_premiere,domaine_dis,domaine_Cl_dis,vit );
337 calcul_derivees_secondes_croisees(derivee_seconde,domaine_dis,domaine_Cl_dis,derivee_premiere);
338
339 // traitement par faces pour avoir la contribution par element de D
340 double val2,val3,val4;
341 for (n_bord=0; n_bord<le_dom.nb_front_Cl(); n_bord++)
342 {
343 const Cond_lim& la_cl = le_dom_Cl.les_conditions_limites(n_bord);
344 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
345 int ndeb = le_bord.num_premiere_face();
346 int nfin = ndeb + le_bord.nb_faces();
347 for (fac=ndeb; fac<nfin; fac++)
348 {
349 poly1 = face_voisins(fac,0);
350 if (poly1 != -1)
351 {
352 if (Objet_U::dimension == 2 )
353 {
354 val2 = derivee_seconde(fac,0,1)*derivee_seconde(fac,0,1);
355 val3 = derivee_seconde(fac,1,0)*derivee_seconde(fac,1,0);
356 if (!is_visco_const)
357 visco=tab_visco[poly1];
358 E[poly1]+= 2*visco*visco_turb(poly1)*(val2+val3);
359 }
360 else
361 {
362 val2 = derivee_seconde(fac,1,0)+derivee_seconde(fac,2,0); //d2v_dx2(fac)+d2w_dx2(fac);
363 val2 *= val2;
364 val3 = derivee_seconde(fac,0,1)+derivee_seconde(fac,2,1); //d2u_dy2(fac)+d2w_dy2(fac);
365 val3 *= val3;
366 val4 = derivee_seconde(fac,0,2)+derivee_seconde(fac,1,2); //d2u_dz2(fac)+d2v_dz2(fac);
367 val4 *= val4;
368 if (!is_visco_const)
369 visco=tab_visco[poly1];
370 E[poly1] += 2*visco*visco_turb(poly1)*(val2+val3+val4);
371 }
372 }
373 else
374 {
375 poly2 = face_voisins(fac,1);
376 if (Objet_U::dimension == 2 )
377 {
378 val2 = derivee_seconde(fac,0,1)*derivee_seconde(fac,0,1);
379 val3 = derivee_seconde(fac,1,0)*derivee_seconde(fac,1,0);
380 if (!is_visco_const)
381 visco=tab_visco[poly2];
382 E[poly2]+= 2*visco*visco_turb(poly2)*(val2+val3);
383 }
384 else
385 {
386 val2 = derivee_seconde(fac,1,0)+derivee_seconde(fac,2,0); //d2v_dx2(fac)+d2w_dx2(fac);
387 val2 *= val2;
388 val3 = derivee_seconde(fac,0,1)+derivee_seconde(fac,2,1); //d2u_dy2(fac)+d2w_dy2(fac);
389 val3 *= val3;
390 val4 = derivee_seconde(fac,0,2)+derivee_seconde(fac,1,2); //d2u_dz2(fac)+d2v_dz2(fac);
391 val4 *= val4;
392 if (!is_visco_const)
393 visco=tab_visco[poly2];
394 E[poly2] += 2*visco*visco_turb(poly2)*(val2+val3+val4);
395 }
396 }
397 }
398 }
399
400 // Traitement des faces internes
401 for (; fac<nb_faces_tot; fac++)
402 {
403 poly1=face_voisins(fac,0);
404 poly2=face_voisins(fac,1);
405 if (Objet_U::dimension == 2 )
406 {
407 val2 = derivee_seconde(fac,0,1)*derivee_seconde(fac,0,1);
408 val3 = derivee_seconde(fac,1,0)*derivee_seconde(fac,1,0);
409
410 if (poly1>=0)
411 {
412 if (!is_visco_const)
413 visco=tab_visco[poly1];
414 E[poly1]+= 2*visco*visco_turb(poly1)*(val2+val3);
415 }
416 if (poly2>=0)
417 {
418 if (!is_visco_const)
419 visco=tab_visco[poly2];
420 E[poly2]+= 2*visco*visco_turb(poly2)*(val2+val3);
421 }
422 }
423 else
424 {
425 val2 = derivee_seconde(fac,1,0)+derivee_seconde(fac,2,0); //d2v_dx2(fac)+d2w_dx2(fac);
426 val2 *= val2;
427 val3 = derivee_seconde(fac,0,1)+derivee_seconde(fac,2,1); //d2u_dy2(fac)+d2w_dy2(fac);
428 val3 *= val3;
429 val4 = derivee_seconde(fac,0,2)+derivee_seconde(fac,1,2); //d2u_dz2(fac)+d2v_dz2(fac);
430 val4 *= val4;
431
432 if (poly1>=0)
433 {
434 if (!is_visco_const)
435 visco=tab_visco[poly1];
436 E[poly1] += 2*visco*visco_turb(poly1)*(val2+val3+val4);
437 }
438 if (poly2>=0)
439 {
440 if (!is_visco_const)
441 visco=tab_visco[poly1];
442 E[poly2] += 2*visco*visco_turb(poly2)*(val2+val3+val4);
443 }
444 }
445 }
446 return E;
447}
448
449DoubleTab& Modele_Jones_Launder_VDF::calcul_derivees_premieres_croisees(DoubleTab& derivee_premiere, const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis, const DoubleTab& vit ) const
450{
451 const Domaine_VDF& le_dom = ref_cast(Domaine_VDF,domaine_dis);
452 const Domaine_Cl_VDF& le_dom_Cl = ref_cast(Domaine_Cl_VDF,domaine_Cl_dis);
453 // Cerr<<le_dom_Cl.equation().le_nom()<<finl;exit();
454 const Champ_Face_VDF& vitesse = ref_cast(Champ_Face_VDF,eq_hydraulique->inconnue());
455 const Domaine_Cl_VDF& zcl_hydro = ref_cast(Domaine_Cl_VDF,eq_hydraulique->domaine_Cl_dis());
456 // int nb_faces = le_dom.nb_faces();
457 const IntTab& Qdm = le_dom.Qdm();
458 const IntVect& orientation = le_dom.orientation();
459 IntTrav num(4);
460 int i,num_arete;
461 double pond1,pond2;
462
463 int signe,ori0,ori1;
464 // DEBUT CALCUL DES DERIVEES PREMIERES
465
466 if (Objet_U::dimension == 2)
467 {
468 DoubleVect coef(2);
469 int ndeb,nfin;
470 // Boucle sur les aretes internes pour le calcul
471 // des moyennes des derivees croisees
472
473 ndeb = le_dom.premiere_arete_interne();
474 nfin = le_dom.nb_aretes_internes() + ndeb;
475
476 for (num_arete =ndeb; num_arete<nfin; num_arete++)
477 {
478 for (i=0; i<4; i++)
479 num[i] = Qdm(num_arete,i);
480
481 ori0 = orientation(num[0]);
482 ori1 = orientation(num[2]);
483 coef[0] = 0.5*(vit(num[1])-vit(num[0]))/le_dom.dist_face(num[0],num[1],ori1);
484 coef[1] = 0.5*(vit(num[3])-vit(num[2]))/le_dom.dist_face(num[2],num[3],ori0);
485
486 derivee_premiere(num[0],ori0,ori1) += coef[0];
487 derivee_premiere(num[1],ori0,ori1) += coef[0];
488 derivee_premiere(num[2],ori1,ori0) += coef[1];
489 derivee_premiere(num[3],ori1,ori0) += coef[1];
490 }
491
492 // Boucle sur les aretes_mixte
493 ndeb = le_dom.premiere_arete_mixte();
494 nfin = ndeb + le_dom.nb_aretes_mixtes();
495
496 for (num_arete=ndeb; num_arete<nfin; num_arete++)
497 {
498 for (i=0; i<4; i++)
499 num[i] = Qdm(num_arete,i);
500
501 ori0 = orientation(num[0]);
502 ori1 = orientation(num[2]);
503 coef[0] = 0.5*(vit(num[1])-vit(num[0]))/le_dom.dist_face(num[0],num[1],ori1);
504 coef[1] = 0.5*(vit(num[3])-vit(num[2]))/le_dom.dist_face(num[2],num[3],ori0);
505
506 derivee_premiere(num[0],ori0,ori1) += coef[0];
507 derivee_premiere(num[1],ori0,ori1) += coef[0];
508 derivee_premiere(num[2],ori1,ori0) += coef[1];
509 derivee_premiere(num[3],ori1,ori0) += coef[1];
510 }
511
512 //*******************************
513 //Prise en compte des CL
514 //*******************************
515
516 //*******************************
517 //On parcourt les aretes bords
518 //*******************************
519
520 ndeb = le_dom.premiere_arete_bord();
521 nfin = ndeb + le_dom.nb_aretes_bord();
522 int n_type;
523
524 for (num_arete=ndeb; num_arete<nfin; num_arete++)
525 {
526 n_type=le_dom_Cl.type_arete_bord(num_arete-ndeb);
527
528 for (i=0; i<4; i++)
529 num[i] = Qdm(num_arete,i);
530
531 if (n_type == TypeAreteBordVDF::PERIO_PERIO)
532 {
533 ori0 = orientation(num[0]);
534 ori1 = orientation(num[2]);
535 coef[0] = 0.5*(vit(num[1])-vit(num[0]))/le_dom.dist_face(num[0],num[1],ori1);
536 coef[1] = 0.5*(vit(num[3])-vit(num[2]))/le_dom.dist_face(num[2],num[3],ori0);
537 derivee_premiere(num[0],ori0,ori1) += coef[0];
538 derivee_premiere(num[1],ori0,ori1) += coef[0];
539 derivee_premiere(num[2],ori1,ori0) += coef[1];
540 derivee_premiere(num[3],ori1,ori0) += coef[1];
541 }
542 else
543 {
544 ori0 = orientation(num[0]);
545 ori1 = orientation(num[2]);
546 signe = num[3];
547
548 pond1 = le_dom.face_normales(num[0],ori0);
549 pond2 = le_dom.face_normales(num[1],ori0);
550 double tps=vitesse.temps();
551 double vit_imp = pond2*Champ_Face_get_val_imp_face_bord_sym(vitesse.valeurs(),tps,num[0],ori1,zcl_hydro)+
552 Champ_Face_get_val_imp_face_bord_sym(vitesse.valeurs(),tps,num[1],ori1,zcl_hydro)*pond1; // val tangentielle
553 vit_imp /= pond1+pond2;
554 // double vit_imp = 0.5*(vitesse.val_imp_face_bord(num[0],ori1)+
555 // vitesse.val_imp_face_bord(num[1],ori1)); // val tangentielle
556 coef[0] = 0.5*(vit(num[1])-vit(num[0]))/le_dom.dist_face(num[0],num[1],ori1);
557 ////coef[1] = 0.5*(vit_imp-vit(num[2]))/le_dom.dist_norm_bord(num[2])*signe;
558 coef[1] = 0.5*(vit_imp-vit(num[2]))/le_dom.dist_norm_bord(num[0])*signe;
559 derivee_premiere(num[0],ori0,ori1) += coef[0];
560 derivee_premiere(num[1],ori0,ori1) += coef[0];
561 derivee_premiere(num[2],ori1,ori0) += coef[1];
562 }
563 }
564
565
566 //*******************************
567 //On parcourt les aretes coins
568 //*******************************
569
570 ndeb = le_dom.premiere_arete_coin();
571 nfin = ndeb + le_dom.nb_aretes_coin();
572
573 for (num_arete=ndeb; num_arete<nfin; num_arete++)
574 {
575 for (i=0; i<4; i++)
576 num[i] = Qdm(num_arete,i);
577
578 n_type=le_dom_Cl.type_arete_coin(num_arete-ndeb);
579 //***************************************
580 // Traitement des aretes coin perio-perio
581 //***************************************
582
583 if (n_type == TypeAreteCoinVDF::PERIO_PERIO) // arete de type periodicite-periodicite
584 {
585 ori0 = orientation(num[0]);
586 ori1 = orientation(num[2]);
587 coef[0] = 0.5*(vit(num[1])-vit(num[0]))/le_dom.dist_face(num[0],num[1],ori1);
588 coef[1] = 0.5*(vit(num[3])-vit(num[2]))/le_dom.dist_face(num[2],num[3],ori0);
589 derivee_premiere(num[0],ori0,ori1) += coef[0];
590 derivee_premiere(num[1],ori0,ori1) += coef[0];
591 derivee_premiere(num[2],ori1,ori0) += coef[1];
592 derivee_premiere(num[3],ori1,ori0) += coef[1];
593 }
594
595 //***************************************
596 // Traitement des aretes coin perio-paroi
597 //***************************************
598 else if (n_type == TypeAreteCoinVDF::PERIO_PAROI) // arete de type periodicite-paroi
599 {
600 ori0 = orientation(num[0]);
601 ori1 = orientation(num[2]);
602 signe = num[3];
603
604 pond1 = le_dom.face_normales(num[0],ori0);
605 pond2 = le_dom.face_normales(num[1],ori0);
606 double tps=vitesse.temps();
607 double vit_imp = pond2*Champ_Face_get_val_imp_face_bord_sym(vitesse.valeurs(),tps,num[0],ori1,zcl_hydro)+
608 Champ_Face_get_val_imp_face_bord_sym(vitesse.valeurs(),tps,num[1],ori1,zcl_hydro)*pond1;
609 // val tangentielle
610 vit_imp /= pond1+pond2;
611 // double vit_imp = 0.5*(vitesse.val_imp_face_bord(num[0],ori1)+
612 // vitesse.val_imp_face_bord(num[1],ori1)); // val tangentielle
613 coef[0] = 0.5*(vit(num[1])-vit(num[0]))/le_dom.dist_face(num[0],num[1],ori1);
614 coef[1] = 0.5*(vit_imp-vit(num[2]))/le_dom.dist_norm_bord(num[0])*signe;
615 derivee_premiere(num[0],ori0,ori1) += coef[0];
616 derivee_premiere(num[1],ori0,ori1) += coef[0];
617 derivee_premiere(num[2],ori1,ori0) += coef[1];
618 }
619 else
620 {
621 //Cerr << "Oh!!! Un coin!!!" << finl;
622 }
623
624 }
625
626 // Cerr << "Il y a : " << compt_coin << " coins!!" << finl;
627 }
628
629 else if (Objet_U::dimension == 3)
630 {
631 DoubleVect coef(3);
632
633 // Boucle sur les aretes internes pour le calcul
634 // des moyennes des derivees croisees
635 int ndeb = le_dom.premiere_arete_interne();
636 int nfin = le_dom.nb_aretes_internes() + ndeb;
637
638 for (num_arete =ndeb; num_arete<nfin; num_arete++)
639 {
640 for (i=0; i<4; i++)
641 num[i] = Qdm(num_arete,i);
642
643 ori0 = orientation(num[0]);
644 ori1 = orientation(num[2]);
645 coef[0] = 0.5*(vit(num[1])-vit(num[0]))/le_dom.dist_face(num[0],num[1],ori1);
646 coef[1] = 0.5*(vit(num[3])-vit(num[2]))/le_dom.dist_face(num[2],num[3],ori0);
647
648 derivee_premiere(num[0],ori0,ori1) += coef[0];
649 derivee_premiere(num[1],ori0,ori1) += coef[0];
650 derivee_premiere(num[2],ori1,ori0) += coef[1];
651 derivee_premiere(num[3],ori1,ori0) += coef[1];
652 }
653 // Boucle sur les aretes_mixte
654 ndeb = le_dom.premiere_arete_mixte();
655 nfin = ndeb + le_dom.nb_aretes_mixtes();
656
657 for (num_arete=ndeb; num_arete<nfin; num_arete++)
658 {
659 for (i=0; i<4; i++)
660 num[i] = Qdm(num_arete,i);
661
662 ori0 = orientation(num[0]);
663 ori1 = orientation(num[2]);
664 coef[0] = 0.5*(vit(num[1])-vit(num[0]))/le_dom.dist_face(num[0],num[1],ori1);
665 coef[1] = 0.5*(vit(num[3])-vit(num[2]))/le_dom.dist_face(num[2],num[3],ori0);
666
667 derivee_premiere(num[0],ori0,ori1) += coef[0];
668 derivee_premiere(num[1],ori0,ori1) += coef[0];
669 derivee_premiere(num[2],ori1,ori0) += coef[1];
670 derivee_premiere(num[3],ori1,ori0) += coef[1];
671 }
672
673 //*******************************
674 //Prise en compte des CL
675 //*******************************
676
677 //*******************************
678 //On parcourt les aretes bords
679 //*******************************
680
681 ndeb = le_dom.premiere_arete_bord();
682 nfin = ndeb + le_dom.nb_aretes_bord();
683 int n_type;
684
685 for (num_arete=ndeb; num_arete<nfin; num_arete++)
686 {
687 n_type=le_dom_Cl.type_arete_bord(num_arete-ndeb);
688 for (i=0; i<4; i++)
689 num[i] = Qdm(num_arete,i);
690
691 if (n_type == TypeAreteBordVDF::PERIO_PERIO)
692 {
693 ori0 = orientation(num[0]);
694 ori1 = orientation(num[2]);
695 if (num[0]==num[1] || num[2]==num[3])
696 {
697 Cerr << "2 bords avec une condition limite de periodicite ne sont separees que d'une maille !" << finl;
698 Cerr << "Cela n'est pas valide pour le Modele Jones Launder et son calcul de derivees croisees..." << finl;
699 Process::exit();
700 }
701 coef[0] = 0.5*(vit(num[1])-vit(num[0]))/le_dom.dist_face(num[0],num[1],ori1);
702 coef[1] = 0.5*(vit(num[3])-vit(num[2]))/le_dom.dist_face(num[2],num[3],ori0);
703
704 derivee_premiere(num[0],ori0,ori1) += coef[0];
705 derivee_premiere(num[1],ori0,ori1) += coef[0];
706 derivee_premiere(num[2],ori1,ori0) += coef[1];
707 derivee_premiere(num[3],ori1,ori0) += coef[1];
708 }
709 else
710 {
711 ori0 = orientation(num[0]);
712 ori1 = orientation(num[2]);
713 signe = num[3];
714
715 // double vit_imp = 0.5*(vitesse.val_imp_face_bord(num[0],ori1)+
716 // vitesse.val_imp_face_bord(num[1],ori1)); // val tangentielle
717 pond1 = le_dom.face_normales(num[0],ori0);
718 pond2 = le_dom.face_normales(num[1],ori0);
719 double tps=vitesse.temps();
720 double vit_imp = pond2*Champ_Face_get_val_imp_face_bord_sym(vitesse.valeurs(),tps,num[0],ori1,zcl_hydro)+
721 Champ_Face_get_val_imp_face_bord_sym(vitesse.valeurs(),tps,num[1],ori1,zcl_hydro)*pond1; // val tangentielle
722 // val tangentielle
723 vit_imp /= pond1+pond2;
724
725 coef[0] = 0.5*(vit(num[1])-vit(num[0]))/le_dom.dist_face(num[0],num[1],ori1);
726 coef[1] = 0.5*(vit_imp-vit(num[2]))/le_dom.dist_norm_bord(num[0])*signe;
727 derivee_premiere(num[0],ori0,ori1) += coef[0];
728 derivee_premiere(num[1],ori0,ori1) += coef[0];
729 derivee_premiere(num[2],ori1,ori0) += coef[1];
730 }
731 }
732
733 //*******************************
734 //On parcourt les aretes coins
735 //*******************************
736
737 ndeb = le_dom.premiere_arete_coin();
738 nfin = ndeb + le_dom.nb_aretes_coin();
739
740 for (num_arete=ndeb; num_arete<nfin; num_arete++)
741 {
742 for (i=0; i<4; i++)
743 num[i] = Qdm(num_arete,i);
744
745 n_type=le_dom_Cl.type_arete_coin(num_arete-ndeb);
746
747 //***************************************
748 // Traitement des aretes coin perio-perio
749 //***************************************
750
751 if (n_type == TypeAreteCoinVDF::PERIO_PERIO) // arete de type periodicite-periodicite
752 {
753 ori0 = orientation(num[0]);
754 ori1 = orientation(num[2]);
755 coef[0] = 0.5*(vit(num[1])-vit(num[0]))/le_dom.dist_face(num[0],num[1],ori1);
756 coef[1] = 0.5*(vit(num[3])-vit(num[2]))/le_dom.dist_face(num[2],num[3],ori0);
757
758 derivee_premiere(num[0],ori0,ori1) += coef[0];
759 derivee_premiere(num[1],ori0,ori1) += coef[0];
760 derivee_premiere(num[2],ori1,ori0) += coef[1];
761 derivee_premiere(num[3],ori1,ori0) += coef[1];
762 }
763
764 //***************************************
765 // Traitement des aretes coin perio-paroi
766 //***************************************
767 else if (n_type == TypeAreteCoinVDF::PERIO_PAROI) // arete de type periodicite-paroi
768 {
769 ori0 = orientation(num[0]);
770 ori1 = orientation(num[2]);
771 signe = num[3];
772
773 pond1 = le_dom.face_normales(num[0],ori0);
774 pond2 = le_dom.face_normales(num[1],ori0);
775 double tps=vitesse.temps();
776 double vit_imp = pond2*Champ_Face_get_val_imp_face_bord_sym(vitesse.valeurs(),tps,num[0],ori1,zcl_hydro)+
777 Champ_Face_get_val_imp_face_bord_sym(vitesse.valeurs(),tps,num[1],ori1,zcl_hydro)*pond1;
778
779 // val tangentielle
780 vit_imp /= pond1+pond2;
781 // double vit_imp = 0.5*(vitesse.val_imp_face_bord(num[0],ori1)+
782 // vitesse.val_imp_face_bord(num[1],ori1)); // val tangentielle
783
784 coef[0] = 0.5*(vit(num[1])-vit(num[0]))/le_dom.dist_face(num[0],num[1],ori1);
785 coef[1] = 0.5*(vit_imp-vit(num[2]))/le_dom.dist_norm_bord(num[0])*signe;
786 derivee_premiere(num[0],ori0,ori1) += coef[0];
787 derivee_premiere(num[1],ori0,ori1) += coef[0];
788 derivee_premiere(num[2],ori1,ori0) += coef[1];
789 }
790 else
791 {
792 //Cerr << "Oh!!! Un coin!!!" << finl;
793 }
794 }
795 // Cerr << "Il y a : " << compt_coin << " coins!!" << finl;
796 }
797 derivee_premiere.echange_espace_virtuel();
798 return derivee_premiere;
799}
800
801DoubleTab& Modele_Jones_Launder_VDF::calcul_derivees_secondes_croisees(DoubleTab& derivee_seconde, const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis, const DoubleTab& derivee_premiere ) const
802{
803 const Domaine_VDF& le_dom = ref_cast(Domaine_VDF,domaine_dis);
804 const Domaine_Cl_VDF& le_dom_Cl = ref_cast(Domaine_Cl_VDF,domaine_Cl_dis);
805 // const Champ_Face_VDF& vitesse = ref_cast(Champ_Face_VDF,eq_hydraulique->inconnue());
806
807 // int nb_faces = le_dom.nb_faces();
808 const IntTab& Qdm = le_dom.Qdm();
809 const IntVect& orientation = le_dom.orientation();
810 // const int nb_cond_lim = le_dom_Cl.nb_cond_lim();
811 IntTrav num(4);
812 int i,num_arete,ori0,ori1;
813
814 int ndeb=-10000,nfin=100000;//,signe;
815 // DEBUT CALCUL DES DERIVEES SECONDES
816
817 if (Objet_U::dimension == 2)
818 {
819 DoubleVect coef(2);
820
821 ndeb = le_dom.premiere_arete_interne();
822 nfin = le_dom.nb_aretes_internes() + ndeb;
823 // Boucle sur les aretes internes pour le calcul
824 // des moyennes des derivees croisees
825
826 for (num_arete =ndeb; num_arete<nfin; num_arete++)
827 {
828 for (i=0; i<4; i++)
829 num[i] = Qdm(num_arete,i);
830
831 ori0 = orientation(num[0]);
832 ori1 = orientation(num[2]);
833 coef[0] = 0.5*(derivee_premiere(num[1],ori0,ori1)-derivee_premiere(num[0],ori0,ori1))/le_dom.dist_face(num[0],num[1],ori1);
834 coef[1] = 0.5*(derivee_premiere(num[3],ori1,ori0)-derivee_premiere(num[2],ori1,ori0))/le_dom.dist_face(num[2],num[3],ori0);
835
836 derivee_seconde(num[0],ori0,ori1) += coef[0];
837 derivee_seconde(num[1],ori0,ori1) += coef[0];
838 derivee_seconde(num[2],ori1,ori0) += coef[1];
839 derivee_seconde(num[3],ori1,ori0) += coef[1];
840 }
841
842 // Boucle sur les aretes_mixte
843 ndeb = le_dom.premiere_arete_mixte();
844 nfin = ndeb + le_dom.nb_aretes_mixtes();
845
846 for (num_arete=ndeb; num_arete<nfin; num_arete++)
847 {
848 for (i=0; i<4; i++)
849 num[i] = Qdm(num_arete,i);
850
851 ori0 = orientation(num[0]);
852 ori1 = orientation(num[2]);
853 coef[0] = 0.5*(derivee_premiere(num[1],ori0,ori1)-derivee_premiere(num[0],ori0,ori1))/le_dom.dist_face(num[0],num[1],ori1);
854 coef[1] = 0.5*(derivee_premiere(num[3],ori1,ori0)-derivee_premiere(num[2],ori1,ori0))/le_dom.dist_face(num[2],num[3],ori0);
855
856 derivee_seconde(num[0],ori0,ori1) += coef[0];
857 derivee_seconde(num[1],ori0,ori1) += coef[0];
858 derivee_seconde(num[2],ori1,ori0) += coef[1];
859 derivee_seconde(num[3],ori1,ori0) += coef[1];
860 }
861
862 //*******************************
863 //Prise en compte des CL
864 //*******************************
865
866 //*******************************
867 //On parcourt les aretes bords
868 //*******************************
869 ndeb = le_dom.premiere_arete_bord();
870 nfin = ndeb + le_dom.nb_aretes_bord();
871 int n_type;
872
873 for (num_arete=ndeb; num_arete<nfin; num_arete++)
874 {
875 n_type=le_dom_Cl.type_arete_bord(num_arete-ndeb);
876
877 for (i=0; i<4; i++)
878 num[i] = Qdm(num_arete,i);
879
880 if (n_type == TypeAreteBordVDF::PERIO_PERIO)
881 {
882 ori0 = orientation(num[0]);
883 ori1 = orientation(num[2]);
884 coef[0] = 0.5*(derivee_premiere(num[1],ori0,ori1)-derivee_premiere(num[0],ori0,ori1))/le_dom.dist_face(num[0],num[1],ori1);
885 coef[1] = 0.5*(derivee_premiere(num[3],ori1,ori0)-derivee_premiere(num[2],ori1,ori0))/le_dom.dist_face(num[2],num[3],ori0);
886 derivee_seconde(num[0],ori0,ori1) += coef[0];
887 derivee_seconde(num[1],ori0,ori1) += coef[0];
888 derivee_seconde(num[2],ori1,ori0) += coef[1];
889 derivee_seconde(num[3],ori1,ori0) += coef[1];
890 }
891 else
892 {
893 ori0 = orientation(num[0]);
894 ori1 = orientation(num[2]);
895 // signe = num[3];
896
897 // ///???????????????
898 // double val_deriv_prem = 0.5*(vit.val_imp_face_bord(num[0],ori1)+
899 // vit.val_imp_face_bord(num[1],ori1)); // val tangentielle
900 // ///???????????????
901 coef[0] = 0.5*(derivee_premiere(num[1],ori0,ori1)-derivee_premiere(num[0],ori0,ori1))/le_dom.dist_face(num[0],num[1],ori1);
902 //coef[1] = 0.5*(val_deriv_prem-derivee_premiere(num[2],ori1,ori0))/le_dom.dist_norm_bord(num[0])*signe;
903 derivee_seconde(num[0],ori0,ori1) += coef[0];
904 derivee_seconde(num[1],ori0,ori1) += coef[0];
905 derivee_seconde(num[2],ori1,ori0) = 2.*derivee_seconde(num[2],ori1,ori0);
906 }
907 }
908
909 //*******************************
910 //On parcourt les aretes coins
911 //*******************************
912
913 ndeb = le_dom.premiere_arete_coin();
914 nfin = ndeb + le_dom.nb_aretes_coin();
915
916 for (num_arete=ndeb; num_arete<nfin; num_arete++)
917 {
918 for (i=0; i<4; i++)
919 num[i] = Qdm(num_arete,i);
920
921 n_type=le_dom_Cl.type_arete_coin(num_arete-ndeb);
922 //***************************************
923 // Traitement des aretes coin perio-perio
924 //***************************************
925
926 if (n_type == TypeAreteCoinVDF::PERIO_PERIO) // arete de type periodicite-periodicite
927 {
928 ori0 = orientation(num[0]);
929 ori1 = orientation(num[2]);
930 coef[0] = 0.5*(derivee_premiere(num[1],ori0,ori1)-derivee_premiere(num[0],ori0,ori1))/le_dom.dist_face(num[0],num[1],ori1);
931 coef[1] = 0.5*(derivee_premiere(num[3],ori1,ori0)-derivee_premiere(num[2],ori1,ori0))/le_dom.dist_face(num[2],num[3],ori0);
932 derivee_seconde(num[0],ori0,ori1) += coef[0];
933 derivee_seconde(num[1],ori0,ori1) += coef[0];
934 derivee_seconde(num[2],ori1,ori0) += coef[1];
935 derivee_seconde(num[3],ori1,ori0) += coef[1];
936 }
937
938 //***************************************
939 // Traitement des aretes coin perio-paroi
940 //***************************************
941 else if (n_type == TypeAreteCoinVDF::PERIO_PAROI) // arete de type periodicite-paroi
942 {
943 ori0 = orientation(num[0]);
944 ori1 = orientation(num[2]);
945 //signe = num[3];
946
947 // ///???????????????
948 // double val_deriv_prem = 0.5*(vit.val_imp_face_bord(num[0],ori1)+
949 // vit.val_imp_face_bord(num[1],ori1)); // val tangentielle
950 // ///???????????????
951 coef[0] = 0.5*(derivee_premiere(num[1],ori0,ori1)-derivee_premiere(num[0],ori0,ori1))/le_dom.dist_face(num[0],num[1],ori1);
952 //coef[1] = 0.5*(val_deriv_prem-derivee_premiere(num[2],ori1,ori0))/le_dom.dist_norm_bord(num[0])*signe;
953 derivee_seconde(num[0],ori0,ori1) += coef[0];
954 derivee_seconde(num[1],ori0,ori1) += coef[0];
955 // derivee_seconde(num[2],ori1,ori0) += coef[1];
956 derivee_seconde(num[2],ori1,ori0) = 2.*derivee_seconde(num[2],ori1,ori0);
957 }
958 else
959 {
960 //Cerr << "Oh!!! Un coin!!!" << finl;
961 }
962
963 }
964 // Cerr << "Il y a : " << compt_coin << " coins!!" << finl;
965 }
966
967 else if (Objet_U::dimension == 3)
968 {
969 // Cerr << " tableau des derivee_premiere " << derivee_premiere << finl;
970 DoubleVect coef(3);
971
972 // Boucle sur les aretes internes pour le calcul
973 // des moyennes des derivees croisees
974 ndeb = le_dom.premiere_arete_interne();
975 nfin = le_dom.nb_aretes_internes() + ndeb;
976
977 for (num_arete =ndeb; num_arete<nfin; num_arete++)
978 {
979 for (i=0; i<4; i++)
980 num[i] = Qdm(num_arete,i);
981
982 ori0 = orientation(num[0]);
983 ori1 = orientation(num[2]);
984 coef[0] = 0.5*(derivee_premiere(num[1],ori0,ori1)-derivee_premiere(num[0],ori0,ori1))/le_dom.dist_face(num[0],num[1],ori1);
985 coef[1] = 0.5*(derivee_premiere(num[3],ori1,ori0)-derivee_premiere(num[2],ori1,ori0))/le_dom.dist_face(num[2],num[3],ori0);
986
987 derivee_seconde(num[0],ori0,ori1) += coef[0];
988 derivee_seconde(num[1],ori0,ori1) += coef[0];
989 derivee_seconde(num[2],ori1,ori0) += coef[1];
990 derivee_seconde(num[3],ori1,ori0) += coef[1];
991 }
992
993 // Boucle sur les aretes_mixte
994 ndeb = le_dom.premiere_arete_mixte();
995 nfin = ndeb + le_dom.nb_aretes_mixtes();
996
997 for (num_arete=ndeb; num_arete<nfin; num_arete++)
998 {
999 for (i=0; i<4; i++)
1000 num[i] = Qdm(num_arete,i);
1001
1002 ori0 = orientation(num[0]);
1003 ori1 = orientation(num[2]);
1004 coef[0] = 0.5*(derivee_premiere(num[1],ori0,ori1)-derivee_premiere(num[0],ori0,ori1))/le_dom.dist_face(num[0],num[1],ori1);
1005 coef[1] = 0.5*(derivee_premiere(num[3],ori1,ori0)-derivee_premiere(num[2],ori1,ori0))/le_dom.dist_face(num[2],num[3],ori0);
1006
1007 derivee_seconde(num[0],ori0,ori1) += coef[0];
1008 derivee_seconde(num[1],ori0,ori1) += coef[0];
1009 derivee_seconde(num[2],ori1,ori0) += coef[1];
1010 derivee_seconde(num[3],ori1,ori0) += coef[1];
1011 }
1012
1013 //*******************************
1014 //Prise en compte des CL
1015 //*******************************
1016
1017 //*******************************
1018 //On parcourt les aretes bords
1019 //*******************************
1020
1021 ndeb = le_dom.premiere_arete_bord();
1022 nfin = ndeb + le_dom.nb_aretes_bord();
1023 int n_type;
1024
1025 for (num_arete=ndeb; num_arete<nfin; num_arete++)
1026 {
1027 n_type=le_dom_Cl.type_arete_bord(num_arete-ndeb);
1028 for (i=0; i<4; i++)
1029 num[i] = Qdm(num_arete,i);
1030
1031 if (n_type == TypeAreteBordVDF::PERIO_PERIO)
1032 {
1033 ori0 = orientation(num[0]);
1034 ori1 = orientation(num[2]);
1035 coef[0] = 0.5*(derivee_premiere(num[1],ori0,ori1)-derivee_premiere(num[0],ori0,ori1))/le_dom.dist_face(num[0],num[1],ori1);
1036 coef[1] = 0.5*(derivee_premiere(num[3],ori1,ori0)-derivee_premiere(num[2],ori1,ori0))/le_dom.dist_face(num[2],num[3],ori0);
1037
1038 derivee_seconde(num[0],ori0,ori1) += coef[0];
1039 derivee_seconde(num[1],ori0,ori1) += coef[0];
1040 derivee_seconde(num[2],ori1,ori0) += coef[1];
1041 derivee_seconde(num[3],ori1,ori0) += coef[1];
1042 }
1043 else
1044 {
1045 ori0 = orientation(num[0]);
1046 ori1 = orientation(num[2]);
1047 //signe = num[3];
1048
1049 // ///???????????????
1050 // double val_deriv_prem = 0.5*(vit.val_imp_face_bord(num[0],ori1)+
1051 // vit.val_imp_face_bord(num[1],ori1)); // val tangentielle
1052 // ///???????????????
1053
1054 coef[0] = 0.5*(derivee_premiere(num[1],ori0,ori1)-derivee_premiere(num[0],ori0,ori1))/le_dom.dist_face(num[0],num[1],ori1);
1055 //coef[1] = 0.5*(val_deriv_prem-derivee_premiere(num[2],ori1,ori0))/le_dom.dist_norm_bord(num[0])*signe;
1056 derivee_seconde(num[0],ori0,ori1) += coef[0];
1057 derivee_seconde(num[1],ori0,ori1) += coef[0];
1058 // derivee_seconde(num[2],ori1,ori0) += coef[1];
1059 derivee_seconde(num[2],ori1,ori0) = 2.*derivee_seconde(num[2],ori1,ori0);
1060 }
1061 }
1062
1063 //*******************************
1064 //On parcourt les aretes coins
1065 //*******************************
1066
1067 ndeb = le_dom.premiere_arete_coin();
1068 nfin = ndeb + le_dom.nb_aretes_coin();
1069
1070
1071 for (num_arete=ndeb; num_arete<nfin; num_arete++)
1072 {
1073 for (i=0; i<4; i++)
1074 num[i] = Qdm(num_arete,i);
1075
1076 n_type=le_dom_Cl.type_arete_coin(num_arete-ndeb);
1077 //***************************************
1078 // Traitement des aretes coin perio-perio
1079 //***************************************
1080 if (n_type == TypeAreteCoinVDF::PERIO_PERIO) // arete de type periodicite-periodicite
1081 {
1082 ori0 = orientation(num[0]);
1083 ori1 = orientation(num[2]);
1084 coef[0] = 0.5*(derivee_premiere(num[1],ori0,ori1)-derivee_premiere(num[0],ori0,ori1))/le_dom.dist_face(num[0],num[1],ori1);
1085 coef[1] = 0.5*(derivee_premiere(num[3],ori1,ori0)-derivee_premiere(num[2],ori1,ori0))/le_dom.dist_face(num[2],num[3],ori0);
1086
1087 derivee_seconde(num[0],ori0,ori1) += coef[0];
1088 derivee_seconde(num[1],ori0,ori1) += coef[0];
1089 derivee_seconde(num[2],ori1,ori0) += coef[1];
1090 derivee_seconde(num[3],ori1,ori0) += coef[1];
1091 }
1092
1093 //***************************************
1094 // Traitement des aretes coin perio-paroi
1095 //***************************************
1096 else if (n_type == TypeAreteCoinVDF::PERIO_PAROI) // arete de type periodicite-paroi
1097 {
1098 ori0 = orientation(num[0]);
1099 ori1 = orientation(num[2]);
1100 //signe = num[3];
1101
1102 // ///???????????????
1103 // double val_deriv_prem = 0.5*(vit.val_imp_face_bord(num[0],ori1)+
1104 // vit.val_imp_face_bord(num[1],ori1)); // val tangentielle
1105 // ///???????????????
1106
1107 coef[0] = 0.5*(derivee_premiere(num[1],ori0,ori1)-derivee_premiere(num[0],ori0,ori1))/le_dom.dist_face(num[0],num[1],ori1);
1108 //coef[1] = 0.5*(val_deriv_prem-derivee_premiere(num[2],ori1,ori0))/le_dom.dist_norm_bord(num[0])*signe;
1109 derivee_seconde(num[0],ori0,ori1) += coef[0];
1110 derivee_seconde(num[1],ori0,ori1) += coef[0];
1111 derivee_seconde(num[2],ori1,ori0) = 2.*derivee_seconde(num[2],ori1,ori0);
1112 // derivee_seconde(num[2],ori1,ori0) += coef[1];
1113 }
1114 else
1115 {
1116 //Cerr << "Oh!!! Un coin!!!" << finl;
1117 }
1118 }
1119 // Cerr << "Il y a : " << compt_coin << " coins!!" << finl;
1120 }
1121 // Cerr << "derivee_seconde " << derivee_seconde << finl;
1122 derivee_seconde.echange_espace_virtuel();
1123 return derivee_seconde;
1124}
1125
1126DoubleTab& Modele_Jones_Launder_VDF::Calcul_F1( DoubleTab& F1, const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis, const DoubleTab& P, const DoubleTab& K_eps_Bas_Re,const Champ_base& ch_visco) const
1127{
1128 const Domaine_VDF& le_dom = ref_cast(Domaine_VDF,domaine_dis);
1129 int nb_elem = le_dom.nb_elem();
1130 for (int elem=0; elem <nb_elem; elem ++ )
1131 F1[elem] = 1.;
1132 return F1;
1133}
1134
1135DoubleTab& Modele_Jones_Launder_VDF::Calcul_F2( DoubleTab& F2, DoubleTab& Deb, const Domaine_dis_base& domaine_dis,const DoubleTab& K_eps_Bas_Re,const Champ_base& ch_visco ) const
1136{
1137 double visco=-1;
1138 const DoubleTab& tab_visco=ch_visco.valeurs();
1139 int is_visco_const=sub_type(Champ_Uniforme,ch_visco);
1140 if (is_visco_const)
1141 visco=tab_visco(0,0);
1142 const Domaine_VDF& le_dom = ref_cast(Domaine_VDF,domaine_dis);
1143 int nb_elem = le_dom.nb_elem();
1144 double Re;
1145 int elem;
1146
1147 for (elem=0; elem< nb_elem ; elem++)
1148 {
1149 if (!is_visco_const)
1150 visco=tab_visco[elem];
1151 if (visco>0)
1152 {
1153 Re = K_eps_Bas_Re(elem,0)*K_eps_Bas_Re(elem,0)/(K_eps_Bas_Re(elem,1)+DMINFLOAT)/visco;
1154 F2[elem] = 1. - (0.3*exp(-1.*carre(Re)));
1155 }
1156 else
1157 F2[elem] = 1.;
1158 }
1159 // F2=1;
1160// Cerr<<F2.mp_min_vect()<<" F2 "<<F2.mp_max_vect()<<finl;
1161 return F2;
1162}
1163/*
1164 DoubleTab& Modele_Jones_Launder_VDF::Calcul_F2( DoubleTab& F2, DoubleTab& D, const Domaine_dis_base& domaine_dis,const DoubleTab& K_eps_Bas_Re, const DoubleTab& tab_visco ) const
1165 {
1166 Process::exit();
1167 const Domaine_VDF& le_dom = ref_cast(Domaine_VDF,domaine_dis);
1168 int nb_elem = le_dom.nb_elem();
1169 DoubleTab Re(nb_elem);
1170 int elem;
1171
1172 for (elem=0; elem< nb_elem ; elem++) {
1173 if (K_eps_Bas_Re(elem,1)>0) {
1174 Re(elem) = (K_eps_Bas_Re(elem,0)*K_eps_Bas_Re(elem,0))/(tab_visco(elem)*K_eps_Bas_Re(elem,1));
1175 F2[elem] = 1. - (0.3*exp(-1*carre(Re(elem))));
1176 } else {
1177 F2[elem] = 1.;
1178 }
1179 }
1180 return F2;
1181 }
1182*/
1183
1184DoubleTab& Modele_Jones_Launder_VDF::Calcul_Fmu( DoubleTab& Fmu,const Domaine_dis_base& domaine_dis,const Domaine_Cl_dis_base& domaine_Cl_dis,const DoubleTab& K_eps_Bas_Re,const Champ_Don_base& ch_visco ) const
1185{
1186 double visco=-1;
1187 const DoubleTab& tab_visco=ch_visco.valeurs();
1188 int is_visco_const=sub_type(Champ_Uniforme,ch_visco);
1189 if (is_visco_const)
1190 visco=tab_visco(0,0);
1191 const Domaine_VDF& le_dom = ref_cast(Domaine_VDF,domaine_dis);
1192 Fmu = 0;
1193 int nb_elem = le_dom.nb_elem();
1194 double Re;
1195 int elem;
1196 // Cerr << " Calc Fmu " << finl;
1197 for (elem=0; elem< nb_elem ; elem++)
1198 {
1199 if (!is_visco_const)
1200 visco=tab_visco[elem];
1201 if (visco>0)
1202 {
1203 Re = K_eps_Bas_Re(elem,0)*K_eps_Bas_Re(elem,0)/(K_eps_Bas_Re(elem,1)+DMINFLOAT)/visco;
1204 Fmu[elem] = exp(-2.5/(1.+Re/50.));
1205 }
1206 else
1207 Fmu[elem] = 1;
1208 // provisoire
1209 // Fmu[elem] = exp(-3.4/((1.+Re/50.)*(1.+Re/50.)));
1210 }
1211 //Cerr<<Fmu.mp_min_vect()<<" fmuuuuuuuuuuuuuuu " <<Fmu.mp_max_vect()<<finl;
1212 /*
1213 Cerr<<K_eps_Bas_Re(0,0)<<" ke "<<K_eps_Bas_Re(0,1)<<finl;
1214 Cerr<<"re "<< K_eps_Bas_Re(0,0)*K_eps_Bas_Re(0,0)/K_eps_Bas_Re(0,1)/visco;
1215 Cerr<<"visco "<<visco<<finl;
1216 */
1217 // Fmu=1;
1218 return Fmu;
1219}
1220/*
1221 DoubleTab& Modele_Jones_Launder_VDF::Calcul_Fmu( DoubleTab& Fmu,const Domaine_dis_base& domaine_dis,const DoubleTab& K_eps_Bas_Re,const DoubleTab& tab_visco ) const
1222 {
1223 Process::exit();
1224 const Domaine_VDF& le_dom = ref_cast(Domaine_VDF,domaine_dis);
1225 double Re;
1226 Fmu = 0;
1227 int nb_elem = le_dom.nb_elem();
1228 int elem;
1229 for (elem=0; elem< nb_elem ; elem++) {
1230 // Fmu[elem] = exp(-2.5/(1.+K_eps_Bas_Re(elem,0)*K_eps_Bas_Re(elem,0)/(tab_visco(elem)*K_eps_Bas_Re(elem,1))));
1231 Re = K_eps_Bas_Re(elem,0)*K_eps_Bas_Re(elem,0)/K_eps_Bas_Re(elem,1)/tab_visco(elem);
1232 Fmu[elem] = exp(-2.5/(1.+Re/50.));
1233 Fmu[elem] = exp(-3.4/((1.+Re/50.)*(1.+Re/50.)));
1234
1235 }
1236 return Fmu;
1237 }
1238*/
1239void Modele_Jones_Launder_VDF::mettre_a_jour(double temps)
1240{
1241 ;
1242}
1243
1244DoubleTab& Modele_Jones_Launder_VDF::Calcul_Fmu_BiK( DoubleTab& Fmu,const Domaine_dis_base& domaine_dis,const Domaine_Cl_dis_base& domaine_Cl_dis,const DoubleTab& K_Bas_Re,const DoubleTab& eps_Bas_Re,const Champ_Don_base& ch_visco ) const
1245{
1246 double visco=-1;
1247 const DoubleTab& tab_visco=ch_visco.valeurs();
1248 int is_visco_const=sub_type(Champ_Uniforme,ch_visco);
1249 if (is_visco_const)
1250 visco=tab_visco(0,0);
1251 const Domaine_VDF& le_dom = ref_cast(Domaine_VDF,domaine_dis);
1252 Fmu = 0;
1253 int nb_elem = le_dom.nb_elem();
1254 double Re;
1255 int elem;
1256 // Cerr << " Calc Fmu " << finl;
1257 for (elem=0; elem< nb_elem ; elem++)
1258 {
1259 if (!is_visco_const)
1260 visco=tab_visco[elem];
1261 if (visco>0)
1262 {
1263 Re = K_Bas_Re(elem)*K_Bas_Re(elem)/(eps_Bas_Re(elem)+DMINFLOAT)/visco;
1264 Fmu[elem] = exp(-2.5/(1.+Re/50.));
1265 }
1266 else
1267 Fmu[elem] = 1;
1268 // provisoire
1269 // Fmu[elem] = exp(-3.4/((1.+Re/50.)*(1.+Re/50.)));
1270 }
1271 //Cerr<<Fmu.mp_min_vect()<<" fmuuuuuuuuuuuuuuu " <<Fmu.mp_max_vect()<<finl;
1272 /*
1273 Cerr<<K_eps_Bas_Re(0,0)<<" ke "<<K_eps_Bas_Re(0,1)<<finl;
1274 Cerr<<"re "<< K_eps_Bas_Re(0,0)*K_eps_Bas_Re(0,0)/K_eps_Bas_Re(0,1)/visco;
1275 Cerr<<"visco "<<visco<<finl;
1276 */
1277 // Fmu=1;
1278 return Fmu;
1279}
1280
1281
1282DoubleTab& Modele_Jones_Launder_VDF::Calcul_F2_BiK( DoubleTab& F2, DoubleTab& Deb, const Domaine_dis_base& domaine_dis,const DoubleTab& K_Bas_Re,const DoubleTab& eps_Bas_Re,const Champ_base& ch_visco ) const
1283{
1284 double visco=-1;
1285 const DoubleTab& tab_visco=ch_visco.valeurs();
1286 int is_visco_const=sub_type(Champ_Uniforme,ch_visco);
1287 if (is_visco_const)
1288 visco=tab_visco(0,0);
1289 const Domaine_VDF& le_dom = ref_cast(Domaine_VDF,domaine_dis);
1290 int nb_elem = le_dom.nb_elem();
1291 double Re;
1292 int elem;
1293
1294 for (elem=0; elem< nb_elem ; elem++)
1295 {
1296 if (!is_visco_const)
1297 visco=tab_visco[elem];
1298 if (visco>0)
1299 {
1300 Re = K_Bas_Re(elem)*K_Bas_Re(elem)/(eps_Bas_Re(elem)+DMINFLOAT)/visco;
1301 F2[elem] = 1. - (0.3*exp(-1.*carre(Re)));
1302 }
1303 else
1304 F2[elem] = 1.;
1305 }
1306 // F2=1;
1307// Cerr<<F2.mp_min_vect()<<" F2 "<<F2.mp_max_vect()<<finl;
1308 return F2;
1309}
1310
1311
1312
1313DoubleTab& Modele_Jones_Launder_VDF::Calcul_F1_BiK( DoubleTab& F1, const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis, const DoubleTab& P, const DoubleTab& K_Bas_Re, const DoubleTab& eps_Bas_Re,const Champ_base& ch_visco) const
1314{
1315 const Domaine_VDF& le_dom = ref_cast(Domaine_VDF,domaine_dis);
1316 int nb_elem = le_dom.nb_elem();
1317 for (int elem=0; elem <nb_elem; elem ++ )
1318 F1[elem] = 1.;
1319 return F1;
1320}
1321
1322
1323DoubleTab& Modele_Jones_Launder_VDF::Calcul_E_BiK(DoubleTab& E,const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis, const DoubleTab& vit,const DoubleTab& K_Bas_Re,const DoubleTab& eps_Bas_Re,const Champ_Don_base& ch_visco, const DoubleTab& visco_turb ) const
1324{
1325 return Calcul_E( E, domaine_dis, domaine_Cl_dis, vit, K_Bas_Re, ch_visco, visco_turb );
1326}
1327
1328
1329DoubleTab& Modele_Jones_Launder_VDF::Calcul_D_BiK(DoubleTab& D,const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis,
1330 const DoubleTab& vitesse,const DoubleTab& K_Bas_Re,const DoubleTab& eps_Bas_Re, const Champ_Don_base& ch_visco ) const
1331{
1332 double visco=-1;
1333 const DoubleTab& tab_visco=ch_visco.valeurs();
1334 int is_visco_const=sub_type(Champ_Uniforme,ch_visco);
1335 if (is_visco_const)
1336 visco=tab_visco(0,0);
1337
1338 const Domaine_VDF& le_dom = ref_cast(Domaine_VDF,domaine_dis);
1339 const Domaine_Cl_VDF& le_dom_Cl = ref_cast(Domaine_Cl_VDF,domaine_Cl_dis);
1340 D = 0;
1341 // return D;
1342 // const DoubleVect& volumes = le_dom.volumes();
1343 const DoubleVect& porosite_surf = domaine_Cl_dis.equation().milieu().porosite_face();
1344 const DoubleVect& volume_entrelaces = le_dom.volumes_entrelaces();
1345 // int nb_elem = le_dom.nb_elem();
1346 int nb_elem_tot = le_dom.nb_elem_tot();
1347 //const Domaine& domaine=le_dom.domaine();
1348
1349 //int nb_faces_elem = domaine.nb_faces_elem();
1350 //IntTrav numfa(nb_faces_elem);
1351 double coef;
1352 // const IntTab& elem_faces = le_dom.elem_faces();
1353 const IntTab& face_voisins = le_dom.face_voisins();
1354 int nb_faces = le_dom.nb_faces();
1355
1356 double gradk;
1357 int num_face,poly1,poly2,ori, ndeb, nfin;
1358
1359 // Calcul de Gradient de racine de K.
1360 if (mp_min_vect(K_Bas_Re)<0)
1361 {
1362 Cerr << "Il y'a des valeurs negatives dans les valeurs de K" << finl;
1363 Cerr << "dans Modele_Jones_Launder_VDF::Calcul_D" << finl;
1364 Cerr << "On arrete le calcul." << finl;
1365 Process::exit();
1366 }
1367 // Boucle sur les bords pour traiter les conditions aux limites
1368 for (int n_bord=0; n_bord<le_dom.nb_front_Cl(); n_bord++)
1369
1370 {
1371 const Cond_lim& la_cl = le_dom_Cl.les_conditions_limites(n_bord);
1372
1373 if ( sub_type(Dirichlet,la_cl.valeur())||
1374 sub_type(Dirichlet_homogene,la_cl.valeur()) ||
1375 sub_type(Dirichlet_paroi_defilante,la_cl.valeur()) ||
1376 sub_type(Echange_externe_impose,la_cl.valeur())
1377 )
1378
1379
1380 {
1381 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1382 ndeb = le_bord.num_premiere_face();
1383 nfin = ndeb + le_bord.nb_faces();
1384
1385 for (num_face=ndeb; num_face<nfin; num_face++)
1386 {
1387 gradk = 0;
1388 poly1 = face_voisins(num_face,0);
1389 if (poly1 != -1)
1390 {
1391 // coef = 0.5;
1392 coef = volume_entrelaces(num_face)*porosite_surf(num_face)*0.5;
1393 gradk = ( - sqrt(K_Bas_Re(poly1)))/le_dom.dist_norm_bord(num_face);
1394 if (!is_visco_const)
1395 visco=tab_visco[poly1];
1396 D[poly1] += 2*visco*(gradk*gradk)*coef;
1397 }
1398 else
1399 {
1400 poly2 = face_voisins(num_face,1);
1401 // coef = 0.5;
1402 coef = volume_entrelaces(num_face)*porosite_surf(num_face)*0.5;
1403 gradk = ((sqrt(K_Bas_Re(poly2)) ))/le_dom.dist_norm_bord(num_face);
1404 //
1405 if (!is_visco_const)
1406 visco=tab_visco[poly2];
1407 D[poly2] += 2*visco*(gradk*gradk)*coef;
1408 }
1409 }
1410 }
1411 else if (sub_type(Periodique,la_cl.valeur()))
1412 {
1413 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1414 ndeb = le_bord.num_premiere_face();
1415 nfin = ndeb + le_bord.nb_faces();
1416
1417 for (num_face=ndeb; num_face<nfin; num_face++)
1418 {
1419 gradk = 0;
1420 poly1 = face_voisins(num_face,0);
1421 poly2 = face_voisins(num_face,1);
1422 ori = le_dom.orientation(num_face);
1423 // coef = 0.5;
1424
1425 coef = volume_entrelaces(num_face)*porosite_surf(num_face);
1426
1427
1428 gradk = (sqrt(K_Bas_Re(poly2))-sqrt(K_Bas_Re(poly1)))/le_dom.dist_elem_period(poly1,poly2,ori);
1429 if (!is_visco_const)
1430 visco=tab_visco[poly1];
1431 D[poly1] += 2*visco*(gradk*gradk)*coef;
1432 if (!is_visco_const)
1433 visco=tab_visco[poly2];
1434 D[poly2] += 2*visco*(gradk*gradk)*coef;
1435 }
1436
1437 }
1438 else if (sub_type(Symetrie,la_cl.valeur()))
1439 ;
1440 else if ( (sub_type(Neumann,la_cl.valeur()))
1441 ||
1442 (sub_type(Neumann_homogene,la_cl.valeur()))
1443 )
1444 {
1445 // do nothing
1446 ;
1447 }
1448 else
1449 {
1450 Cerr<<la_cl->que_suis_je()<< "not implemented in calculer_D"<<finl;
1451 Process::exit();
1452 }
1453 }
1454
1455 // Traitement des faces internes
1456 for (num_face=le_dom.premiere_face_int(); num_face<nb_faces; num_face++)
1457 {
1458 poly1 = face_voisins(num_face,0);
1459 poly2 = face_voisins(num_face,1);
1460 ori = le_dom.orientation(num_face);
1461 // coef = 0.5;
1462
1463 coef = volume_entrelaces(num_face)*porosite_surf(num_face);
1464
1465 gradk = (sqrt(K_Bas_Re(poly2))-sqrt(K_Bas_Re(poly1)))/(le_dom.xp(poly2,ori)- le_dom.xp(poly1,ori));
1466 // Cerr<<" ici "<< num_face<< " "<<gradk*gradk/K_eps_Bas_Re(poly2,0)<<" K "<<K_eps_Bas_Re(poly2,0)/K_eps_Bas_Re(0,0)<<finl;
1467 if (num_face==-396)
1468 Cerr << "K_eps_Bas_Re(poly2,0)=" << K_Bas_Re(poly2)/K_Bas_Re(0) << " K_eps_Bas_Re(poly1,0)=" << K_Bas_Re(poly1)/K_Bas_Re(0) << " test "<<sqrt( K_Bas_Re(poly2,0))/sqrt( K_Bas_Re(0))<<" "<<sqrt( K_Bas_Re(poly1))/sqrt( K_Bas_Re(0))<< " "<<(sqrt(K_Bas_Re(poly2))-sqrt(K_Bas_Re(poly1)))/sqrt( K_Bas_Re(0))<<" "<< K_Bas_Re(0)<<finl;
1469 if (!is_visco_const)
1470 visco=tab_visco[poly1];
1471 D[poly1] += 2*visco*(gradk*gradk)*coef;
1472 if (!is_visco_const)
1473 visco=tab_visco[poly2];
1474 D[poly2] += 2*visco*(gradk*gradk)*coef;
1475 }
1476 // GF on a nb_face_elem contributions par elem ....
1477 // mais 2 contributions dans chaque direction
1478
1479 // provisoire
1480 D/=2;
1481 const DoubleVect& volumes= le_dom.volumes();
1482 for (int i=0; i<nb_elem_tot; i++)
1483 D(i)/=volumes(i);
1484 //Cerr<<D.mp_min_vect()<<" DDDDDDDDDDDDDD "<<D.mp_max_vect()<<finl;
1485 //D=0;
1486 return D;
1487 // D abord sans le 1/3 2/3
1488}
1489
Champ_Don_base
classe Champ_Don_base classe de base des Champs donnes (non calcules)
Definition Champ_Don_base.h:32
Champ_Don_base::valeurs
DoubleTab & valeurs() override
Surcharge Champ_base::valeurs() Renvoie le tableau des valeurs.
Definition Champ_Don_base.h:49
Champ_Face_VDF
class Champ_Face_VDF Cette classe sert a representer un champ vectoriel dont on ne calcule
Definition Champ_Face_VDF.h:42
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
Definition Champ_Inc_base.h:77
Champ_Proto::valeurs
virtual DoubleTab & valeurs()=0
Champ_Uniforme
classe Champ_Uniforme Represente un champ constant dans l'espace et dans le temps.
Definition Champ_Uniforme.h:26
Champ_base
classe Champ_base Cette classe est la base de la hierarchie des champs.
Definition Champ_base.h:43
Champ_base::temps
double temps() const
Renvoie le temps du champ.
Definition Champ_base.cpp:1004
Cond_lim
classe Cond_lim Classe generique servant a representer n'importe quelle classe
Definition Cond_lim.h:31
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_paroi_defilante
classe Dirichlet_paroi_defilante Impose la vitesse de paroi dnas une equation de type Navier_Stokes.
Definition Dirichlet_paroi_defilante.h:26
Dirichlet
classe Dirichlet Cette classe est la classe de base de la hierarchie des conditions aux limites de ty...
Definition Dirichlet.h:31
Domaine_32_64::creer_tableau_elements
virtual void creer_tableau_elements(Array_base &, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT) const
creation d'un tableau parallele de valeurs aux elements.
Definition Domaine.cpp:851
Domaine_Cl_VDF
class Domaine_Cl_VDF
Definition Domaine_Cl_VDF.h:63
Domaine_Cl_VDF::type_arete_bord
int type_arete_bord(int num_arete) const
Definition Domaine_Cl_VDF.h:76
Domaine_Cl_VDF::type_arete_coin
const int & type_arete_coin(int num_arete) const
Definition Domaine_Cl_VDF.h:78
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::nb_aretes_coin
int nb_aretes_coin() const
Definition Domaine_VDF.h:514
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_VDF::premiere_arete_bord
int premiere_arete_bord() const
Definition Domaine_VDF.h:538
Domaine_VDF::premiere_arete_coin
int premiere_arete_coin() const
Definition Domaine_VDF.h:522
Domaine_VDF::nb_aretes_internes
int nb_aretes_internes() const
Definition Domaine_VDF.h:562
Domaine_VDF::premiere_arete_interne
int premiere_arete_interne() const
Definition Domaine_VDF.h:578
Domaine_VDF::nb_aretes_bord
int nb_aretes_bord() const
Definition Domaine_VDF.h:530
Domaine_VDF::nb_aretes_mixtes
int nb_aretes_mixtes() const
Definition Domaine_VDF.h:546
Domaine_VDF::dist_face
double dist_face(int, int, int k) const
Definition Domaine_VDF.h:308
Domaine_VDF::dist_elem_period
double dist_elem_period(int, int, int) const
Definition Domaine_VDF.h:618
Domaine_VDF::premiere_arete_mixte
int premiere_arete_mixte() const
Definition Domaine_VDF.h:554
Domaine_VDF::Qdm
int Qdm(int num_arete, int) const
Definition Domaine_VDF.h:231
Domaine_VDF::face_normales
double face_normales(int, int) const override
Definition Domaine_VDF.h:187
Domaine_VDF::dist_norm_bord
double dist_norm_bord(int num_face) const override
Definition Domaine_VDF.h:382
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::nb_faces_tot
int nb_faces_tot() const
renvoie le nombre total de faces.
Definition Domaine_VF.h:481
Domaine_VF::creer_tableau_faces
void creer_tableau_faces(Array_base &, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT) const
Definition Domaine_VF.cpp:881
Domaine_VF::volumes
double volumes(int i) const
Definition Domaine_VF.h:113
Domaine_VF::elem_faces
int elem_faces(int i, int j) const
renvoie le numero de le ieme face de la maille num_elem la facon dont ces faces sont numerotees est
Definition Domaine_VF.h:543
Domaine_VF::xp
double xp(int num_elem, int k) const
Definition Domaine_VF.h:77
Domaine_VF::premiere_face_int
int premiere_face_int() const
une face est interne ssi elle separe deux elements.
Definition Domaine_VF.h:463
Domaine_VF::face_voisins
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
Definition Domaine_VF.h:418
Domaine_dis_base
classe Domaine_dis_base Cette classe est la base de la hierarchie des domaines discretisees.
Definition Domaine_dis_base.h:39
Domaine_dis_base::nb_elem_tot
int nb_elem_tot() const
Definition Domaine_dis_base.h:50
Domaine_dis_base::nb_elem
int nb_elem() const
Definition Domaine_dis_base.h:49
Domaine_dis_base::nb_front_Cl
int nb_front_Cl() const
Definition Domaine_dis_base.h:53
Domaine_dis_base::domaine
const Domaine & domaine() const
Definition Domaine_dis_base.h:46
Echange_externe_impose
Classe Echange_externe_impose: Cette classe represente le cas particulier de la classe.
Definition Echange_externe_impose.h:46
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
Modele_Fonc_Bas_Reynolds_Base
Definition Modele_Fonc_Bas_Reynolds_Base.h:38
Modele_Jones_Launder_VDF
Definition Modele_Jones_Launder_VDF.h:31
Modele_Jones_Launder_VDF::Calcul_F2_BiK
DoubleTab & Calcul_F2_BiK(DoubleTab &, DoubleTab &, const Domaine_dis_base &, const DoubleTab &, const DoubleTab &, const Champ_base &) const override
Definition Modele_Jones_Launder_VDF.cpp:1282
Modele_Jones_Launder_VDF::Calcul_E
DoubleTab & Calcul_E(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &, const DoubleTab &, const Champ_Don_base &, const DoubleTab &) const override
Definition Modele_Jones_Launder_VDF.cpp:252
Modele_Jones_Launder_VDF::Calcul_D_BiK
DoubleTab & Calcul_D_BiK(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const Champ_Don_base &) const override
Definition Modele_Jones_Launder_VDF.cpp:1329
Modele_Jones_Launder_VDF::Calcul_F2
DoubleTab & Calcul_F2(DoubleTab &, DoubleTab &, const Domaine_dis_base &, const DoubleTab &, const Champ_base &) const override
Definition Modele_Jones_Launder_VDF.cpp:1135
Modele_Jones_Launder_VDF::Calcul_Fmu
DoubleTab & Calcul_Fmu(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &, const Champ_Don_base &) const override
Definition Modele_Jones_Launder_VDF.cpp:1184
Modele_Jones_Launder_VDF::Calcul_F1
DoubleTab & Calcul_F1(DoubleTab &F1, const Domaine_dis_base &domaine_dis, const Domaine_Cl_dis_base &domaine_Cl_dis, const DoubleTab &P, const DoubleTab &K_eps_Bas_Re, const Champ_base &ch_visco) const override
Definition Modele_Jones_Launder_VDF.cpp:1126
Modele_Jones_Launder_VDF::Calcul_F1_BiK
DoubleTab & Calcul_F1_BiK(DoubleTab &F1, const Domaine_dis_base &domaine_dis, const Domaine_Cl_dis_base &domaine_Cl_dis, const DoubleTab &P, const DoubleTab &K_Bas_Re, const DoubleTab &eps_Bas_Re, const Champ_base &ch_visco) const override
Definition Modele_Jones_Launder_VDF.cpp:1313
Modele_Jones_Launder_VDF::Calcul_E_BiK
DoubleTab & Calcul_E_BiK(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const Champ_Don_base &, const DoubleTab &) const override
Definition Modele_Jones_Launder_VDF.cpp:1323
Modele_Jones_Launder_VDF::calcul_derivees_premieres_croisees
DoubleTab & calcul_derivees_premieres_croisees(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &) const
Definition Modele_Jones_Launder_VDF.cpp:449
Modele_Jones_Launder_VDF::Calcul_D
DoubleTab & Calcul_D(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &, const DoubleTab &, const Champ_Don_base &) const override
Definition Modele_Jones_Launder_VDF.cpp:83
Modele_Jones_Launder_VDF::associer
void associer(const Domaine_dis_base &, const Domaine_Cl_dis_base &) override
Definition Modele_Jones_Launder_VDF.cpp:74
Modele_Jones_Launder_VDF::calcul_derivees_secondes_croisees
DoubleTab & calcul_derivees_secondes_croisees(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &) const
Definition Modele_Jones_Launder_VDF.cpp:801
Modele_Jones_Launder_VDF::Calcul_Fmu_BiK
DoubleTab & Calcul_Fmu_BiK(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &, const DoubleTab &, const Champ_Don_base &) const override
Definition Modele_Jones_Launder_VDF.cpp:1244
Modele_Jones_Launder_VDF::mettre_a_jour
void mettre_a_jour(double) override
Definition Modele_Jones_Launder_VDF.cpp:1239
Modele_Jones_Launder_VDF::lire
Entree & lire(const Motcle &, Entree &)
Definition Modele_Jones_Launder_VDF.cpp:66
MorEqn::equation
const Equation_base & equation() const
Renvoie la reference sur l'equation pointe par MorEqn::mon_equation.
Definition MorEqn.h:62
Motcle
Une chaine de caractere (Nom) en majuscules.
Definition Motcle.h:26
Neumann_homogene
Classe Neumann_homogene Cette classe est la classe de base de la hierarchie des conditions aux limite...
Definition Neumann_homogene.h:31
Neumann
Classe Neumann Cette classe est la classe de base de la hierarchie des conditions aux limites de type...
Definition Neumann.h:31
Objet_U::Entree
friend class Entree
Definition Objet_U.h:76
Objet_U::dimension
static int dimension
Definition Objet_U.h:99
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
Periodique
classe Periodique Cette classe represente une condition aux limites periodique.
Definition Periodique.h:31
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
Symetrie
classe Symetrie Sur les faces de symetrie on a les proprietes suivantes:
Definition Symetrie.h:37
TRUSTVect::line_size
int line_size() const
Definition TRUSTVect.tpp:67
TRUSTVect::echange_espace_virtuel
virtual void echange_espace_virtuel(IsExchangeBlocking exchange_type=IsExchangeBlocking::DefaultBlocking, const std::string kernel_name="noname")
Definition TRUSTVect.tpp:282