TRUST 1.9.8
HPC thermohydraulic platform
Loading...
Searching...
No Matches
Op_Dift_VEF_Face_Q1.cpp
1/****************************************************************************
2* Copyright (c) 2026, CEA
3* All rights reserved.
4*
5* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
6* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
7* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
8* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
9*
10* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
11* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
12* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
13*
14*****************************************************************************/
15
16#include <Dirichlet_paroi_defilante.h>
17#include <Echange_externe_impose.h>
18#include <Dirichlet_paroi_fixe.h>
19#include <Op_Dift_VEF_Face_Q1.h>
20#include <Champ_Uniforme.h>
21#include <Domaine_Cl_VEF.h>
22#include <Neumann_paroi.h>
23#include <Milieu_base.h>
24#include <Champ_Q1NC.h>
25#include <Periodique.h>
26#include <TRUSTTrav.h>
27
28#include <Debog.h>
29
30Implemente_instanciable(Op_Dift_VEF_Face_Q1, "Op_Dift_VEF_Q1NC", Op_Dift_VEF_base);
31
32Sortie& Op_Dift_VEF_Face_Q1::printOn(Sortie& s) const { return s << que_suis_je(); }
33
34Entree& Op_Dift_VEF_Face_Q1::readOn(Entree& s) { return s; }
35
36static double viscA_Q1(const Domaine_VEF& le_dom, int num_face, int num2, int dimension, int num_elem, double diffu)
37{
38 double DSiSj = le_dom.face_normales(num_face, 0) * le_dom.face_normales(num2, 1) - le_dom.face_normales(num_face, 1) * le_dom.face_normales(num2, 0);
39
40 if (dimension == 3)
41 DSiSj += ((le_dom.face_normales(num_face, 1) * le_dom.face_normales(num2, 2) - le_dom.face_normales(num_face, 2) * le_dom.face_normales(num2, 1))
42 + (le_dom.face_normales(num_face, 2) * le_dom.face_normales(num2, 0) - le_dom.face_normales(num_face, 0) * le_dom.face_normales(num2, 2)));
43
44 DSiSj = diffu / le_dom.volumes(num_elem) * std::fabs(DSiSj);
45 return DSiSj;
46}
47
48DoubleTab& Op_Dift_VEF_Face_Q1::ajouter(const DoubleTab& inconnue, DoubleTab& resu) const
49{
50 const Domaine_Cl_VEF& domaine_Cl_VEF = domaine_cl_vef();
51 const Domaine_VEF& domaine_VEF = domaine_vef();
52 const DoubleVect& porosite_face = equation().milieu().porosite_face();
53 const IntTab& elem_faces = domaine_VEF.elem_faces(), &face_voisins = domaine_VEF.face_voisins();
54 const IntVect& rang_elem_non_std = domaine_VEF.rang_elem_non_std();
55 const int n1 = domaine_VEF.nb_faces(), nb_faces_elem = domaine_VEF.domaine().nb_faces_elem(), nb_comp = resu.line_size();
56 const double mu = diffusivite(0);
57 const DoubleTab& mu_turb = diffusivite_turbulente().valeurs(), &face_normale = domaine_VEF.face_normales();
58 DoubleVect n(dimension);
59 DoubleTrav Tgrad(dimension, dimension);
60
61 // On traite les faces bord
62 if (nb_comp == 1)
63 for (int n_bord = 0; n_bord < domaine_VEF.nb_front_Cl(); n_bord++)
64 {
65 Debog::verifier("Op_Dift_VEF_Face_Q1::ajouter apres nb_compo= 1 nbords, resu", resu);
66
67 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
68 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
69 int nb_faces_bord = le_bord.nb_faces();
70 int num1 = 0;
71 int num2 = nb_faces_bord;
72
73 if (sub_type(Periodique, la_cl.valeur()))
74 {
75 const Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
76 int fac_asso;
77 for (int ind_face = num1; ind_face < num2; ind_face++)
78 {
79 const int num_face = le_bord.num_face(ind_face), elem1 = face_voisins(num_face, 0);
80 fac_asso = le_bord.num_face(la_cl_perio.face_associee(ind_face));
81 double d_mu = mu + mu_turb(elem1);
82 for (int i = 0; i < nb_faces_elem; i++)
83 {
84 const int j = elem_faces(elem1, i);
85 if ((j > num_face) && (j != fac_asso))
86 {
87 const double valA = viscA_Q1(domaine_VEF, num_face, j, dimension, elem1, d_mu);
88 resu(num_face, 0) += valA * inconnue(j, 0) * porosite_face(num_face);
89 resu(num_face, 0) -= valA * inconnue(num_face, 0) * porosite_face(num_face);
90 resu(j, 0) += 0.5 * valA * inconnue(num_face, 0) * porosite_face(j);
91 resu(j, 0) -= 0.5 * valA * inconnue(j, 0) * porosite_face(j);
92 }
93 }
94 const int elem2 = face_voisins(num_face, 1);
95 d_mu = mu + mu_turb(elem2);
96 for (int i = 0; i < nb_faces_elem; i++)
97 {
98 const int j = elem_faces(elem2, i);
99 if ((j > num_face) && (j != fac_asso))
100 {
101 const double valA = viscA_Q1(domaine_VEF, num_face, j, dimension, elem2, d_mu);
102 resu(num_face, 0) += valA * inconnue(j, 0) * porosite_face(num_face);
103 resu(num_face, 0) -= valA * inconnue(num_face, 0) * porosite_face(num_face);
104 resu(j, 0) += 0.5 * valA * inconnue(num_face, 0) * porosite_face(j);
105 resu(j, 0) -= 0.5 * valA * inconnue(j, 0) * porosite_face(j);
106 }
107 }
108 }
109 }
110 else
111 for (int ind_face = num1; ind_face < num2; ind_face++)
112 {
113 const int num_face = le_bord.num_face(ind_face), elem1 = face_voisins(num_face, 0);
114 const double d_mu = mu + mu_turb(elem1);
115 for (int i = 0; i < nb_faces_elem; i++)
116 {
117 const int j = elem_faces(elem1, i);
118 if (j > num_face)
119 {
120 const double valA = viscA_Q1(domaine_VEF, num_face, j, dimension, elem1, d_mu);
121 resu(num_face, 0) += valA * inconnue(j, 0) * porosite_face(num_face);
122 resu(num_face, 0) -= valA * inconnue(num_face, 0) * porosite_face(num_face);
123 resu(j, 0) += valA * inconnue(num_face, 0) * porosite_face(j);
124 resu(j, 0) -= valA * inconnue(j, 0) * porosite_face(j);
125 }
126 }
127 }
128 }
129 else // nb_comp > 1
130 for (int n_bord = 0; n_bord < domaine_VEF.nb_front_Cl(); n_bord++)
131 {
132 Debog::verifier("Op_Dift_VEF_Face_Q1::ajouter apres nb_compo> 1 nbords, resu", resu);
133
134 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
135 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
136 //const IntTab& elem_faces = domaine_VEF.elem_faces();
137 int nb_faces_bord_tot = le_bord.nb_faces_tot();
138 int nb_faces_bord = le_bord.nb_faces();
139 int num1 = 0;
140 int num2 = nb_faces_bord;
141
142 if (sub_type(Periodique, la_cl.valeur()))
143 {
144 const Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
145 int fac_asso;
146 for (int ind_face = num1; ind_face < num2; ind_face++)
147 {
148 const int num_face = le_bord.num_face(ind_face), elem1 = face_voisins(num_face, 0);
149 fac_asso = le_bord.num_face(la_cl_perio.face_associee(ind_face));
150 double d_mu = mu + mu_turb(elem1);
151 for (int i = 0; i < nb_faces_elem; i++)
152 {
153 const int j = elem_faces(elem1, i);
154 if ((j > num_face) && (j != fac_asso))
155 {
156 const double valA = viscA_Q1(domaine_VEF, num_face, j, dimension, elem1, d_mu);
157 for (int nc = 0; nc < nb_comp; nc++)
158 {
159 resu(num_face, nc) += valA * inconnue(j, nc) * porosite_face(num_face);
160 resu(num_face, nc) -= valA * inconnue(num_face, nc) * porosite_face(num_face);
161 resu(j, nc) += 0.5 * valA * inconnue(num_face, nc) * porosite_face(j);
162 resu(j, nc) -= 0.5 * valA * inconnue(j, nc) * porosite_face(j);
163 }
164 }
165 }
166 const int elem2 = face_voisins(num_face, 1);
167 d_mu = mu + mu_turb(elem2);
168 for (int i = 0; i < nb_faces_elem; i++)
169 {
170 const int j = elem_faces(elem2, i);
171 if ((j > num_face) && (j != fac_asso))
172 {
173 const double valA = viscA_Q1(domaine_VEF, num_face, j, dimension, elem2, d_mu);
174 for (int nc = 0; nc < nb_comp; nc++)
175 {
176 resu(num_face, nc) += valA * inconnue(j, nc) * porosite_face(num_face);
177 resu(num_face, nc) -= valA * inconnue(num_face, nc) * porosite_face(num_face);
178 resu(j, nc) += 0.5 * valA * inconnue(num_face, nc) * porosite_face(j);
179 resu(j, nc) -= 0.5 * valA * inconnue(j, nc) * porosite_face(j);
180 }
181 }
182 }
183 }
184 Debog::verifier("Op_Dift_VEF_Face_Q1::ajouter apres periodique , resu", resu);
185
186 }
187 else if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) || sub_type(Dirichlet_paroi_defilante, la_cl.valeur()))
188 {
189 if (le_modele_turbulence->utiliser_loi_paroi())
190 {
191 if (dimension == 2)
192 {
193 for (int ind_face = num1; ind_face < num2; ind_face++)
194 {
195 const int num_face = le_bord.num_face(ind_face);
196 int elem = face_voisins(num_face, 0);
197 // Calcul du gradient gradU=tau*Ny*N
198 // Au bord, n toujours oriente vers l'exterieur
199 n[0] = face_normale(num_face, 0);
200 n[1] = face_normale(num_face, 1);
201 n /= norme_array(n);
202 double signe = 0;
203 Tgrad(0, 0) = tau_tan_(num_face, 0) * n[1] * n[0];
204 Tgrad(0, 1) = tau_tan_(num_face, 0) * n[1] * n[1];
205 Tgrad(1, 0) = -tau_tan_(num_face, 0) * n[0] * n[0];
206 Tgrad(1, 1) = -tau_tan_(num_face, 0) * n[0] * n[1];
207 // On determine le sens du frottement parietal :
208 signe = 0;
209 for (int i = 0; i < nb_faces_elem; i++)
210 {
211 const int j = elem_faces(elem, i);
212 // Calcul du signe de (Utan-Uparoi)*tan :
213 signe += (inconnue(j, 1) - inconnue(num_face, 1)) * face_normale(num_face, 0) - (inconnue(j, 0) - inconnue(num_face, 0)) * face_normale(num_face, 1);
214 }
215 if (signe < 0)
216 signe = -1.;
217 else
218 signe = 1.;
219 // On calcule la contribution sur chaque face de l'element :
220 for (int i = 0; i < nb_faces_elem; i++)
221 {
222 const int j = elem_faces(elem, i);
223 double ori = 1.;
224 if (domaine_VEF.face_voisins(j, 0) != elem)
225 ori = -1;
226 if (j != num_face)
227 for (int nc = 0; nc < nb_comp; nc++)
228 resu(j, nc) -= ori * signe * (Tgrad(nc, 0) * face_normale(j, 0) + Tgrad(nc, 1) * face_normale(j, 1)) * porosite_face(j);
229 }
230 }
231 Debog::verifier("Op_Dift_VEF_Face_Q1::ajouter apres dirichlet dim2, resu", resu);
232 }
233 else if (dimension == 3)
234 {
235 for (int ind_face = num1; ind_face < num2; ind_face++)
236 {
237 const int num_face = le_bord.num_face(ind_face);
238 int elem = face_voisins(num_face, 0);
239 // Calcul du gradient gradU
240 // Au bord, n toujours oriente vers l'exterieur
241 n[0] = face_normale(num_face, 0);
242 n[1] = face_normale(num_face, 1);
243 n[2] = face_normale(num_face, 2);
244 n /= norme_array(n);
245 // Difference avec le 2D car ici tau_tan est signe !
246 // A tester malgre tout !
247 Tgrad(0, 0) = tau_tan_(num_face, 0) * n[0];
248 Tgrad(0, 1) = tau_tan_(num_face, 0) * n[1];
249 Tgrad(0, 2) = tau_tan_(num_face, 0) * n[2];
250 Tgrad(1, 0) = tau_tan_(num_face, 1) * n[0];
251 Tgrad(1, 1) = tau_tan_(num_face, 1) * n[1];
252 Tgrad(1, 2) = tau_tan_(num_face, 1) * n[2];
253 Tgrad(2, 0) = tau_tan_(num_face, 2) * n[0];
254 Tgrad(2, 1) = tau_tan_(num_face, 2) * n[1];
255 Tgrad(2, 2) = tau_tan_(num_face, 2) * n[2];
256 // On calcule la contribution sur chaque face de l'element :
257 for (int i = 0; i < nb_faces_elem; i++)
258 {
259 const int j = elem_faces(elem, i);
260 double ori = 1.;
261 if (domaine_VEF.face_voisins(j, 0) != elem)
262 ori = -1;
263 if (j != num_face)
264 for (int nc = 0; nc < nb_comp; nc++)
265 resu(j, nc) -= ori * (Tgrad(nc, 0) * face_normale(j, 0) + Tgrad(nc, 1) * face_normale(j, 1) + Tgrad(nc, 2) * face_normale(j, 2)) * porosite_face(j);
266 }
267 }
268 Debog::verifier("Op_Dift_VEF_Face_Q1::ajouter apres dirichlet dim3, resu", resu);
269 }
270 }
271 Debog::verifier("Op_Dift_VEF_Face_Q1::ajouter apres dirichlet (je crois), resu", resu);
272 }
273 else
274 {
275 Debog::verifier("Op_Dift_VEF_Face_Q1::ajouter avant else des bords, resu", resu);
276 num1 = 0;
277 num2 = nb_faces_bord_tot;
278 for (int ind_face = num1; ind_face < num2; ind_face++)
279 {
280 const int num_face = le_bord.num_face(ind_face);
281 int elem = face_voisins(num_face, 0);
282 const double d_mu = mu + mu_turb(elem);
283 // Boucle sur les faces :
284 for (int ii = 0; ii < nb_faces_elem; ii++)
285 {
286 const int j = elem_faces(elem, ii);
287 if (((j > num_face) && (ind_face < nb_faces_bord)) || ((j != num_face) && (ind_face >= nb_faces_bord)))
288 {
289 const double valA = viscA_Q1(domaine_VEF, num_face, j, dimension, elem, d_mu);
290 for (int nc = 0; nc < nb_comp; nc++)
291 {
292 resu(num_face, nc) += valA * inconnue(j, nc) * porosite_face(num_face);
293 resu(num_face, nc) -= valA * inconnue(num_face, nc) * porosite_face(num_face);
294 resu(j, nc) += valA * inconnue(num_face, nc) * porosite_face(j);
295 resu(j, nc) -= valA * inconnue(j, nc) * porosite_face(j);
296 }
297 }
298 }
299 }
300 Debog::verifier("Op_Dift_VEF_Face_Q1::ajouter apres else des bords, resu", resu);
301 }
302 }
303
304 Debog::verifier("Op_Dift_VEF_Face_Q1::ajouter apres bords, resu", resu);
305
306 // On traite les faces internes
307
308 for (int num_face = domaine_VEF.premiere_face_int(); num_face < n1; num_face++)
309 for (int kk = 0; kk < 2; kk++)
310 {
311 const int elem0 = face_voisins(num_face, kk);
312 const double d_mu = mu + mu_turb(elem0);
313 // On elimine les elements avec CL de paroi (rang>=1)
314 int rang = rang_elem_non_std(elem0);
315 if (rang < 1)
316 {
317 for (int i = 0; i < nb_faces_elem; i++)
318 {
319 const int j = elem_faces(elem0, i);
320 if (j > num_face)
321 {
322 const double valA = viscA_Q1(domaine_VEF, num_face, j, dimension, elem0, d_mu);
323 for (int nc = 0; nc < nb_comp; nc++)
324 {
325 resu(num_face, nc) += valA * inconnue(j, nc) * porosite_face(num_face);
326 resu(num_face, nc) -= valA * inconnue(num_face, nc) * porosite_face(num_face);
327 resu(j, nc) += valA * inconnue(num_face, nc) * porosite_face(j);
328 resu(j, nc) -= valA * inconnue(j, nc) * porosite_face(j);
329 }
330 }
331 }
332 }
333 }
334
335 Debog::verifier("Op_Dift_VEF_Face_Q1::ajouter apres interne, resu", resu);
336
337 for (int n_bord = 0; n_bord < domaine_VEF.nb_front_Cl(); n_bord++)
338 {
339 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
340 if (sub_type(Neumann_paroi, la_cl.valeur()))
341 {
342 const Neumann_paroi& la_cl_paroi = ref_cast(Neumann_paroi, la_cl.valeur());
343 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
344 int ndeb = le_bord.num_premiere_face();
345 int nfin = ndeb + le_bord.nb_faces();
346 for (int face = ndeb; face < nfin; face++)
347 for (int comp = 0; comp < nb_comp; comp++)
348 resu(face, comp) += la_cl_paroi.flux_impose(face - ndeb, comp) * domaine_VEF.surface(face) * porosite_face(face);
349 }
350 else if (sub_type(Echange_externe_impose, la_cl.valeur()))
351 {
352 const Echange_externe_impose& la_cl_paroi = ref_cast(Echange_externe_impose, la_cl.valeur());
353 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
354 int ndeb = le_bord.num_premiere_face();
355 int nfin = ndeb + le_bord.nb_faces();
356 for (int face = ndeb; face < nfin; face++)
357 for (int comp = 0; comp < nb_comp; comp++)
358 resu(face, comp) += la_cl_paroi.h_imp(face - ndeb, comp) * (la_cl_paroi.T_ext(face - ndeb) - inconnue(face, comp)) * domaine_VEF.surface(face) * porosite_face(face);
359 }
360 }
361
362 Debog::verifier("Op_Dift_VEF_Face_Q1::ajouter apres fin, resu", resu);
363 modifier_flux(*this);
364 return resu;
365}
366
367void Op_Dift_VEF_Face_Q1::contribuer_a_avec(const DoubleTab& transporte, Matrice_Morse& matrice) const
368{
369 const Domaine_Cl_VEF& domaine_Cl_VEF = domaine_cl_vef();
370 const Domaine_VEF& domaine_VEF = domaine_vef();
371 const DoubleVect& porosite_face = equation().milieu().porosite_face();
372 const IntTab& elem_faces = domaine_VEF.elem_faces(), &face_voisins = domaine_VEF.face_voisins();
373 const IntVect& rang_elem_non_std = domaine_VEF.rang_elem_non_std();
374 const int nb_comp = transporte.line_size(), nb_faces_elem = domaine_VEF.domaine().nb_faces_elem(), n1 = domaine_VEF.nb_faces();
375 const DoubleTab& mu_turb = diffusivite_turbulente().valeurs();
376 const double mu = diffusivite(0);
377
378 //DoubleVect n(dimension);
379 //DoubleTrav Tgrad(dimension, dimension);
380 auto& tab1 = matrice.get_set_tab1();
381 auto& tab2 = matrice.get_set_tab2();
382 auto& coeff = matrice.get_set_coeff();
383
384 // On traite les faces bord
385 for (int n_bord = 0; n_bord < domaine_VEF.nb_front_Cl(); n_bord++)
386 {
387 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
388 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
389 const int num1 = le_bord.num_premiere_face(), num2 = num1 + le_bord.nb_faces();
390
391 if (sub_type(Periodique, la_cl.valeur()))
392 {
393 const Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
394 for (int num_face0 = num1; num_face0 < num2; num_face0++)
395 {
396 const int fac_asso = la_cl_perio.face_associee(num_face0 - num1) + num1, elem1 = face_voisins(num_face0, 0);
397 double d_mu = mu + mu_turb(elem1);
398 for (int i = 0; i < nb_faces_elem; i++)
399 {
400 const int j = elem_faces(elem1, i);
401 if ((j > num_face0) && (j != fac_asso))
402 {
403 const double valA = viscA_Q1(domaine_VEF, num_face0, j, dimension, elem1, d_mu);
404 for (int nc = 0; nc < nb_comp; nc++)
405 {
406 for (auto kk = tab1[num_face0 * nb_comp + nc] - 1; kk < tab1[num_face0 * nb_comp + nc + 1] - 1; kk++)
407 {
408 if (tab2[kk] - 1 == num_face0 * nb_comp + nc)
409 coeff(kk) += valA * porosite_face(num_face0 * nb_comp + nc);
410 if (tab2[kk] - 1 == j * nb_comp + nc)
411 coeff(kk) -= 0.5 * valA * porosite_face(j * nb_comp + nc);
412 }
413 for (auto kk = tab1[j * nb_comp + nc] - 1; kk < tab1[j * nb_comp + nc + 1] - 1; kk++)
414 {
415 if (tab2[kk] - 1 == num_face0 * nb_comp + nc)
416 coeff(kk) -= valA * porosite_face(num_face0 * nb_comp + nc);
417 if (tab2[kk] - 1 == j * nb_comp + nc)
418 coeff(kk) += 0.5 * valA * porosite_face(j * nb_comp + nc);
419 }
420 }
421 }
422 }
423 const int elem2 = face_voisins(num_face0, 1);
424 d_mu = mu + mu_turb(elem2);
425 for (int i = 0; i < nb_faces_elem; i++)
426 {
427 const int j = elem_faces(elem2, i);
428 if ((j > num_face0) && (j != fac_asso))
429 {
430 const double valA = viscA_Q1(domaine_VEF, num_face0, j, dimension, elem2, d_mu);
431 for (int nc = 0; nc < nb_comp; nc++)
432 {
433
434 for (auto kk = tab1[num_face0 * nb_comp + nc] - 1; kk < tab1[num_face0 * nb_comp + nc + 1] - 1; kk++)
435 {
436 if (tab2[kk] - 1 == num_face0 * nb_comp + nc)
437 coeff(kk) += valA * porosite_face(num_face0 * nb_comp + nc);
438 if (tab2[kk] - 1 == j * nb_comp + nc)
439 coeff(kk) -= 0.5 * valA * porosite_face(j * nb_comp + nc);
440 }
441 for (auto kk = tab1[j * nb_comp + nc] - 1; kk < tab1[j * nb_comp + nc + 1] - 1; kk++)
442 {
443 if (tab2[kk] - 1 == num_face0 * nb_comp + nc)
444 coeff(kk) -= valA * porosite_face(num_face0 * nb_comp + nc);
445 if (tab2[kk] - 1 == j * nb_comp + nc)
446 coeff(kk) += 0.5 * valA * porosite_face(j * nb_comp + nc);
447 }
448 }
449 }
450 }
451 }
452 }
453 else
454 for (int num_face = num1; num_face < num2; num_face++)
455 {
456 int elembis = face_voisins(num_face, 0);
457 const double d_mu = mu + mu_turb(elembis);
458 // Boucle sur les faces :
459 for (int ii = 0; ii < nb_faces_elem; ii++)
460 {
461 const int j = elem_faces(elembis, ii);
462 if (j > num_face)
463 {
464 const double valA = viscA_Q1(domaine_VEF, num_face, j, dimension, elembis, d_mu);
465 for (int nc = 0; nc < nb_comp; nc++)
466 {
467 for (auto kk = tab1[num_face * nb_comp + nc] - 1; kk < tab1[num_face * nb_comp + nc + 1] - 1; kk++)
468 {
469 if (tab2[kk] - 1 == num_face * nb_comp + nc)
470 coeff(kk) += valA * porosite_face(num_face * nb_comp + nc);
471 if (tab2[kk] - 1 == j * nb_comp + nc)
472 coeff(kk) -= valA * porosite_face(j * nb_comp + nc);
473 }
474 for (auto kk = tab1[j * nb_comp + nc] - 1; kk < tab1[j * nb_comp + nc + 1] - 1; kk++)
475 {
476 if (tab2[kk] - 1 == num_face * nb_comp + nc)
477 coeff(kk) -= valA * porosite_face(num_face * nb_comp + nc);
478 if (tab2[kk] - 1 == j * nb_comp + nc)
479 coeff(kk) += valA * porosite_face(j * nb_comp + nc);
480 }
481 }
482 }
483 }
484 }
485 }
486
487 // On traite les faces internes
488 for (int num_face0 = domaine_VEF.premiere_face_int(); num_face0 < n1; num_face0++)
489 for (int ll = 0; ll < 2; ll++)
490 {
491 const int elem = face_voisins(num_face0, ll);
492 const double d_mu = mu + mu_turb(elem);
493 // On elimine les elements avec CL de paroi (rang>=1)
494 int rang = rang_elem_non_std(elem);
495 if (rang < 1)
496 {
497 for (int i = 0; i < nb_faces_elem; i++)
498 {
499 const int j = elem_faces(elem, i);
500 if (j > num_face0)
501 {
502 const double valA = viscA_Q1(domaine_VEF, num_face0, j, dimension, elem, d_mu);
503 for (int nc = 0; nc < nb_comp; nc++)
504 {
505 for (auto kk = tab1[num_face0 * nb_comp + nc] - 1; kk < tab1[num_face0 * nb_comp + nc + 1] - 1; kk++)
506 {
507 if (tab2[kk] - 1 == num_face0 * nb_comp + nc)
508 coeff(kk) += valA * porosite_face(num_face0 * nb_comp + nc);
509 if (tab2[kk] - 1 == j * nb_comp + nc)
510 coeff(kk) -= valA * porosite_face(j * nb_comp + nc);
511 }
512 for (auto kk = tab1[j * nb_comp + nc] - 1; kk < tab1[j * nb_comp + nc + 1] - 1; kk++)
513 {
514 if (tab2[kk] - 1 == num_face0 * nb_comp + nc)
515 coeff(kk) -= valA * porosite_face(num_face0 * nb_comp + nc);
516 if (tab2[kk] - 1 == j * nb_comp + nc)
517 coeff(kk) += valA * porosite_face(j * nb_comp + nc);
518 }
519 }
520 }
521 }
522 }
523 }
524}
Champ_Don_base::valeurs
DoubleTab & valeurs() override
Surcharge Champ_base::valeurs() Renvoie le tableau des valeurs.
Definition Champ_Don_base.h:49
Cond_lim
classe Cond_lim Classe generique servant a representer n'importe quelle classe
Definition Cond_lim.h:31
Debog::verifier
static void verifier(const char *const msg, double)
Definition Debog.cpp:21
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_paroi_fixe
classe Dirichlet_paroi_fixe Represente une paroi immobile dans une equation de type Navier_Stokes.
Definition Dirichlet_paroi_fixe.h:26
Domaine_32_64::nb_faces_elem
int nb_faces_elem(int=0) const
Renvoie le nombre de face de type i des elements geometriques constituants le domaine.
Definition Domaine.h:484
Domaine_Cl_VEF
Definition Domaine_Cl_VEF.h:40
Domaine_Cl_dis_base::les_conditions_limites
const Cond_lim & les_conditions_limites(int) const
Renvoie la i-ieme condition aux limites.
Definition Domaine_Cl_dis_base.cpp:386
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:54
Domaine_VEF::rang_elem_non_std
IntVect & rang_elem_non_std()
Definition Domaine_VEF.h:86
Domaine_VF::nb_faces
int nb_faces() const
renvoie le nombre global de faces.
Definition Domaine_VF.h:471
Domaine_VF::face_normales
virtual double face_normales(int face, int comp) const
Definition Domaine_VF.h:47
Domaine_VF::volumes
double volumes(int i) const
Definition Domaine_VF.h:113
Domaine_VF::surface
virtual double surface(int i) const
Definition Domaine_VF.h:53
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::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
Echange_impose_base::h_imp
virtual double h_imp(int num) const
Renvoie la valeur du coefficient d'echange de chaleur impose sur la i-eme composante.
Definition Echange_impose_base.cpp:112
Echange_impose_base::T_ext
virtual double T_ext(int num) const
Renvoie la valeur de la temperature imposee sur la i-eme composante du champ de frontiere.
Definition Echange_impose_base.cpp:76
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
Front_VF::nb_faces_tot
int nb_faces_tot() const
Definition Front_VF.h:58
Front_VF::num_face
int num_face(const int) const
Definition Front_VF.h:68
Matrice_Morse
Classe Matrice_Morse Represente une matrice M (creuse), non necessairement carree.
Definition Matrice_Morse.h:50
Matrice_Morse::get_set_tab2
auto & get_set_tab2()
Definition Matrice_Morse.h:103
Matrice_Morse::get_set_coeff
auto & get_set_coeff()
Definition Matrice_Morse.h:108
Matrice_Morse::get_set_tab1
auto & get_set_tab1()
Definition Matrice_Morse.h:98
Milieu_base::porosite_face
DoubleVect & porosite_face()
Definition Milieu_base.h:62
MorEqn::equation
const Equation_base & equation() const
Renvoie la reference sur l'equation pointe par MorEqn::mon_equation.
Definition MorEqn.h:62
Neumann_paroi
Classe Neumann_paroi Cette condition limite correspond a un flux impose pour l'equation de.
Definition Neumann_paroi.h:29
Neumann::flux_impose
virtual double flux_impose(int i) const
Renvoie la valeur du flux impose sur la i-eme composante du champ representant le flux a la frontiere...
Definition Neumann.cpp:35
Objet_U::dimension
static int dimension
Definition Objet_U.h:99
Objet_U::que_suis_je
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
Definition Objet_U.cpp:104
Objet_U::readOn
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
Definition Objet_U.cpp:293
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
Op_Diff_Turbulent_base::diffusivite_turbulente
const Champ_Fonc_base & diffusivite_turbulente() const
Definition Op_Diff_Turbulent_base.h:32
Op_Diff_VEF_base::domaine_vef
const Domaine_VEF & domaine_vef() const
Definition Op_Diff_VEF_base.h:57
Op_Diff_VEF_base::domaine_cl_vef
const Domaine_Cl_VEF & domaine_cl_vef() const
Definition Op_Diff_VEF_base.h:58
Op_Dift_VEF_Face_Q1
Definition Op_Dift_VEF_Face_Q1.h:23
Op_Dift_VEF_Face_Q1::ajouter
DoubleTab & ajouter(const DoubleTab &, DoubleTab &) const override
Definition Op_Dift_VEF_Face_Q1.cpp:48
Op_Dift_VEF_Face_Q1::contribuer_a_avec
void contribuer_a_avec(const DoubleTab &, Matrice_Morse &) const override
DOES NOTHING - to override in derived classes.
Definition Op_Dift_VEF_Face_Q1.cpp:367
Op_Dift_VEF_base
Definition Op_Dift_VEF_base.h:25
Op_Dift_VEF_base::diffusivite
const Champ_base & diffusivite() const override
Definition Op_Dift_VEF_base.h:33
Op_Dift_VEF_base::tau_tan_
DoubleTab tau_tan_
Definition Op_Dift_VEF_base.h:62
Op_VEF_Face::modifier_flux
void modifier_flux(const Operateur_base &) const
Definition Op_VEF_Face.cpp:338
Periodique
classe Periodique Cette classe represente une condition aux limites periodique.
Definition Periodique.h:31
Periodique::face_associee
int face_associee(int i) const
Definition Periodique.h:35
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52
TRUSTVect::line_size
int line_size() const
Definition TRUSTVect.tpp:67