TrioCFD 1.9.9_beta
TrioCFD documentation
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 // Process boundary faces
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 // Compute the gradient gradU=tau*Ny*N
198 // At the boundary, n is always oriented outward
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 // Determine the direction of the wall friction:
208 signe = 0;
209 for (int i = 0; i < nb_faces_elem; i++)
210 {
211 const int j = elem_faces(elem, i);
212 // Compute the sign of (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 // Compute the contribution on each face of the 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 // Compute the gradient gradU
240 // At the boundary, n is always oriented outward
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 with 2D because here tau_tan is signed!
246 // To test nevertheless!
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 // Compute the contribution on each face of the 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 // Loop over 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 // Process the internal faces
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 // Skip elements with wall boundary conditions (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 // Process boundary faces
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 // Loop over 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 // Process the internal faces
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 // Skip elements with wall boundary conditions (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
Overrides Champ_base::valeurs() Returns the array of values.
Definition Champ_Don_base.h:49
Cond_lim
class Cond_lim Generic class used to represent any class
Definition Cond_lim.h:31
Debog::verifier
static void verifier(const char *const msg, double)
Definition Debog.cpp:21
Dirichlet_paroi_defilante
Dirichlet_paroi_defilante Imposes the wall velocity in an equation of type Navier_Stokes.
Definition Dirichlet_paroi_defilante.h:26
Dirichlet_paroi_fixe
Dirichlet_paroi_fixe Represents a fixed wall in a Navier-Stokes type equation.
Definition Dirichlet_paroi_fixe.h:26
Domaine_32_64::nb_faces_elem
int nb_faces_elem(int=0) const
Returns the number of faces of type i of the geometric elements that make up the domain.
Definition Domaine.h:484
Domaine_Cl_VEF
Definition Domaine_Cl_VEF.h:38
Domaine_Cl_dis_base::les_conditions_limites
const Cond_lim & les_conditions_limites(int) const
Returns the i-th boundary condition.
Definition Domaine_Cl_dis_base.cpp:386
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:53
Domaine_VEF::rang_elem_non_std
IntVect & rang_elem_non_std()
Definition Domaine_VEF.h:85
Domaine_VF::nb_faces
int nb_faces() const
Returns the total number of faces.
Definition Domaine_VF.h:471
Domaine_VF::face_normales
virtual double face_normales(int face, int comp) const
Definition Domaine_VF.h:47
Domaine_VF::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
Returns the index of the i-th face of element num_elem; the face numbering convention is.
Definition Domaine_VF.h:542
Domaine_VF::premiere_face_int
int premiere_face_int() const
A face is internal if and only if it separates two elements.
Definition Domaine_VF.h:463
Domaine_VF::face_voisins
int face_voisins(int num_face, int i) const
Returns the neighbouring element of num_face in direction i.
Definition Domaine_VF.h:418
Domaine_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: This class represents the special case of the class.
Definition Echange_externe_impose.h:46
Echange_impose_base::h_imp
virtual double h_imp(int num) const
Returns the value of the imposed heat exchange coefficient on the i-th component.
Definition Echange_impose_base.cpp:112
Echange_impose_base::T_ext
virtual double T_ext(int num) const
Returns the value of the imposed temperature on the i-th component of the boundary field.
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
Matrice_Morse class - Represents a (sparse) matrix M, not necessarily square,.
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
Returns the reference to the equation pointed to by MorEqn::mon_equation.
Definition MorEqn.h:62
Neumann_paroi
Classe Neumann_paroi This boundary condition corresponds to an imposed flux for the.
Definition Neumann_paroi.h:29
Neumann::flux_impose
virtual double flux_impose(int i) const
Returns the value of the imposed flux on the i-th component of the field representing the flux at the...
Definition Neumann.cpp:35
Objet_U::dimension
static int dimension
Definition Objet_U.h:94
Objet_U::que_suis_je
const Nom & que_suis_je() const
Returns the string identifying the class.
Definition Objet_U.cpp:104
Objet_U::readOn
virtual Entree & readOn(Entree &)
Reads an Objet_U from an input stream. Virtual method to override.
Definition Objet_U.cpp:289
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Writes the object to an output stream. Virtual method to override.
Definition Objet_U.cpp:278
Op_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
class Periodique This class represents a periodic boundary condition.
Definition Periodique.h:31
Periodique::face_associee
int face_associee(int i) const
Definition Periodique.h:35
Sortie
Base class for output streams.
Definition Sortie.h:52
TRUSTVect::line_size
int line_size() const
Definition TRUSTVect.tpp:67