TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Eval_Conv_VDF_Face.tpp
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#ifndef Eval_Conv_VDF_Face_TPP_included
17#define Eval_Conv_VDF_Face_TPP_included
18
19/* ************************************** *
20 * ********* POUR L'EXPLICITE ********** *
21 * ************************************** */
22
23template <typename DERIVED_T> template<Type_Flux_Fa7 Fa7_Type, typename Type_Double> inline std::enable_if_t< Fa7_Type == Type_Flux_Fa7::SORTIE_LIBRE, void>
24Eval_Conv_VDF_Face<DERIVED_T>::flux_fa7(const DoubleTab& inco, const DoubleTab* a_r, int face, const Neumann_sortie_libre& la_cl, int num1, Type_Double& flux) const
25{
26 const int elem1 = elem_(face, 0), elem2 = elem_(face,1);
27 for (int k = 0; k < flux.size_array(); k++)
28 {
29 double psc = dt_vitesse(face, k) * surface(face);
30 if (a_r && DERIVED_T::IS_AMONT) psc *= (*a_r)((elem1 != -1) ? elem1 : elem2, k);
31 flux[k] = -psc * inco(face, k) * porosite(face);
32 }
33}
34
35template <typename DERIVED_T> template<Type_Flux_Fa7 Fa7_Type, typename Type_Double> inline std::enable_if_t< Fa7_Type == Type_Flux_Fa7::ELEM, void>
36Eval_Conv_VDF_Face<DERIVED_T>::flux_fa7(const DoubleTab& inco, const DoubleTab* a_r, int num_elem, int fac1, int fac2, Type_Double& flux) const
37{
38 const int ncomp = flux.size_array();
39 double psc = 0.25*(dt_vitesse(fac1)+dt_vitesse(fac2))*(surface(fac1)+surface(fac2));
40 if (DERIVED_T::IS_AMONT)
41 {
42 for (int k = 0; k < ncomp; k++)
43 {
44 psc = 0.25*(dt_vitesse(fac1,k)+dt_vitesse(fac2,k))*(surface(fac1)+surface(fac2));
45 const int f = psc > 0 ? fac1 : fac2;
46
47 if (a_r)
48 {
49 const int elem = elem_(f, 0), elem2 = elem_(f, 1);
50 const int e = dt_vitesse(f,k) > 0 ? (elem > -1 ? elem : elem2) : (elem2 > -1 ? elem2 : elem);
51 psc *= (*a_r)(e, k);
52 }
53
54 flux[k] = -psc * inco(f, k) * porosite(f);
55 }
56 }
57 else if (DERIVED_T::IS_CENTRE)
58 for (int k = 0; k < ncomp; k++) flux[k] = -psc*0.5*(inco(fac1,k)*porosite(fac1)+inco(fac2,k)*porosite(fac2));
59 else
60 {
61 const int num0_0 = face_amont_princ_(fac1,0), num1_1 = face_amont_princ_(fac2,1);
62 if (DERIVED_T::IS_CENTRE4)
63 {
64 const int ori = orientation(fac1);
65 if ( (num0_0 == -1) || (num1_1== -1) )
66 for (int k = 0; k < ncomp; k++) flux[k] = -psc*0.5*(inco(fac1,k)*porosite(fac1)+inco(fac2,k)*porosite(fac2)); // Schema centre 2
67 else // Schema centre 4
68 {
69 Type_Double vit_0(ncomp),vit_0_0(ncomp),vit_1_1(ncomp),vit_1(ncomp);
70 const double dx = dim_elem_(num_elem,ori), dxam = dim_elem_(elem_(fac1,0),ori), dxav = dim_elem_(elem_(fac2,1),ori);
71 double g1, g2, g3, g4;
72 calcul_g_(dxam,dx,dxav,g1,g2,g3,g4);
73 for (int k = 0; k < ncomp; k++)
74 {
75 vit_0_0[k] = inco(num0_0,k)*porosite(num0_0);
76 vit_0[k] = inco(fac1,k)*porosite(fac1);
77 vit_1[k] = inco(fac2,k)*porosite(fac2);
78 vit_1_1[k] = inco(num1_1,k)*porosite(num1_1);
79 flux[k] = -conv_centre_(psc,vit_0_0[k],vit_0[k],vit_1[k],vit_1_1[k],g1,g2,g3,g4);
80 }
81 }
82 }
83 else // QUICK
84 {
85 if (psc > 0)
86 {
87 if (num0_0 == -1)
88 for (int k=0; k<ncomp; k++) flux[k] = -psc*inco(fac1,k)*porosite(fac1); // Schema amont
89 else
90 {
91 Type_Double vit_0(ncomp), vit_0_0(ncomp), vit_1(ncomp);
92 const int ori = orientation(fac1), elem_amont = elem_(fac1,0);
93 const double dx = dim_elem_(num_elem,ori), dm = dist_elem_period_(elem_amont,num_elem,ori), dxam = dim_elem_(elem_amont,ori);
94 for (int k = 0; k < ncomp; k++)
95 {
96 vit_0[k] = inco(fac1,k)*porosite(fac1);
97 vit_1[k] = inco(fac2,k)*porosite(fac2);
98 vit_0_0[k] = inco(num0_0,k)*porosite(num0_0);
99 flux[k] = -conv_quick_sharp_plus_(psc,vit_0[k],vit_1[k],vit_0_0[k],dx,dm,dxam);
100 }
101 }
102 }
103 else // (psc < 0)
104 {
105 if (num1_1 == -1)
106 for (int k = 0; k < ncomp; k++) flux[k] = -psc*inco(fac2,k)*porosite(fac2); // Schema amont
107 else
108 {
109 Type_Double vit_0(ncomp), vit_1(ncomp), vit_1_1(ncomp);
110 const int ori = orientation(fac2), elem_amont = elem_(fac2,1);
111 const double dx = dim_elem_(num_elem,ori), dm = dist_elem_period_(num_elem,elem_amont,ori), dxam = dim_elem_(elem_amont,ori);
112 for (int k = 0; k < ncomp; k++)
113 {
114 vit_0[k] = inco(fac1,k)*porosite(fac1);
115 vit_1[k] = inco(fac2,k)*porosite(fac2);
116 vit_1_1[k] = inco(num1_1,k)*porosite(num1_1);
117 flux[k] = -conv_quick_sharp_moins_(psc,vit_0[k],vit_1[k],vit_1_1[k],dx,dm,dxam);
118 }
119 }
120 }
121 }
122 }
123}
124
125template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double> inline std::enable_if_t< Arete_Type == Type_Flux_Arete::INTERNE, void>
126Eval_Conv_VDF_Face<DERIVED_T>::flux_arete(const DoubleTab& inco, const DoubleTab* a_r, int fac1, int fac2,int fac3, int fac4, Type_Double& flux) const
127{
128 const int ncomp = flux.size_array();
129 double psc = 0.25 * ((dt_vitesse(fac1)*porosite(fac1)+dt_vitesse(fac2)*porosite(fac2))*(surface(fac1)+surface(fac2)));
130 if (DERIVED_T::IS_AMONT)
131 {
132 for (int k = 0; k < ncomp; k++)
133 {
134 psc = 0.25 * ((dt_vitesse(fac1,k)*porosite(fac1)+dt_vitesse(fac2,k)*porosite(fac2))*(surface(fac1)+surface(fac2)));
135 const int f = psc > 0 ? fac3 : fac4;
136
137 if (a_r)
138 {
139 const int elem = elem_(f, 0), elem2 = elem_(f, 1);
140 const int e = dt_vitesse(f, k) > 0 ? (elem > -1 ? elem : elem2) : (elem2 > -1 ? elem2 : elem);
141 psc *= (*a_r)(e, k);
142 }
143
144 flux[k] = -psc * inco(f, k);
145 }
146 }
147 else if (DERIVED_T::IS_CENTRE)
148 for (int k = 0; k < ncomp; k++) flux[k] = -0.5*(inco(fac3,k)+inco(fac4,k))*psc ;
149 else
150 {
151 const int ori = orientation(fac1), num0_0 = face_amont_conj_(fac3,ori,0), num1_1 = face_amont_conj_(fac4,ori,1);
152 if (DERIVED_T::IS_CENTRE4)
153 {
154 if ( (num0_0 == -1)||(num1_1== -1) )
155 for (int k = 0; k < ncomp; k++) flux[k] = -0.5*(inco(fac3,k)+inco(fac4,k))*psc ; // Schema centre 2 (pas assez de faces)
156 else // Schema Centre 4
157 {
158 Type_Double vit_0(ncomp), vit_0_0(ncomp), vit_1_1(ncomp), vit_1(ncomp);
159 // Inutile de prendre dist_face_period pour dx car fac3 et fac4 ne peuvent etre periodiques (arete interne)
160 const double dx = dist_face_(fac3,fac4,ori), dxam = dist_face_period_(num0_0,fac3,ori), dxav = dist_face_period_(fac4,num1_1,ori);
161 double g1, g2, g3, g4;
162 calcul_g_(dxam,dx,dxav,g1,g2,g3,g4);
163 for (int k = 0; k < ncomp; k++)
164 {
165 vit_0_0[k] = inco(num0_0,k);
166 vit_0[k] = inco(fac3,k);
167 vit_1[k] = inco(fac4,k);
168 vit_1_1[k] = inco(num1_1,k);
169 flux[k] = -conv_centre_(psc,vit_0_0[k],vit_0[k],vit_1[k],vit_1_1[k],g1,g2,g3,g4);
170 }
171 }
172 }
173 else // IS_QUICK
174 {
175 if (psc > 0)
176 {
177 if (num0_0 == -1)
178 for (int k = 0; k < ncomp; k++) flux[k] = -psc*inco(fac3,k); // Schema amont
179 else // Schema quick
180 {
181 Type_Double vit_0(ncomp), vit_0_0(ncomp), vit_1(ncomp);
182 const double dx = dist_face_period_(fac3,fac4,ori), dm = dim_face_(fac3,ori), dxam = dist_face_period_(num0_0,fac3,ori);
183 for (int k = 0; k < ncomp; k++)
184 {
185 vit_0[k] = inco(fac3,k);
186 vit_1[k] = inco(fac4,k);
187 vit_0_0[k] = inco(num0_0,k);
188 flux[k] = -conv_quick_sharp_plus_(psc,vit_0[k],vit_1[k],vit_0_0[k],dx,dm,dxam);
189 }
190 }
191 }
192 else // (psc <= 0)
193 {
194 if (num1_1 == -1)
195 for (int k = 0; k < ncomp; k++) flux[k] = -psc*inco(fac4,k); // Schema amont
196 else // Schema quick
197 {
198 Type_Double vit_0(ncomp), vit_1(ncomp), vit_1_1(ncomp);
199 const double dx = dist_face_period_(fac3,fac4,ori), dm = dim_face_(fac4,ori), dxam = dist_face_period_(fac4,num1_1,ori);
200 for (int k = 0; k < ncomp; k++)
201 {
202 vit_0[k] = inco(fac3,k);
203 vit_1[k] = inco(fac4,k);
204 vit_1_1[k] = inco(num1_1,k);
205 flux[k] = -conv_quick_sharp_moins_(psc,vit_0[k],vit_1[k],vit_1_1[k],dx,dm,dxam);
206 }
207 }
208 }
209 }
210 }
211}
212
213template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double> inline std::enable_if_t< Arete_Type == Type_Flux_Arete::MIXTE, void>
214Eval_Conv_VDF_Face<DERIVED_T>::flux_arete(const DoubleTab& inco, const DoubleTab* a_r, int fac1, int fac2,int fac3, int fac4, Type_Double& flux) const
215{
216 double psc = 0.25*((dt_vitesse(fac1)*porosite(fac1)+dt_vitesse(fac2)*porosite(fac2))*(surface(fac1)+surface(fac2)));
217 if (DERIVED_T::IS_CENTRE)
218 for (int k = 0; k < flux.size_array(); k++) flux[k] = -psc*0.5*(inco(fac3,k)+inco(fac4,k));
219 else
220 {
221 for (int k = 0; k < flux.size_array(); k++)
222 {
223 psc = 0.25 * ((dt_vitesse(fac1, k) * porosite(fac1) + dt_vitesse(fac2, k) * porosite(fac2)) * (surface(fac1) + surface(fac2)));
224 if (psc > 0)
225 {
226 const int elem = elem_(fac3, 0) > 0 ? elem_(fac3, 0) : elem_(fac3, 1);
227 if (a_r) psc *= (*a_r)(elem,k);
228 flux[k] = -psc * inco(fac3, k);
229 }
230 else
231 {
232 const int elem = elem_(fac4, 0) > 0 ? elem_(fac4, 0) : elem_(fac4, 1);
233 if (a_r) psc *= (*a_r)(elem,k);
234 flux[k] = -psc * inco(fac4, k);
235 }
236 }
237 }
238}
239
240template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double>
241inline std::enable_if_t<Arete_Type == Type_Flux_Arete::FLUIDE || Arete_Type == Type_Flux_Arete::NAVIER_FLUIDE || Arete_Type == Type_Flux_Arete::PAROI_FLUIDE, void>
242Eval_Conv_VDF_Face<DERIVED_T>::flux_arete(const DoubleTab& inco, const DoubleTab* a_r, int fac1, int fac2, int fac3, int signe, Type_Double& flux3, Type_Double& flux1_2) const
243{
244 assert(flux3.size_array() == flux1_2.size_array());
245 constexpr bool is_SYM = (Arete_Type == Type_Flux_Arete::NAVIER_FLUIDE);
246 if (DERIVED_T::IS_AXI && is_SYM) return;
247 const int ncomp = flux3.size_array();
248
249 const int pfb = premiere_face_bord(), ori = orientation(fac3), rang1 = DERIVED_T::IS_QUICK ? fac1 : (fac1-pfb), rang2 = DERIVED_T::IS_QUICK ? fac2 :(fac2-pfb); // TODO : FIXME : euh ? pourquoi ca ?
250
251 for (int k = 0; k < ncomp; k++)
252 {
253 double psc = 0.25*((dt_vitesse(fac1,k)*porosite(fac1)+dt_vitesse(fac2,k)*porosite(fac2))*(surface(fac1)+surface(fac2)));
254 if ((psc*signe)>0)
255 {
256 const int elem = elem_(fac3, 0), elem2 = elem_(fac3, 1);
257 const int e = dt_vitesse(fac3, k) > 0 ? (elem > -1 ? elem : elem2) : (elem2 > -1 ? elem2 : elem);
258 const double aa_r = (a_r && DERIVED_T::IS_AMONT) ? (*a_r)(e, k) : 1.0;
259 flux3[k] = -aa_r*inco(fac3,k)*psc ;
260 }
261 else
262 {
263 const int ind = ncomp*ori+k;
264 const double vf1 = Champ_Face_get_val_imp_face_bord_sym(inco,inconnue->temps(),rang1,ind,la_zcl());
265 const double vf2 = Champ_Face_get_val_imp_face_bord_sym(inco,inconnue->temps(),rang2,ind,la_zcl());
266 const int elem = elem_(fac3, 0), elem2 = elem_(fac3, 1);
267 const int e = dt_vitesse(fac3, k) > 0 ? (elem > -1 ? elem : elem2) : (elem2 > -1 ? elem2 : elem);
268 const double aa_r = (a_r && DERIVED_T::IS_AMONT) ? (*a_r)(e, k) : 1.0;
269 flux3[k] = -aa_r * 0.5 * (vf1 + vf2) * psc ;
270 }
271 }
272
273 for (int k = 0; k < ncomp; k++)
274 {
275 double psc = 0.5*dt_vitesse(fac3,k)*surface(fac3)*porosite(fac3);
276 if (psc>0)
277 {
278 const int elem = elem_(fac1, 0), elem2 = elem_(fac1, 1);
279 const int e = dt_vitesse(fac1, k) > 0 ? (elem > -1 ? elem : elem2) : (elem2 > -1 ? elem2 : elem);
280 const double aa_r = (a_r && DERIVED_T::IS_AMONT) ? (*a_r)(e, k) : 1.0;
281 flux1_2[k] = -aa_r * psc * inco(fac1, k);
282 }
283 else
284 {
285 const int elem = elem_(fac2, 0), elem2 = elem_(fac2, 1);
286 const int e = dt_vitesse(fac2, k) > 0 ? (elem > -1 ? elem : elem2) : (elem2 > -1 ? elem2 : elem);
287 const double aa_r = (a_r && DERIVED_T::IS_AMONT) ? (*a_r)(e, k) : 1.0;
288 flux1_2[k] = (DERIVED_T::IS_CENTRE || DERIVED_T::IS_CENTRE4) ? -psc*0.5*(inco(fac1,k)+inco(fac2,k)) : -aa_r * psc * inco(fac2, k);
289 }
290 }
291}
292
293template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double> inline std::enable_if_t< Arete_Type == Type_Flux_Arete::PERIODICITE, void>
294Eval_Conv_VDF_Face<DERIVED_T>::flux_arete(const DoubleTab& inco, const DoubleTab* a_r, int fac1, int fac2 , int fac3, int fac4, Type_Double& flux3_4, Type_Double& flux1_2) const
295{
296 assert(flux3_4.size_array() == flux1_2.size_array());
297 if (DERIVED_T::IS_QUICK) // XXX : LOL
298 {
299 if (DERIVED_T::IS_AXI) return;
300 else
301 {
302 flux_arete < Type_Flux_Arete::INTERNE > (inco, a_r, fac1, fac2, fac3, fac4, flux3_4);
303 flux_arete < Type_Flux_Arete::INTERNE > (inco, a_r, fac3, fac4, fac1, fac2, flux1_2);
304 return;
305 }
306 }
307 const int ncomp = flux3_4.size_array();
308
309 // FIXME : pb_multi !
310 if (ncomp > 1) throw;
311
312 double psc = 0.25*(dt_vitesse(fac1)*porosite(fac1)+dt_vitesse(fac2)*porosite(fac2))*(surface(fac1) +surface(fac2));
313 if (DERIVED_T::IS_CENTRE)
314 for (int k = 0; k < ncomp; k++) flux3_4[k] = -psc*0.5*(inco(fac3,k)+inco(fac4,k));
315 else if (DERIVED_T::IS_CENTRE4)
316 {
317 const int ori = orientation(fac1), num0_0 = face_amont_conj_(fac3,ori,0),num1_1 = face_amont_conj_(fac4,ori,1);
318 if ( (num0_0 == -1)||(num1_1== -1) )
319 for (int k = 0; k < ncomp; k++) flux3_4[k] = -psc*0.5*(inco(fac3,k)+inco(fac4,k)); // Schema centre 2 (pas assez de faces)
320 else // Schema Centre4
321 {
322 Type_Double vit_0(ncomp), vit_0_0(ncomp), vit_1_1(ncomp), vit_1(ncomp);
323 const double dx = dist_face_period_(fac3,fac4,ori), dxam = dist_face_period_(num0_0,fac3,ori), dxav = dist_face_period_(fac4,num1_1,ori);
324 double g1, g2, g3, g4;
325 calcul_g_(dxam,dx,dxav,g1,g2,g3,g4);
326 for (int k = 0; k < ncomp; k++)
327 {
328 vit_0_0[k] = inco(num0_0,k);
329 vit_0[k] = inco(fac3,k);
330 vit_1[k] = inco(fac4,k);
331 vit_1_1[k] = inco(num1_1,k);
332 flux3_4[k] = -conv_centre_(psc,vit_0_0[k],vit_0[k],vit_1[k],vit_1_1[k],g1,g2,g3,g4);
333 }
334 }
335 }
336 else
337 {
338 if (psc>0)
339 for (int k = 0; k < ncomp; k++)
340 {
341// if (a_r && DERIVED_T::IS_AMONT) psc *= (*a_r)(elem_(fac3,0),0); // FIXME
342 flux3_4[k] = -psc*inco(fac3,k);
343 }
344 else for (int k = 0; k < ncomp; k++)
345 {
346// if (a_r && DERIVED_T::IS_AMONT) psc *= (*a_r)(elem_(fac4,0),0); // FIXME
347 flux3_4[k] = -psc*inco(fac4,k);
348 }
349 }
350
351 psc = 0.25*(dt_vitesse(fac3)*porosite(fac3)+dt_vitesse(fac4)*porosite(fac4))*(surface(fac3)+surface(fac4));
352 if (DERIVED_T::IS_CENTRE)
353 for (int k = 0; k < ncomp; k++) flux1_2[k] = -psc*0.5*(inco(fac1,k)+inco(fac2,k));
354 else if (DERIVED_T::IS_CENTRE4)
355 {
356 const int ori = orientation(fac3), num0_0 = face_amont_conj_(fac1,ori,0), num1_1 = face_amont_conj_(fac2,ori,1);
357
358 if ( (num0_0 == -1)||(num1_1== -1) )
359 for (int k=0; k<ncomp; k++) flux1_2[k] = -psc*0.5*(inco(fac1,k)+inco(fac2,k)); // Schema centre 2 (pas assez de faces)
360 else // Schema Centre4
361 {
362 Type_Double vit_0(ncomp), vit_0_0(ncomp), vit_1_1(ncomp), vit_1(ncomp);
363 const double dx = dist_face_period_(fac1,fac2,ori),dxam = dist_face_period_(num0_0,fac1,ori), dxav = dist_face_period_(fac2,num1_1,ori);
364 double g1, g2, g3, g4;
365 calcul_g_(dxam,dx,dxav,g1,g2,g3,g4);
366 for (int k = 0; k < ncomp; k++)
367 {
368 vit_0_0[k] = inco(num0_0,k);
369 vit_0[k] = inco(fac1,k);
370 vit_1[k] = inco(fac2,k);
371 vit_1_1[k]=inco(num1_1,k);
372 flux1_2[k] = -conv_centre_(psc,vit_0_0[k],vit_0[k],vit_1[k],vit_1_1[k],g1,g2,g3,g4);
373 }
374 }
375 }
376 else
377 {
378 if (psc>0)
379 for (int k = 0; k < ncomp; k++)
380 {
381// if (a_r && DERIVED_T::IS_AMONT) psc *= (*a_r)(elem_(fac1,0),0); // FIXME
382 flux1_2[k] = -psc*inco(fac1,k);
383 }
384 else for (int k = 0; k < ncomp; k++)
385 {
386// if (a_r && DERIVED_T::IS_AMONT) psc *= (*a_r)(elem_(fac2,0),0); // FIXME
387 flux1_2[k] = -psc*inco(fac2,k);
388 }
389 }
390}
391
392template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double> inline std::enable_if_t< Arete_Type == Type_Flux_Arete::COIN_FLUIDE, void>
393Eval_Conv_VDF_Face<DERIVED_T>::flux_arete(const DoubleTab& inco, const DoubleTab* a_r, int fac1, int , int fac3, int signe, Type_Double& flux3, Type_Double& flux1_2) const
394{
395 assert(flux3.size_array() == flux1_2.size_array());
396 if (!DERIVED_T::IS_AMONT)
397 {
398 Cerr << "Flux_arete with Type_Flux_Arete::COIN_FLUIDE is only coded for amont scheme !" <<finl;
399 Process::exit();
400 }
401
402 const int ncomp = flux3.size_array();
403
404 // FIXME : pb_multi !
405 if (ncomp > 1) throw;
406
407 double psc = 0.5 * dt_vitesse(fac1) * porosite(fac1) * surface(fac1);
408 if ((psc * signe) > 0)
409 for (int k = 0; k < ncomp; k++)
410 {
411// if (a_r) psc *= (*a_r)(elem_(fac3,0),0); // FIXME
412 flux3[k] = -inco(fac3,k) * psc;
413 }
414 else
415 {
416 const int pfb = premiere_face_bord(), rang1 = (fac1 - pfb), ori = orientation(fac3);
417 for (int k = 0; k < ncomp; k++) flux3[k] = -Champ_Face_get_val_imp_face_bord(inconnue->temps(), rang1, ori, la_zcl()) * psc;
418 }
419
420 psc = 0.5 * dt_vitesse(fac3) * surface(fac3) * porosite(fac3);
421 if (psc > 0)
422 for (int k = 0; k < ncomp; k++)
423 {
424// if (a_r) psc *= (*a_r)(elem_(fac1,0),0); // FIXME
425 flux1_2[k] = -psc * inco(fac1,k);
426 }
427 else
428 {
429 const int pfb = premiere_face_bord(), rang3 = (fac3 - pfb), ori = orientation(fac1);
430 for (int k = 0; k < ncomp; k++) flux1_2[k] = -psc * Champ_Face_get_val_imp_face_bord(inconnue->temps(), rang3, ori, la_zcl());
431 }
432}
433
434/* ************************************** *
435 * ********* POUR L'IMPLICITE ********** *
436 * ************************************** */
437
438template <typename DERIVED_T> template<Type_Flux_Fa7 Fa7_Type, typename Type_Double> inline std::enable_if_t< Fa7_Type == Type_Flux_Fa7::SORTIE_LIBRE, void>
439Eval_Conv_VDF_Face<DERIVED_T>::coeffs_fa7(const DoubleTab* a_r, int face,const Neumann_sortie_libre& la_cl, Type_Double& aii, Type_Double& ajj) const
440{
441 assert(aii.size_array() == ajj.size_array());
442 if (DERIVED_T::IS_CENTRE || DERIVED_T::IS_AXI || DERIVED_T::IS_CENTRE4) return;
443 const int elem1 = elem_(face,0), elem2 = elem_(face,1);
444
445 for (int k = 0; k < aii.size_array(); k++)
446 {
447 double psc = dt_vitesse(face,k) * surface(face) * porosite(face);
448 if (a_r) psc *= (*a_r)((elem1 != -1) ? elem1 : elem2, k);
449 fill_coeffs_proto < Type_Double > (k, psc, psc, aii, ajj);
450 }
451}
452
453template <typename DERIVED_T> template<Type_Flux_Fa7 Fa7_Type, typename Type_Double> inline std::enable_if_t< Fa7_Type == Type_Flux_Fa7::ELEM, void>
454Eval_Conv_VDF_Face<DERIVED_T>::coeffs_fa7(const DoubleTab* a_r, int num_elem, int fac1, int fac2, Type_Double& aii, Type_Double& ajj ) const
455{
456 assert(aii.size_array() == ajj.size_array());
457 if (DERIVED_T::IS_CENTRE || DERIVED_T::IS_AXI || DERIVED_T::IS_CENTRE4) return;
458 else
459 {
460 for (int k = 0; k < aii.size_array(); k++)
461 {
462 double psc = 0.25 * (dt_vitesse(fac1, k) + dt_vitesse(fac2, k)) * (surface(fac1) + surface(fac2));
463 const int f = psc > 0 ? fac1 : fac2;
464
465 if (a_r)
466 {
467 const int elem = elem_(f, 0), elem2 = elem_(f, 1);
468 const int e = dt_vitesse(f, k) > 0 ? (elem > -1 ? elem : elem2) : (elem2 > -1 ? elem2 : elem);
469 psc *= (*a_r)(e, k);
470 }
471
472 const double psc1 = psc * porosite(fac1), psc2 = psc * porosite(fac2);
473 fill_coeffs_proto<Type_Double>(k, psc1, psc2, aii, ajj);
474 }
475 }
476}
477
478template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double> inline
479std::enable_if_t<Arete_Type == Type_Flux_Arete::INTERNE || Arete_Type == Type_Flux_Arete::MIXTE || Arete_Type == Type_Flux_Arete::PERIODICITE, void>
480Eval_Conv_VDF_Face<DERIVED_T>::coeffs_arete(const DoubleTab* a_r, int fac1, int fac2 , int fac3, int fac4, Type_Double& aii, Type_Double& ajj) const
481{
482 assert(aii.size_array() == ajj.size_array());
483 if (DERIVED_T::IS_CENTRE || DERIVED_T::IS_AXI || DERIVED_T::IS_CENTRE4) return;
484 else
485 {
486 for (int k = 0; k < aii.size_array(); k++)
487 {
488 double psc = 0.25 * ((dt_vitesse(fac1, k) * porosite(fac1) + dt_vitesse(fac2, k) * porosite(fac2)) *
489 (surface(fac1) + surface(fac2)));
490 if (a_r)
491 {
492 if (psc > 0)
493 {
494 const int elem = elem_(fac3, 0) > 0 ? elem_(fac3, 0) : elem_(fac3, 1);
495 psc *= (*a_r)(elem, k);
496 }
497 else
498 {
499 const int elem = elem_(fac4, 0) > 0 ? elem_(fac4, 0) : elem_(fac4, 1);
500 psc *= (*a_r)(elem, k);
501 }
502 }
503 fill_coeffs_proto<Type_Double>(k, psc, psc, aii, ajj);
504 }
505 }
506}
507
508template <typename DERIVED_T> template<Type_Flux_Arete Arete_Type, typename Type_Double> inline
509std::enable_if_t<Arete_Type == Type_Flux_Arete::FLUIDE || Arete_Type == Type_Flux_Arete::NAVIER_FLUIDE || Arete_Type == Type_Flux_Arete::PAROI_FLUIDE || Arete_Type == Type_Flux_Arete::COIN_FLUIDE, void>
510Eval_Conv_VDF_Face<DERIVED_T>::coeffs_arete(const DoubleTab* a_r, int fac1, int fac2,int fac3,int signe,Type_Double& aii1_2, Type_Double& aii3_4, Type_Double& ajj1_2) const
511{
512 assert(aii1_2.size_array() == aii3_4.size_array() && aii1_2.size_array() == ajj1_2.size_array());
513 if (DERIVED_T::IS_CENTRE || DERIVED_T::IS_AXI || DERIVED_T::IS_CENTRE4) return;
514
515 constexpr bool is_COIN = (Arete_Type == Type_Flux_Arete::COIN_FLUIDE);
516 const int ncomp = aii1_2.size_array();
517
518 for (int k = 0; k < ncomp; k++)
519 {
520 double psc = is_COIN ? 0.5 * dt_vitesse(fac1,k) * porosite(fac1) * surface(fac1) : 0.25 * ((dt_vitesse(fac1,k) * porosite(fac1) + dt_vitesse(fac2,k) * porosite(fac2)) * (surface(fac1) + surface(fac2)));
521 if ((psc * signe) > 0)
522 {
523 const int elem = elem_(fac3, 0), elem2 = elem_(fac3, 1);
524 const int e = elem > -1 ? elem : elem2;
525 const double aa_r = (a_r && DERIVED_T::IS_AMONT) ? (*a_r)(e, k) : 1.0;
526 aii3_4[k] = aa_r * psc;
527 }
528 else
529 aii3_4[k] = 0.;
530
531 int elem = elem_(fac1, 0), elem2 = elem_(fac1, 1);
532 int e = elem > -1 ? elem : elem2;
533 double aa_r = (a_r && DERIVED_T::IS_AMONT) ? (*a_r)(e, k) : 1.0;
534 const double psc1 = aa_r * 0.5 * dt_vitesse(fac3,k) * surface(fac3) * porosite(fac3);
535
536 elem = elem_(fac2, 0), elem2 = elem_(fac2, 1);
537 e = elem > -1 ? elem : elem2;
538 aa_r = (a_r && DERIVED_T::IS_AMONT) ? (*a_r)(e, k) : 1.0;
539 const double psc2 = aa_r * 0.5 * dt_vitesse(fac3,k) * surface(fac3) * porosite(fac3);
540 fill_coeffs_proto < Type_Double > (k, psc1, psc2, aii1_2, ajj1_2);
541 }
542}
543
544template <typename DERIVED_T> template <typename Type_Double>
545inline void Eval_Conv_VDF_Face<DERIVED_T>::fill_coeffs_proto(const int k, const double psc1, const double psc2, Type_Double& A, Type_Double& B) const
546{
547 if (psc1 > 0)
548 {
549 A[k] = psc1;
550 B[k] = 0.;
551 }
552 else
553 {
554 A[k] = 0.;
555 B[k] = -psc2;
556 }
557}
558
559#endif /* Eval_Conv_VDF_Face_TPP_included */
Eval_Conv_VDF_Face::flux_fa7
std::enable_if_t< Fa7_Type==Type_Flux_Fa7::SORTIE_LIBRE, void > flux_fa7(const DoubleTab &, const DoubleTab *, int, const Neumann_sortie_libre &, int, Type_Double &) const
Definition Eval_Conv_VDF_Face.tpp:24
Eval_Conv_VDF_Face::coeffs_arete
std::enable_if_t< Arete_Type==Type_Flux_Arete::INTERNE||Arete_Type==Type_Flux_Arete::MIXTE||Arete_Type==Type_Flux_Arete::PERIODICITE, void > coeffs_arete(const DoubleTab *, int, int, int, int, Type_Double &, Type_Double &) const
Definition Eval_Conv_VDF_Face.tpp:480
Eval_Conv_VDF_Face::flux_arete
std::enable_if_t< Arete_Type==Type_Flux_Arete::INTERNE, void > flux_arete(const DoubleTab &, const DoubleTab *, int, int, int, int, Type_Double &) const
Definition Eval_Conv_VDF_Face.tpp:126
Eval_Conv_VDF_Face::coeffs_fa7
std::enable_if_t< Fa7_Type==Type_Flux_Fa7::SORTIE_LIBRE, void > coeffs_fa7(const DoubleTab *, int, const Neumann_sortie_libre &, Type_Double &, Type_Double &) const
Definition Eval_Conv_VDF_Face.tpp:439
Neumann_sortie_libre
classe Neumann_sortie_libre Cette classe represente une frontiere ouverte sans vitesse imposee
Definition Neumann_sortie_libre.h:34
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455