TRUST 1.9.8
HPC thermohydraulic platform
Loading...
Searching...
No Matches
distances_VDF.cpp
1/****************************************************************************
2* Copyright (c) 2023, 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 <distances_VDF.h>
17
18void moy_2D_vit(const DoubleVect& vit, int elem, int iori, const Domaine_VDF& domaine, double& u)
19{
20 int num1, num2;
21 if (iori == 0)
22 {
23 num1 = domaine.elem_faces(elem, 1);
24 num2 = domaine.elem_faces(elem, 3);
25 }
26 else if (iori == 1)
27 {
28 num1 = domaine.elem_faces(elem, 0);
29 num2 = domaine.elem_faces(elem, 2);
30 }
31 else
32 {
33 Cerr << "valeur de iori " << iori << " impossible en 2D" << finl;
34 Process::exit();
35 num1 = num2 = -1;
36 }
37 u = 0.5 * (vit(num1) + vit(num2));
38}
39
40double norm_2D_vit(const DoubleVect& vit, int elem, int iori, const Domaine_VDF& domaine, double& u)
41{
42 double v;
43 moy_2D_vit(vit, elem, iori, domaine, v);
44 v = std::fabs(v);
45 if (v == 0)
46 u = 0;
47 else
48 u = 1;
49 return v;
50}
51
52double norm_2D_vit(const DoubleVect& vit, int elem, int iori, const Domaine_VDF& domaine, double u_paroi, double v_paroi, double& u)
53{
54 double vit_paroi;
55 if (iori == 0)
56 vit_paroi = v_paroi;
57 else if (iori == 1)
58 vit_paroi = u_paroi;
59 else
60 {
61 Cerr << "valeur de iori " << iori << " impossible en 2D" << finl;
62 Process::exit();
63 vit_paroi = 0;
64 }
65 double v;
66 double n_v;
67 moy_2D_vit(vit, elem, iori, domaine, v);
68
69 //YB:30/01/04:
70 //Les valeurs du cisaillement parietal sont maintenant signees.
71 //En considerant les valeurs signees de la projection de la vitesse sur le plan parietal
72 v = v - vit_paroi;
73 n_v = std::fabs(v);
74
75 //Fin modif YB
76
77 if (v == 0)
78 u = 0;
79 else if (v > 0)
80 u = 1;
81 else
82 u = -1;
83 return n_v;
84}
85
86void moy_3D_vit(const DoubleVect& vit, int elem, int iori, const Domaine_VDF& domaine, double& val1, double& val2)
87{
88 int num1, num2, num3, num4;
89 if (iori == 0)
90 {
91 num1 = domaine.elem_faces(elem, 1);
92 num2 = domaine.elem_faces(elem, 4);
93 num3 = domaine.elem_faces(elem, 2);
94 num4 = domaine.elem_faces(elem, 5);
95 }
96 else if (iori == 1)
97 {
98 num1 = domaine.elem_faces(elem, 0);
99 num2 = domaine.elem_faces(elem, 3);
100 num3 = domaine.elem_faces(elem, 2);
101 num4 = domaine.elem_faces(elem, 5);
102 }
103 else if (iori == 2)
104 {
105 num1 = domaine.elem_faces(elem, 0);
106 num2 = domaine.elem_faces(elem, 3);
107 num3 = domaine.elem_faces(elem, 1);
108 num4 = domaine.elem_faces(elem, 4);
109 }
110 else
111 {
112 Cerr << "valeur de iori " << iori << " impossible en 3D" << finl;
113 Process::exit();
114 num1 = num2 = num3 = num4 = -1;
115 }
116 val1 = 0.5 * (vit(num1) + vit(num2));
117 val2 = 0.5 * (vit(num3) + vit(num4));
118}
119
120double norm_3D_vit(const DoubleVect& vit, int elem, int iori, const Domaine_VDF& domaine, double& val1, double& val2)
121{
122 moy_3D_vit(vit, elem, iori, domaine, val1, val2);
123 double v1 = std::fabs(val1);
124 double v2 = std::fabs(val2);
125 double norm_vit = sqrt(v1 * v1 + v2 * v2);
126 val1 = v1 / (norm_vit + DMINFLOAT);
127 val2 = v2 / (norm_vit + DMINFLOAT);
128 return norm_vit;
129}
130
131double norm_3D_vit(const DoubleVect& vit, int elem, int iori, const Domaine_VDF& domaine, double u_paroi, double v_paroi, double w_paroi, double& val1, double& val2)
132{
133 double v1, v2, norm_vit;
134 moy_3D_vit(vit, elem, iori, domaine, val1, val2);
135 if (iori == 0)
136 {
137 v1 = val1 - v_paroi; // EB 28/08/25 : for a wall of normal x, val1 is the velocity in y direction and val2 in z direction
138 v2 = val2 - w_paroi;
139 }
140 else if (iori == 1)
141 {
142 v1 = val1 - u_paroi;
143 v2 = val2 - w_paroi;
144 }
145 else if (iori == 2)
146 {
147 v1 = val1 - u_paroi;
148 v2 = val2 - v_paroi;
149 }
150 //Fin modif YB
151
152 else
153 {
154 Cerr << "valeur de iori " << iori << " impossible en 3D" << finl;
155 Process::exit();
156 v1 = v2 = 0;
157 }
158 norm_vit = sqrt(v1 * v1 + v2 * v2);
159 val1 = v1 / (norm_vit + DMINFLOAT);
160 val2 = v2 / (norm_vit + DMINFLOAT);
161 return norm_vit;
162}
163
164double norm_vit(const DoubleVect& vit, int elem, int ori, const Domaine_VDF& domaine, const ArrOfDouble& vit_paroi, ArrOfDouble& val)
165{
166 if (Objet_U::dimension == 3)
167 return norm_3D_vit(vit, elem, ori, domaine, vit_paroi[0], vit_paroi[1], vit_paroi[2], val[0], val[1]);
168 else
169 return norm_2D_vit(vit, elem, ori, domaine, vit_paroi[0], vit_paroi[1], val[0]);
170}
171
172
173void calcul_interne2D(int num_elem, int elx0, int elx1, int ely0, int ely1, const Domaine_VDF& domaine_VDF, const DoubleTab& val, DoubleTab& rot)
174{
175 const IntTab& elem_faces = domaine_VDF.elem_faces();
176 double delta_x_0, delta_x_1, delta_y_0, delta_y_1;
177 double delta_x, delta_y;
178 double deriv_vx, deriv_uy;
179
180 int N = val.line_size(), n;
181
182 deriv_vx = 0;
183 deriv_uy = 0;
184
185 delta_x_0 = domaine_VDF.dist_elem_period(num_elem, elx0, 0);
186 delta_x_1 = domaine_VDF.dist_elem_period(elx1, num_elem, 0);
187 delta_y_0 = domaine_VDF.dist_elem_period(num_elem, ely0, 1);
188 delta_y_1 = domaine_VDF.dist_elem_period(ely1, num_elem, 1);
189
190 delta_x = (delta_x_1 - delta_x_0) * (delta_x_1 + delta_x_0) / (delta_x_1 * delta_x_0);
191 delta_y = (delta_y_1 - delta_y_0) * (delta_y_1 + delta_y_0) / (delta_y_1 * delta_y_0);
192
193 for (n=0; n<N; n++)
194 {
195 deriv_vx = (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 1), n) + delta_x * val(elem_faces(num_elem, 1), n) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 1), n));
196 deriv_vx += (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 3), n) + delta_x * val(elem_faces(num_elem, 3), n) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 3), n));
197 deriv_vx *= 0.5 / (delta_x_0 + delta_x_1);
198
199 deriv_uy = (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 0), n) + delta_y * val(elem_faces(num_elem, 0), n) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 0), n));
200 deriv_uy += (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 2), n) + delta_y * val(elem_faces(num_elem, 2), n) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 2), n));
201 deriv_uy *= 0.5 / (delta_y_0 + delta_y_1);
202
203 rot(num_elem,n) = deriv_vx - deriv_uy;
204 }
205}
206
207void calcul_bord2D(int num_elem, int elx0, int elx1, int ely0, int ely1, const Domaine_VDF& domaine_VDF, const DoubleTab& val, DoubleTab& rot)
208{
209 const IntTab& elem_faces = domaine_VDF.elem_faces();
210 double delta_x_0, delta_x_1, delta_y_0, delta_y_1;
211 double delta_x, delta_y;
212 double deriv_vx, deriv_uy;
213
214 int N = val.line_size(), n;
215
216 deriv_vx = 0;
217 deriv_uy = 0;
218
219 for (n=0; n<N; n++)
220 {
221 // Traitement des elements bord
222 if ((elx0 == -1) || (elx1 == -1))
223 {
224 if (elx0 == -1)
225 {
226 if (elx1 != -1)
227 {
228 delta_x_1 = domaine_VDF.dist_elem_period(elx1, num_elem, 0);
229 deriv_vx = (val(elem_faces(elx1, 1), n) - val(elem_faces(num_elem, 1), n) + val(elem_faces(elx1, 3), n) - val(elem_faces(num_elem, 3), n));
230 deriv_vx *= 0.5 / delta_x_1;
231
232 }
233 else
234 deriv_vx = 0;
235 }
236 else // elx1 = -1 et elx0 != -1
237 {
238 delta_x_0 = domaine_VDF.dist_elem_period(num_elem, elx0, 0);
239 deriv_vx = (val(elem_faces(num_elem, 1), n) - val(elem_faces(elx0, 1), n) + val(elem_faces(num_elem, 3), n) - val(elem_faces(elx0, 3), n));
240 deriv_vx *= 0.5 / delta_x_0;
241
242 }
243 }
244 else // elx0 != -1 et elx1 != -1
245 {
246 delta_x_0 = domaine_VDF.dist_elem_period(num_elem, elx0, 0);
247 delta_x_1 = domaine_VDF.dist_elem_period(elx1, num_elem, 0);
248 delta_x = (delta_x_1 - delta_x_0) * (delta_x_1 + delta_x_0) / (delta_x_1 * delta_x_0);
249 deriv_vx = (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 1), n) + delta_x * val(elem_faces(num_elem, 1), n) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 1), n));
250 deriv_vx += (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 3), n) + delta_x * val(elem_faces(num_elem, 3), n) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 3), n));
251 deriv_vx *= 0.5 / (delta_x_0 + delta_x_1);
252
253 }
254
255 if ((ely0 == -1) || (ely1 == -1))
256 {
257 if (ely0 == -1)
258 {
259 if (ely1 != -1)
260 {
261 delta_y_1 = domaine_VDF.dist_elem_period(ely1, num_elem, 1);
262 deriv_uy = (val(elem_faces(ely1, 0), n) - val(elem_faces(num_elem, 0), n) + val(elem_faces(ely1, 2), n) - val(elem_faces(num_elem, 2), n));
263 deriv_uy *= 0.5 / delta_y_1;
264
265 }
266 else
267 deriv_uy = 0;
268 }
269 else // ely1 = -1 et ely0 != -1
270 {
271 delta_y_0 = domaine_VDF.dist_elem_period(num_elem, ely0, 1);
272 deriv_uy = (val(elem_faces(num_elem, 0), n) - val(elem_faces(ely0, 0), n) + val(elem_faces(num_elem, 2), n) - val(elem_faces(ely0, 2), n));
273 deriv_uy *= 0.5 / delta_y_0;
274
275 }
276 }
277 else // ely0 != -1 et ely1 != -1
278 {
279 delta_y_0 = domaine_VDF.dist_elem_period(num_elem, ely0, 1);
280 delta_y_1 = domaine_VDF.dist_elem_period(ely1, num_elem, 1);
281 delta_y = (delta_y_1 - delta_y_0) * (delta_y_1 + delta_y_0) / (delta_y_1 * delta_y_0);
282
283 deriv_uy = (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 0), n) + delta_y * val(elem_faces(num_elem, 0), n) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 0), n));
284 deriv_uy += (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 2), n) + delta_y * val(elem_faces(num_elem, 2), n) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 2), n));
285 deriv_uy *= 0.5 / (delta_y_0 + delta_y_1);
286
287 }
288
289 rot(num_elem,n) = deriv_vx - deriv_uy;
290 }
291}
292
293void calrotord2centelemdim2(DoubleTab& rot, const DoubleTab& val, const Domaine_VDF& domaine_VDF, int nb_elem, const IntTab& face_voisins, const IntTab& elem_faces)
294{
295 if (rot.dimension(0) != nb_elem)
296 rot.resize(nb_elem);
297 int elx0, elx1, ely0, ely1;
298
299 for (int num_elem = 0; num_elem < nb_elem; num_elem++)
300 {
301 elx0 = face_voisins(elem_faces(num_elem, 0), 0);
302 elx1 = face_voisins(elem_faces(num_elem, 2), 1);
303 ely0 = face_voisins(elem_faces(num_elem, 1), 0);
304 ely1 = face_voisins(elem_faces(num_elem, 3), 1);
305
306 if ((elx0 != -1) && (elx1 != -1) && (ely0 != -1) && (ely1 != -1))
307 // Cas d'un element interne
308
309 calcul_interne2D(num_elem, elx0, elx1, ely0, ely1, domaine_VDF, val, rot);
310 else
311 calcul_bord2D(num_elem, elx0, elx1, ely0, ely1, domaine_VDF, val, rot);
312 }
313}
314
315// Traitement des elements internes
316void calcul_interne3D(int num_elem, int elx0, int elx1, int ely0, int ely1, int elz0, int elz1, const Domaine_VDF& domaine_VDF, const DoubleTab& val, DoubleTab& rot)
317{
318 const IntTab& elem_faces = domaine_VDF.elem_faces();
319 double delta_x_0, delta_x_1, delta_y_0, delta_y_1, delta_z_0, delta_z_1;
320 double delta_x, delta_y, delta_z;
321 double deriv_wy, deriv_vz, deriv_uz, deriv_wx, deriv_vx, deriv_uy;
322
323 int N = val.line_size();
324
325 deriv_wy = 0;
326 deriv_vz = 0;
327 deriv_uz = 0;
328 deriv_wx = 0;
329 deriv_vx = 0;
330 deriv_uy = 0;
331
332 delta_x_0 = domaine_VDF.dist_elem_period(num_elem, elx0, 0);
333 delta_x_1 = domaine_VDF.dist_elem_period(elx1, num_elem, 0);
334 delta_y_0 = domaine_VDF.dist_elem_period(num_elem, ely0, 1);
335 delta_y_1 = domaine_VDF.dist_elem_period(ely1, num_elem, 1);
336 delta_z_0 = domaine_VDF.dist_elem_period(num_elem, elz0, 2);
337 delta_z_1 = domaine_VDF.dist_elem_period(elz1, num_elem, 2);
338
339 delta_x = (delta_x_1 - delta_x_0) * (delta_x_1 + delta_x_0) / (delta_x_1 * delta_x_0);
340 delta_y = (delta_y_1 - delta_y_0) * (delta_y_1 + delta_y_0) / (delta_y_1 * delta_y_0);
341 delta_z = (delta_z_1 - delta_z_0) * (delta_z_1 + delta_z_0) / (delta_z_1 * delta_z_0);
342
343 for (int n=0; n<N; n++)
344 {
345 deriv_vz = (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 1), n) + delta_z * val(elem_faces(num_elem, 1), n) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 1), n));
346 deriv_vz += (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 4), n) + delta_z * val(elem_faces(num_elem, 4), n) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 4), n));
347 deriv_vz *= 0.5 / (delta_z_0 + delta_z_1);
348
349 deriv_wy = (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 2), n) + delta_y * val(elem_faces(num_elem, 2), n) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 2), n));
350 deriv_wy += (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 5), n) + delta_y * val(elem_faces(num_elem, 5), n) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 5), n));
351 deriv_wy *= 0.5 / (delta_y_0 + delta_y_1);
352
353 deriv_uz = (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 0), n) + delta_z * val(elem_faces(num_elem, 0), n) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 0), n));
354 deriv_uz += (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 3), n) + delta_z * val(elem_faces(num_elem, 3), n) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 3), n));
355 deriv_uz *= 0.5 / (delta_z_0 + delta_z_1);
356
357 deriv_wx = (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 2), n) + delta_x * val(elem_faces(num_elem, 2), n) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 2), n));
358 deriv_wx += (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 5), n) + delta_x * val(elem_faces(num_elem, 5), n) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 5), n));
359 deriv_wx *= 0.5 / (delta_x_0 + delta_x_1);
360
361 deriv_vx = (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 1), n) + delta_x * val(elem_faces(num_elem, 1), n) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 1), n));
362 deriv_vx += (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 4), n) + delta_x * val(elem_faces(num_elem, 4), n) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 4), n));
363 deriv_vx *= 0.5 / (delta_x_0 + delta_x_1);
364
365 deriv_uy = (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 0), n) + delta_y * val(elem_faces(num_elem, 0), n) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 0), n));
366 deriv_uy += (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 3), n) + delta_y * val(elem_faces(num_elem, 3), n) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 3), n));
367 deriv_uy *= 0.5 / (delta_y_0 + delta_y_1);
368
369 rot(num_elem, N*0+n) = deriv_wy - deriv_vz;
370 rot(num_elem, N*1+n) = deriv_uz - deriv_wx;
371 rot(num_elem, N*2+n) = deriv_vx - deriv_uy;
372 }
373}
374
375// Traitement des elements bord
376void calcul_bord3D(int num_elem, int elx0, int elx1, int ely0, int ely1, int elz0, int elz1, const Domaine_VDF& domaine_VDF, const DoubleTab& val, DoubleTab& rot)
377{
378 const IntTab& elem_faces = domaine_VDF.elem_faces();
379 double delta_x_0, delta_x_1, delta_y_0, delta_y_1, delta_z_0, delta_z_1;
380 double delta_x, delta_y, delta_z;
381 double deriv_wy, deriv_vz, deriv_uz, deriv_wx, deriv_vx, deriv_uy;
382
383 int N = val.line_size();
384
385 deriv_wy = 0;
386 deriv_vz = 0;
387 deriv_uz = 0;
388 deriv_wx = 0;
389 deriv_vx = 0;
390 deriv_uy = 0;
391
392 for (int n=0; n<N; n++)
393 {
394
395 // Traitement des elements bord
396
397 if ((elx0 == -1) || (elx1 == -1))
398 {
399 if (elx0 == -1)
400 {
401 if (elx1 != -1)
402 {
403 delta_x_1 = domaine_VDF.dist_elem_period(elx1, num_elem, 0);
404 deriv_vx = (val(elem_faces(elx1, 1), n) - val(elem_faces(num_elem, 1), n) + val(elem_faces(elx1, 4), n) - val(elem_faces(num_elem, 4), n));
405 deriv_vx *= 0.5 / delta_x_1;
406 deriv_wx = (val(elem_faces(elx1, 2), n) - val(elem_faces(num_elem, 2), n) + val(elem_faces(elx1, 5), n) - val(elem_faces(num_elem, 5), n));
407 deriv_wx *= 0.5 / delta_x_1;
408 }
409 else
410 {
411 deriv_vx = 0;
412 deriv_wx = 0;
413 }
414 }
415 else // elx1 = -1 et elx0 != -1
416 {
417 delta_x_0 = domaine_VDF.dist_elem_period(num_elem, elx0, 0);
418 deriv_vx = (val(elem_faces(num_elem, 1), n) - val(elem_faces(elx0, 1), n) + val(elem_faces(num_elem, 4), n) - val(elem_faces(elx0, 4), n));
419 deriv_vx *= 0.5 / delta_x_0;
420 deriv_wx = (val(elem_faces(num_elem, 2), n) - val(elem_faces(elx0, 2), n) + val(elem_faces(num_elem, 5), n) - val(elem_faces(elx0, 5), n));
421 deriv_wx *= 0.5 / delta_x_0;
422 }
423
424 if ((ely0 == -1) || (ely1 == -1))
425 {
426 if (ely0 == -1)
427 {
428 if (ely1 != -1)
429 {
430 delta_y_1 = domaine_VDF.dist_elem_period(ely1, num_elem, 1);
431 deriv_uy = (val(elem_faces(ely1, 0), n) - val(elem_faces(num_elem, 0), n) + val(elem_faces(ely1, 3), n) - val(elem_faces(num_elem, 3), n));
432 deriv_uy *= 0.5 / delta_y_1;
433 deriv_wy = (val(elem_faces(ely1, 2), n) - val(elem_faces(num_elem, 2), n) + val(elem_faces(ely1, 5), n) - val(elem_faces(num_elem, 5), n));
434 deriv_wy *= 0.5 / delta_y_1;
435 }
436 else
437 {
438 deriv_uy = 0;
439 deriv_wy = 0;
440 }
441 }
442 else // ely1 = -1 et ely0 != -1
443 {
444 delta_y_0 = domaine_VDF.dist_elem_period(num_elem, ely0, 1);
445 deriv_uy = (val(elem_faces(num_elem, 0), n) - val(elem_faces(ely0, 0), n) + val(elem_faces(num_elem, 3), n) - val(elem_faces(ely0, 3), n));
446 deriv_uy *= 0.5 / delta_y_0;
447 deriv_wy = (val(elem_faces(num_elem, 2), n) - val(elem_faces(ely0, 2), n) + val(elem_faces(num_elem, 5), n) - val(elem_faces(ely0, 5), n));
448 deriv_wy *= 0.5 / delta_y_0;
449 }
450 }
451 else
452 {
453 delta_y_0 = domaine_VDF.dist_elem_period(num_elem, ely0, 1);
454 delta_y_1 = domaine_VDF.dist_elem_period(ely1, num_elem, 1);
455 delta_y = (delta_y_1 - delta_y_0) * (delta_y_1 + delta_y_0) / (delta_y_1 * delta_y_0);
456
457 deriv_uy = (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 0), n) + delta_y * val(elem_faces(num_elem, 0), n) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 0), n));
458 deriv_uy += (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 3), n) + delta_y * val(elem_faces(num_elem, 3), n) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 3), n));
459 deriv_uy *= 0.5 / (delta_y_0 + delta_y_1);
460
461 deriv_wy = (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 2), n) + delta_y * val(elem_faces(num_elem, 2), n) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 2), n));
462 deriv_wy += (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 5), n) + delta_y * val(elem_faces(num_elem, 5), n) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 5), n));
463 deriv_wy *= 0.5 / (delta_y_0 + delta_y_1);
464 }
465
466 if ((elz0 == -1) || (elz1 == -1))
467 {
468 if (elz0 == -1)
469 {
470 if (elz1 != -1)
471 {
472 delta_z_1 = domaine_VDF.dist_elem_period(elz1, num_elem, 2);
473 deriv_uz = (val(elem_faces(elz1, 0), n) - val(elem_faces(num_elem, 0), n) + val(elem_faces(elz1, 3), n) - val(elem_faces(num_elem, 3), n));
474 deriv_uz *= 0.5 / delta_z_1;
475 deriv_vz = (val(elem_faces(elz1, 1), n) - val(elem_faces(num_elem, 1), n) + val(elem_faces(elz1, 4), n) - val(elem_faces(num_elem, 4), n));
476 deriv_vz *= 0.5 / delta_z_1;
477 }
478 else
479 {
480 deriv_uz = 0;
481 deriv_vz = 0;
482 }
483 }
484 else // elz1 = -1 et elz0 != -1
485 {
486 delta_z_0 = domaine_VDF.dist_elem_period(num_elem, elz0, 2);
487 deriv_uz = (val(elem_faces(num_elem, 0), n) - val(elem_faces(elz0, 0), n) + val(elem_faces(num_elem, 3), n) - val(elem_faces(elz0, 3), n));
488 deriv_uz *= 0.5 / delta_z_0;
489 deriv_vz = (val(elem_faces(num_elem, 1), n) - val(elem_faces(elz0, 1), n) + val(elem_faces(num_elem, 4), n) - val(elem_faces(elz0, 4), n));
490 deriv_vz *= 0.5 / delta_z_0;
491 }
492 }
493 else
494 {
495 delta_z_0 = domaine_VDF.dist_elem_period(num_elem, elz0, 2);
496 delta_z_1 = domaine_VDF.dist_elem_period(elz1, num_elem, 2);
497 delta_z = (delta_z_1 - delta_z_0) * (delta_z_1 + delta_z_0) / (delta_z_1 * delta_z_0);
498
499 deriv_uz = (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 0), n) + delta_z * val(elem_faces(num_elem, 0), n) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 0), n));
500 deriv_uz += (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 3), n) + delta_z * val(elem_faces(num_elem, 3), n) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 3), n));
501 deriv_uz *= 0.5 / (delta_z_0 + delta_z_1);
502
503 deriv_vz = (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 1), n) + delta_z * val(elem_faces(num_elem, 1), n) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 1), n));
504 deriv_vz += (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 4), n) + delta_z * val(elem_faces(num_elem, 4), n) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 4), n));
505 deriv_vz *= 0.5 / (delta_z_0 + delta_z_1);
506 }
507
508 }
509 else if ((ely0 == -1) || (ely1 == -1))
510 {
511
512 if (ely0 == -1)
513 {
514 if (ely1 != -1)
515 {
516 delta_y_1 = domaine_VDF.dist_elem_period(ely1, num_elem, 1);
517 deriv_uy = (val(elem_faces(ely1, 0), n) - val(elem_faces(num_elem, 0), n) + val(elem_faces(ely1, 3), n) - val(elem_faces(num_elem, 3), n));
518 deriv_uy *= 0.5 / delta_y_1;
519 deriv_wy = (val(elem_faces(ely1, 2), n) - val(elem_faces(num_elem, 2), n) + val(elem_faces(ely1, 5), n) - val(elem_faces(num_elem, 5), n));
520 deriv_wy *= 0.5 / delta_y_1;
521 }
522 else
523 {
524 deriv_uy = 0;
525 deriv_wy = 0;
526 }
527 }
528 else // ely1 = -1 et ely0 != -1
529 {
530 delta_y_0 = domaine_VDF.dist_elem_period(num_elem, ely0, 1);
531 deriv_uy = (val(elem_faces(num_elem, 0), n) - val(elem_faces(ely0, 0), n) + val(elem_faces(num_elem, 3), n) - val(elem_faces(ely0, 3), n));
532 deriv_uy *= 0.5 / delta_y_0;
533 deriv_wy = (val(elem_faces(num_elem, 2), n) - val(elem_faces(ely0, 2), n) + val(elem_faces(num_elem, 5), n) - val(elem_faces(ely0, 5), n));
534 deriv_wy *= 0.5 / delta_y_0;
535 }
536
537 delta_x_0 = domaine_VDF.dist_elem_period(num_elem, elx0, 0);
538 delta_x_1 = domaine_VDF.dist_elem_period(elx1, num_elem, 0);
539 delta_x = (delta_x_1 - delta_x_0) * (delta_x_1 + delta_x_0) / (delta_x_1 * delta_x_0);
540
541 deriv_wx = (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 2), n) + delta_x * val(elem_faces(num_elem, 2), n) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 2), n));
542 deriv_wx += (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 5), n) + delta_x * val(elem_faces(num_elem, 5), n) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 5), n));
543 deriv_wx *= 0.5 / (delta_x_0 + delta_x_1);
544
545 deriv_vx = (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 1), n) + delta_x * val(elem_faces(num_elem, 1), n) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 1), n));
546 deriv_vx += (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 4), n) + delta_x * val(elem_faces(num_elem, 4), n) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 4), n));
547 deriv_vx *= 0.5 / (delta_x_0 + delta_x_1);
548
549 if ((elz0 == -1) || (elz1 == -1))
550 {
551 if (elz0 == -1)
552 {
553 if (elz1 != -1)
554 {
555 delta_z_1 = domaine_VDF.dist_elem_period(elz1, num_elem, 2);
556 deriv_uz = (val(elem_faces(elz1, 0), n) - val(elem_faces(num_elem, 0), n) + val(elem_faces(elz1, 3), n) - val(elem_faces(num_elem, 3), n));
557 deriv_uz *= 0.5 / delta_z_1;
558 deriv_vz = (val(elem_faces(elz1, 1), n) - val(elem_faces(num_elem, 1), n) + val(elem_faces(elz1, 4), n) - val(elem_faces(num_elem, 4), n));
559 deriv_vz *= 0.5 / delta_z_1;
560 }
561 else
562 {
563 deriv_uz = 0;
564 deriv_vz = 0;
565 }
566 }
567 else
568 {
569 delta_z_0 = domaine_VDF.dist_elem_period(num_elem, elz0, 2);
570 deriv_uz = (val(elem_faces(num_elem, 0), n) - val(elem_faces(elz0, 0), n) + val(elem_faces(num_elem, 3), n) - val(elem_faces(elz0, 3), n));
571 deriv_uz *= 0.5 / delta_z_0;
572 deriv_vz = (val(elem_faces(num_elem, 1), n) - val(elem_faces(elz0, 1), n) + val(elem_faces(num_elem, 4), n) - val(elem_faces(elz0, 4), n));
573 deriv_vz *= 0.5 / delta_z_0;
574 }
575 }
576 else
577 {
578 delta_z_0 = domaine_VDF.dist_elem_period(num_elem, elz0, 2);
579 delta_z_1 = domaine_VDF.dist_elem_period(elz1, num_elem, 2);
580 delta_z = (delta_z_1 - delta_z_0) * (delta_z_1 + delta_z_0) / (delta_z_1 * delta_z_0);
581
582 deriv_uz = (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 0), n) + delta_z * val(elem_faces(num_elem, 0), n) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 0), n));
583 deriv_uz += (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 3), n) + delta_z * val(elem_faces(num_elem, 3), n) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 3), n));
584 deriv_uz *= 0.5 / (delta_z_0 + delta_z_1);
585
586 deriv_vz = (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 1), n) + delta_z * val(elem_faces(num_elem, 1), n) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 1), n));
587 deriv_vz += (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 4), n) + delta_z * val(elem_faces(num_elem, 4), n) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 4), n));
588 deriv_vz *= 0.5 / (delta_z_0 + delta_z_1);
589 }
590 }
591 else if ((elz0 == -1) || (elz1 == -1))
592 {
593 if (elz0 == -1)
594 {
595 if (elz1 != -1)
596 {
597 delta_z_1 = domaine_VDF.dist_elem_period(elz1, num_elem, 2);
598 deriv_uz = (val(elem_faces(elz1, 0), n) - val(elem_faces(num_elem, 0), n) + val(elem_faces(elz1, 3), n) - val(elem_faces(num_elem, 3), n));
599 deriv_uz *= 0.5 / delta_z_1;
600 deriv_vz = (val(elem_faces(elz1, 1), n) - val(elem_faces(num_elem, 1), n) + val(elem_faces(elz1, 4), n) - val(elem_faces(num_elem, 4), n));
601 deriv_vz *= 0.5 / delta_z_1;
602 }
603 else
604 {
605 deriv_uz = 0;
606 deriv_vz = 0;
607 }
608 }
609 else // elz1 = -1 et elz0 != -1
610 {
611 delta_z_0 = domaine_VDF.dist_elem_period(num_elem, elz0, 2);
612 deriv_uz = (val(elem_faces(num_elem, 0), n) - val(elem_faces(elz0, 0), n) + val(elem_faces(num_elem, 3), n) - val(elem_faces(elz0, 3), n));
613 deriv_uz *= 0.5 / delta_z_0;
614 deriv_vz = (val(elem_faces(num_elem, 1), n) - val(elem_faces(elz0, 1), n) + val(elem_faces(num_elem, 4), n) - val(elem_faces(elz0, 4), n));
615 deriv_vz *= 0.5 / delta_z_0;
616 }
617
618 delta_x_0 = domaine_VDF.dist_elem_period(num_elem, elx0, 0);
619 delta_x_1 = domaine_VDF.dist_elem_period(elx1, num_elem, 0);
620 delta_x = (delta_x_1 - delta_x_0) * (delta_x_1 + delta_x_0) / (delta_x_1 * delta_x_0);
621
622 deriv_wx = (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 2), n) + delta_x * val(elem_faces(num_elem, 2), n) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 2), n));
623 deriv_wx += (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 5), n) + delta_x * val(elem_faces(num_elem, 5), n) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 5), n));
624 deriv_wx *= 0.5 / (delta_x_0 + delta_x_1);
625
626 deriv_vx = (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 1), n) + delta_x * val(elem_faces(num_elem, 1), n) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 1), n));
627 deriv_vx += (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 4), n) + delta_x * val(elem_faces(num_elem, 4), n) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 4), n));
628 deriv_vx *= 0.5 / (delta_x_0 + delta_x_1);
629
630 delta_y_0 = domaine_VDF.dist_elem_period(num_elem, ely0, 1);
631 delta_y_1 = domaine_VDF.dist_elem_period(ely1, num_elem, 1);
632 delta_y = (delta_y_1 - delta_y_0) * (delta_y_1 + delta_y_0) / (delta_y_1 * delta_y_0);
633
634 deriv_wy = (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 2), n) + delta_y * val(elem_faces(num_elem, 2), n) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 2), n));
635 deriv_wy += (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 5), n) + delta_y * val(elem_faces(num_elem, 5), n) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 5), n));
636 deriv_wy *= 0.5 / (delta_y_0 + delta_y_1);
637
638 deriv_uy = (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 0), n) + delta_y * val(elem_faces(num_elem, 0), n) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 0), n));
639 deriv_uy += (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 3), n) + delta_y * val(elem_faces(num_elem, 3), n) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 3), n));
640 deriv_uy *= 0.5 / (delta_y_0 + delta_y_1);
641 }
642
643 rot(num_elem, N*0+n) = deriv_wy - deriv_vz;
644 rot(num_elem, N*1+n) = deriv_uz - deriv_wx;
645 rot(num_elem, N*2+n) = deriv_vx - deriv_uy;
646 }
647}
648
649void calrotord2centelemdim3(DoubleTab& rot, const DoubleTab& val, const Domaine_VDF& domaine_VDF, int nb_elem, const IntTab& face_voisins, const IntTab& elem_faces)
650{
651 if (rot.dimension(0) != nb_elem)
652 rot.resize(nb_elem, 3);
653 int elx0, elx1, ely0, ely1, elz0, elz1;
654
655 for (int num_elem = 0; num_elem < nb_elem; num_elem++)
656 {
657
658 elx0 = face_voisins(elem_faces(num_elem, 0), 0);
659 elx1 = face_voisins(elem_faces(num_elem, 3), 1);
660 ely0 = face_voisins(elem_faces(num_elem, 1), 0);
661 ely1 = face_voisins(elem_faces(num_elem, 4), 1);
662 elz0 = face_voisins(elem_faces(num_elem, 2), 0);
663 elz1 = face_voisins(elem_faces(num_elem, 5), 1);
664
665 if ((elx0 != -1) && (elx1 != -1) && (ely0 != -1) && (ely1 != -1) && (elz0 != -1) && (elz1 != -1))
666 // Cas d'un element interne
667
668 calcul_interne3D(num_elem, elx0, elx1, ely0, ely1, elz0, elz1, domaine_VDF, val, rot);
669 else
670 calcul_bord3D(num_elem, elx0, elx1, ely0, ely1, elz0, elz1, domaine_VDF, val, rot);
671 }
672}
673
674// Calcul du produit scalaire du tenseur des vitesses de deformation
675// en coordonnees cartesiennes : calcul 2D puis 3D.
676// Traitement des elements internes
677
678// Cas 2D
679void calcul_dscald_interne2D(int num_elem, int elx0, int elx1, int ely0, int ely1, const Domaine_VDF& domaine_VDF, const DoubleTab& val, DoubleTab& dscald)
680{
681 const IntTab& elem_faces = domaine_VDF.elem_faces();
682 double delta_x_0, delta_x_1, delta_y_0, delta_y_1;
683 double delta_x, delta_y;
684 double delta_xbis, delta_ybis;
685 double deriv_ux, deriv_uy, deriv_vx, deriv_vy;
686
687 deriv_ux = 0;
688 deriv_uy = 0;
689 deriv_vx = 0;
690 deriv_vy = 0;
691
692 delta_x_0 = domaine_VDF.dist_elem(num_elem, elx0, 0);
693 delta_x_1 = domaine_VDF.dist_elem(elx1, num_elem, 0);
694 delta_y_0 = domaine_VDF.dist_elem(num_elem, ely0, 1);
695 delta_y_1 = domaine_VDF.dist_elem(ely1, num_elem, 1);
696
697 delta_x = (delta_x_1 - delta_x_0) * (delta_x_1 + delta_x_0) / (delta_x_1 * delta_x_0);
698 delta_y = (delta_y_1 - delta_y_0) * (delta_y_1 + delta_y_0) / (delta_y_1 * delta_y_0);
699
700 delta_xbis = domaine_VDF.dim_elem(num_elem, 0);
701 delta_ybis = domaine_VDF.dim_elem(num_elem, 1);
702
703 deriv_vx = (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 1)) + delta_x * val(elem_faces(num_elem, 1)) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 1)));
704 deriv_vx += (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 3)) + delta_x * val(elem_faces(num_elem, 3)) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 3)));
705 deriv_vx *= 0.5 / (delta_x_0 + delta_x_1);
706
707 deriv_uy = (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 0)) + delta_y * val(elem_faces(num_elem, 0)) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 0)));
708 deriv_uy += (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 2)) + delta_y * val(elem_faces(num_elem, 2)) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 2)));
709 deriv_uy *= 0.5 / (delta_y_0 + delta_y_1);
710
711 deriv_ux = (1. / delta_xbis) * (val(elem_faces(num_elem, 0)) - val(elem_faces(num_elem, 2)));
712 deriv_vy = (1. / delta_ybis) * (val(elem_faces(num_elem, 1)) - val(elem_faces(num_elem, 3)));
713
714 dscald[num_elem] = 4. * (deriv_ux * deriv_ux + deriv_vy * deriv_vy) + 2. * ((deriv_uy + deriv_vx) * (deriv_uy + deriv_vx));
715}
716
717// Cas 3D
718void calcul_dscald_interne3D(int num_elem, int elx0, int elx1, int ely0, int ely1, int elz0, int elz1, const Domaine_VDF& domaine_VDF, const DoubleTab& val, DoubleTab& dscald)
719{
720 const IntTab& elem_faces = domaine_VDF.elem_faces();
721 double delta_x_0, delta_x_1, delta_y_0, delta_y_1, delta_z_0, delta_z_1;
722 double delta_x, delta_y, delta_z;
723 double delta_xbis, delta_ybis, delta_zbis;
724 double deriv_ux, deriv_uy, deriv_uz, deriv_vx, deriv_vy, deriv_vz;
725 double deriv_wx, deriv_wy, deriv_wz;
726
727 deriv_ux = 0;
728 deriv_uy = 0;
729 deriv_uz = 0;
730 deriv_vx = 0;
731 deriv_vy = 0;
732 deriv_vz = 0;
733 deriv_wx = 0;
734 deriv_wy = 0;
735 deriv_wz = 0;
736
737 delta_x_0 = domaine_VDF.dist_elem(num_elem, elx0, 0);
738 delta_x_1 = domaine_VDF.dist_elem(elx1, num_elem, 0);
739 delta_y_0 = domaine_VDF.dist_elem(num_elem, ely0, 1);
740 delta_y_1 = domaine_VDF.dist_elem(ely1, num_elem, 1);
741 delta_z_0 = domaine_VDF.dist_elem(num_elem, elz0, 2);
742 delta_z_1 = domaine_VDF.dist_elem(elz1, num_elem, 2);
743
744 delta_x = (delta_x_1 - delta_x_0) * (delta_x_1 + delta_x_0) / (delta_x_1 * delta_x_0);
745 delta_y = (delta_y_1 - delta_y_0) * (delta_y_1 + delta_y_0) / (delta_y_1 * delta_y_0);
746 delta_z = (delta_z_1 - delta_z_0) * (delta_z_1 + delta_z_0) / (delta_z_1 * delta_z_0);
747
748 delta_xbis = domaine_VDF.dim_elem(num_elem, 0);
749 delta_ybis = domaine_VDF.dim_elem(num_elem, 1);
750 delta_zbis = domaine_VDF.dim_elem(num_elem, 2);
751
752 deriv_vz = (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 1)) + delta_z * val(elem_faces(num_elem, 1)) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 1)));
753 deriv_vz += (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 4)) + delta_z * val(elem_faces(num_elem, 4)) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 4)));
754 deriv_vz *= 0.5 / (delta_z_0 + delta_z_1);
755
756 deriv_wy = (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 2)) + delta_y * val(elem_faces(num_elem, 2)) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 2)));
757 deriv_wy += (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 5)) + delta_y * val(elem_faces(num_elem, 5)) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 5)));
758 deriv_wy *= 0.5 / (delta_y_0 + delta_y_1);
759
760 deriv_uz = (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 0)) + delta_z * val(elem_faces(num_elem, 0)) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 0)));
761 deriv_uz += (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 3)) + delta_z * val(elem_faces(num_elem, 3)) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 3)));
762 deriv_uz *= 0.5 / (delta_z_0 + delta_z_1);
763
764 deriv_wx = (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 2)) + delta_x * val(elem_faces(num_elem, 2)) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 2)));
765 deriv_wx += (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 5)) + delta_x * val(elem_faces(num_elem, 5)) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 5)));
766 deriv_wx *= 0.5 / (delta_x_0 + delta_x_1);
767
768 deriv_vx = (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 1)) + delta_x * val(elem_faces(num_elem, 1)) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 1)));
769 deriv_vx += (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 4)) + delta_x * val(elem_faces(num_elem, 4)) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 4)));
770 deriv_vx *= 0.5 / (delta_x_0 + delta_x_1);
771
772 deriv_uy = (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 0)) + delta_y * val(elem_faces(num_elem, 0)) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 0)));
773 deriv_uy += (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 3)) + delta_y * val(elem_faces(num_elem, 3)) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 3)));
774 deriv_uy *= 0.5 / (delta_y_0 + delta_y_1);
775
776 deriv_ux = (1. / delta_xbis) * val(elem_faces(num_elem, 3)) - val(elem_faces(num_elem, 0));
777 deriv_vy = (1. / delta_ybis) * val(elem_faces(num_elem, 4)) - val(elem_faces(num_elem, 1));
778 deriv_wz = (1. / delta_zbis) * val(elem_faces(num_elem, 5)) - val(elem_faces(num_elem, 2));
779
780 dscald[num_elem] = 4. * (deriv_ux * deriv_ux + deriv_vy * deriv_vy + deriv_wz * deriv_wz)
781 + 2. * ((deriv_uy + deriv_vx) * (deriv_uy + deriv_vx) + (deriv_wx + deriv_uz) * (deriv_wx + deriv_uz) + (deriv_wy + deriv_vz) * (deriv_wy + deriv_vz));
782}
783
784// Traitement des elements bord
785
786// Cas 2D
787void calcul_dscald_bord2D(int num_elem, int elx0, int elx1, int ely0, int ely1, const Domaine_VDF& domaine_VDF, const DoubleTab& val, DoubleTab& dscald)
788{
789 const IntTab& elem_faces = domaine_VDF.elem_faces();
790 double delta_x_0, delta_x_1, delta_y_0, delta_y_1;
791 double delta_x, delta_y;
792 double delta_xbis, delta_ybis;
793 double deriv_ux, deriv_vx, deriv_uy, deriv_vy;
794
795 deriv_ux = 0;
796 deriv_vx = 0;
797 deriv_uy = 0;
798 deriv_vy = 0;
799
800 // Traitement des elements bord
801 if ((elx0 == -1) || (elx1 == -1))
802 {
803 if (elx0 == -1)
804 {
805 if (elx1 != -1)
806 {
807 delta_x_1 = domaine_VDF.dist_elem(elx1, num_elem, 0);
808 delta_xbis = domaine_VDF.dim_elem(num_elem, 0);
809
810 deriv_vx = (val(elem_faces(elx1, 1)) - val(elem_faces(num_elem, 1)) + val(elem_faces(elx1, 3)) - val(elem_faces(num_elem, 3)));
811 deriv_vx *= 0.5 / delta_x_1;
812 deriv_ux = (val(elem_faces(num_elem, 2)) - val(elem_faces(num_elem, 0)));
813 deriv_ux *= (1. / delta_xbis);
814
815 }
816 else
817 {
818 deriv_ux = 0;
819 deriv_vx = 0;
820 }
821 }
822 else // elx1 = -1 et elx0 != -1
823 {
824 delta_x_0 = domaine_VDF.dist_elem(num_elem, elx0, 0);
825 delta_xbis = domaine_VDF.dim_elem(num_elem, 0);
826
827 deriv_vx = (val(elem_faces(num_elem, 1)) - val(elem_faces(elx0, 1)) + val(elem_faces(num_elem, 3)) - val(elem_faces(elx0, 3)));
828 deriv_vx *= 0.5 / delta_x_0;
829 deriv_ux = (val(elem_faces(num_elem, 2)) - val(elem_faces(num_elem, 0)));
830 deriv_ux *= (1. / delta_xbis);
831 }
832 }
833 else // elx0 != -1 et elx1 != -1
834 {
835 delta_x_0 = domaine_VDF.dist_elem(num_elem, elx0, 0);
836 delta_x_1 = domaine_VDF.dist_elem(elx1, num_elem, 0);
837 delta_x = (delta_x_1 - delta_x_0) * (delta_x_1 + delta_x_0) / (delta_x_1 * delta_x_0);
838 delta_xbis = domaine_VDF.dim_elem(num_elem, 0);
839
840 deriv_vx = (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 1)) + delta_x * val(elem_faces(num_elem, 1)) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 1)));
841 deriv_vx += (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 3)) + delta_x * val(elem_faces(num_elem, 3)) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 3)));
842 deriv_vx *= 0.5 / (delta_x_0 + delta_x_1);
843 deriv_ux = (val(elem_faces(num_elem, 2)) - val(elem_faces(num_elem, 0)));
844 deriv_ux *= (1. / delta_xbis);
845
846 }
847
848 if ((ely0 == -1) || (ely1 == -1))
849 {
850 if (ely0 == -1)
851 {
852 if (ely1 != -1)
853 {
854 delta_y_1 = domaine_VDF.dist_elem(ely1, num_elem, 1);
855 delta_ybis = domaine_VDF.dim_elem(num_elem, 1);
856
857 deriv_uy = (val(elem_faces(ely1, 0)) - val(elem_faces(num_elem, 0)) + val(elem_faces(ely1, 2)) - val(elem_faces(num_elem, 2)));
858 deriv_uy *= 0.5 / delta_y_1;
859 deriv_vy = (val(elem_faces(num_elem, 3)) - val(elem_faces(num_elem, 1)));
860 deriv_vy *= (1. / delta_ybis);
861
862 }
863 else
864 {
865 deriv_uy = 0;
866 deriv_vy = 0;
867 }
868 }
869 else // ely1 = -1 et ely0 != -1
870 {
871 delta_y_0 = domaine_VDF.dist_elem(num_elem, ely0, 1);
872 delta_ybis = domaine_VDF.dim_elem(num_elem, 1);
873
874 deriv_uy = (val(elem_faces(num_elem, 0)) - val(elem_faces(ely0, 0)) + val(elem_faces(num_elem, 2)) - val(elem_faces(ely0, 2)));
875 deriv_uy *= 0.5 / delta_y_0;
876 deriv_vy = (val(elem_faces(num_elem, 3)) - val(elem_faces(num_elem, 1)));
877 deriv_vy *= (1. / delta_ybis);
878 }
879 }
880 else // ely0 != -1 et ely1 != -1
881 {
882 delta_y_0 = domaine_VDF.dist_elem(num_elem, ely0, 1);
883 delta_y_1 = domaine_VDF.dist_elem(ely1, num_elem, 1);
884 delta_y = (delta_y_1 - delta_y_0) * (delta_y_1 + delta_y_0) / (delta_y_1 * delta_y_0);
885 delta_ybis = domaine_VDF.dim_elem(num_elem, 1);
886
887 deriv_uy = (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 0)) + delta_y * val(elem_faces(num_elem, 0)) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 0)));
888 deriv_uy += (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 2)) + delta_y * val(elem_faces(num_elem, 2)) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 2)));
889 deriv_uy *= 0.5 / (delta_y_0 + delta_y_1);
890 deriv_vy = (val(elem_faces(num_elem, 3)) - val(elem_faces(num_elem, 1)));
891 deriv_vy *= (1. / delta_ybis);
892
893 }
894
895 dscald[num_elem] = 4. * (deriv_ux * deriv_ux + deriv_vy * deriv_vy) + 2. * ((deriv_uy + deriv_vx) * (deriv_uy + deriv_vx));
896}
897
898// Cas 3D
899void calcul_dscald_bord3D(int num_elem, int elx0, int elx1, int ely0, int ely1, int elz0, int elz1, const Domaine_VDF& domaine_VDF, const DoubleTab& val, DoubleTab& dscald)
900{
901 const IntTab& elem_faces = domaine_VDF.elem_faces();
902 double delta_x_0, delta_x_1, delta_y_0, delta_y_1, delta_z_0, delta_z_1;
903 double delta_x, delta_y, delta_z;
904 double delta_xbis, delta_ybis, delta_zbis;
905 double deriv_ux, deriv_uy, deriv_uz, deriv_vx, deriv_vy, deriv_vz;
906 double deriv_wx, deriv_wy, deriv_wz;
907
908 deriv_ux = 0;
909 deriv_uy = 0;
910 deriv_uz = 0;
911 deriv_vx = 0;
912 deriv_vy = 0;
913 deriv_vz = 0;
914 deriv_wx = 0;
915 deriv_wy = 0;
916 deriv_wz = 0;
917
918 // Traitement des elements bord
919
920 if ((elx0 == -1) || (elx1 == -1))
921 {
922 if (elx0 == -1)
923 {
924 if (elx1 != -1)
925 {
926 delta_x_1 = domaine_VDF.dist_elem(elx1, num_elem, 0);
927 delta_xbis = domaine_VDF.dim_elem(num_elem, 0);
928 // A ecrire et a voir sur papier : calcul de deriv_ux pres de la paroi
929 deriv_vx = (val(elem_faces(elx1, 1)) - val(elem_faces(num_elem, 1)) + val(elem_faces(elx1, 4)) - val(elem_faces(num_elem, 4)));
930 deriv_vx *= 0.5 / delta_x_1;
931 deriv_wx = (val(elem_faces(elx1, 2)) - val(elem_faces(num_elem, 2)) + val(elem_faces(elx1, 5)) - val(elem_faces(num_elem, 5)));
932 deriv_wx *= 0.5 / delta_x_1;
933 deriv_ux = (val(elem_faces(num_elem, 3)) - val(elem_faces(num_elem, 0)));
934 deriv_ux *= (1. / delta_xbis);
935 }
936 else
937 {
938 deriv_ux = 0;
939 deriv_vx = 0;
940 deriv_wx = 0;
941 }
942 }
943 else // elx1 = -1 et elx0 != -1
944 {
945 delta_x_0 = domaine_VDF.dist_elem(num_elem, elx0, 0);
946 delta_xbis = domaine_VDF.dim_elem(num_elem, 0);
947
948 deriv_vx = (val(elem_faces(num_elem, 1)) - val(elem_faces(elx0, 1)) + val(elem_faces(num_elem, 4)) - val(elem_faces(elx0, 4)));
949 deriv_vx *= 0.5 / delta_x_0;
950 deriv_wx = (val(elem_faces(num_elem, 2)) - val(elem_faces(elx0, 2)) + val(elem_faces(num_elem, 5)) - val(elem_faces(elx0, 5)));
951 deriv_wx *= 0.5 / delta_x_0;
952 deriv_ux = (val(elem_faces(num_elem, 3)) - val(elem_faces(num_elem, 0)));
953 deriv_ux *= (1. / delta_xbis);
954 }
955
956 if ((ely0 == -1) || (ely1 == -1))
957 {
958 if (ely0 == -1)
959 {
960 if (ely1 != -1)
961 {
962 delta_y_1 = domaine_VDF.dist_elem(ely1, num_elem, 1);
963 delta_ybis = domaine_VDF.dim_elem(num_elem, 1);
964
965 deriv_uy = (val(elem_faces(ely1, 0)) - val(elem_faces(num_elem, 0)) + val(elem_faces(ely1, 3)) - val(elem_faces(num_elem, 3)));
966 deriv_uy *= 0.5 / delta_y_1;
967 deriv_wy = (val(elem_faces(ely1, 2)) - val(elem_faces(num_elem, 2)) + val(elem_faces(ely1, 5)) - val(elem_faces(num_elem, 5)));
968 deriv_wy *= 0.5 / delta_y_1;
969 deriv_vy = (val(elem_faces(num_elem, 4)) - val(elem_faces(num_elem, 1)));
970 deriv_vy *= (1. / delta_ybis);
971 }
972 else
973 {
974 deriv_uy = 0;
975 deriv_vy = 0;
976 deriv_wy = 0;
977 }
978 }
979 else // ely1 = -1 et ely0 != -1
980 {
981 delta_y_0 = domaine_VDF.dist_elem(num_elem, ely0, 1);
982 delta_ybis = domaine_VDF.dim_elem(num_elem, 1);
983
984 deriv_uy = (val(elem_faces(num_elem, 0)) - val(elem_faces(ely0, 0)) + val(elem_faces(num_elem, 3)) - val(elem_faces(ely0, 3)));
985 deriv_uy *= 0.5 / delta_y_0;
986 deriv_wy = (val(elem_faces(num_elem, 2)) - val(elem_faces(ely0, 2)) + val(elem_faces(num_elem, 5)) - val(elem_faces(ely0, 5)));
987 deriv_wy *= 0.5 / delta_y_0;
988 deriv_vy = (val(elem_faces(num_elem, 4)) - val(elem_faces(num_elem, 1)));
989 deriv_vy *= (1. / delta_ybis);
990 }
991 }
992 else
993 {
994 delta_y_0 = domaine_VDF.dist_elem(num_elem, ely0, 1);
995 delta_y_1 = domaine_VDF.dist_elem(ely1, num_elem, 1);
996 delta_y = (delta_y_1 - delta_y_0) * (delta_y_1 + delta_y_0) / (delta_y_1 * delta_y_0);
997 delta_ybis = domaine_VDF.dim_elem(num_elem, 1);
998
999 deriv_uy = (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 0)) + delta_y * val(elem_faces(num_elem, 0)) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 0)));
1000 deriv_uy += (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 3)) + delta_y * val(elem_faces(num_elem, 3)) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 3)));
1001 deriv_uy *= 0.5 / (delta_y_0 + delta_y_1);
1002
1003 deriv_wy = (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 2)) + delta_y * val(elem_faces(num_elem, 2)) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 2)));
1004 deriv_wy += (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 5)) + delta_y * val(elem_faces(num_elem, 5)) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 5)));
1005 deriv_wy *= 0.5 / (delta_y_0 + delta_y_1);
1006 deriv_vy = (val(elem_faces(num_elem, 4)) - val(elem_faces(num_elem, 1)));
1007 deriv_vy *= (1. / delta_ybis);
1008 }
1009
1010 if ((elz0 == -1) || (elz1 == -1))
1011 {
1012 if (elz0 == -1)
1013 {
1014 if (elz1 != -1)
1015 {
1016 delta_z_1 = domaine_VDF.dist_elem(elz1, num_elem, 2);
1017 delta_zbis = domaine_VDF.dim_elem(num_elem, 2);
1018
1019 deriv_uz = (val(elem_faces(elz1, 0)) - val(elem_faces(num_elem, 0)) + val(elem_faces(elz1, 3)) - val(elem_faces(num_elem, 3)));
1020 deriv_uz *= 0.5 / delta_z_1;
1021 deriv_vz = (val(elem_faces(elz1, 1)) - val(elem_faces(num_elem, 1)) + val(elem_faces(elz1, 4)) - val(elem_faces(num_elem, 4)));
1022 deriv_vz *= 0.5 / delta_z_1;
1023 deriv_wz = (val(elem_faces(num_elem, 5)) - val(elem_faces(num_elem, 2)));
1024 deriv_wz *= (1. / delta_zbis);
1025 }
1026 else
1027 {
1028 deriv_uz = 0;
1029 deriv_vz = 0;
1030 deriv_wz = 0;
1031 }
1032 }
1033 else // elz1 = -1 et elz0 != -1
1034 {
1035 delta_z_0 = domaine_VDF.dist_elem(num_elem, elz0, 2);
1036 delta_zbis = domaine_VDF.dim_elem(num_elem, 2);
1037
1038 deriv_uz = (val(elem_faces(num_elem, 0)) - val(elem_faces(elz0, 0)) + val(elem_faces(num_elem, 3)) - val(elem_faces(elz0, 3)));
1039 deriv_uz *= 0.5 / delta_z_0;
1040 deriv_vz = (val(elem_faces(num_elem, 1)) - val(elem_faces(elz0, 1)) + val(elem_faces(num_elem, 4)) - val(elem_faces(elz0, 4)));
1041 deriv_vz *= 0.5 / delta_z_0;
1042 deriv_wz = (val(elem_faces(num_elem, 5)) - val(elem_faces(num_elem, 2)));
1043 deriv_wz *= (1. / delta_zbis);
1044 }
1045 }
1046 else
1047 {
1048 delta_z_0 = domaine_VDF.dist_elem(num_elem, elz0, 2);
1049 delta_z_1 = domaine_VDF.dist_elem(elz1, num_elem, 2);
1050 delta_z = (delta_z_1 - delta_z_0) * (delta_z_1 + delta_z_0) / (delta_z_1 * delta_z_0);
1051 delta_zbis = domaine_VDF.dim_elem(num_elem, 2);
1052
1053 deriv_uz = (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 0)) + delta_z * val(elem_faces(num_elem, 0)) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 0)));
1054 deriv_uz += (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 3)) + delta_z * val(elem_faces(num_elem, 3)) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 3)));
1055 deriv_uz *= 0.5 / (delta_z_0 + delta_z_1);
1056
1057 deriv_vz = (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 1)) + delta_z * val(elem_faces(num_elem, 1)) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 1)));
1058 deriv_vz += (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 4)) + delta_z * val(elem_faces(num_elem, 4)) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 4)));
1059 deriv_vz *= 0.5 / (delta_z_0 + delta_z_1);
1060 deriv_wz = (val(elem_faces(num_elem, 5)) - val(elem_faces(num_elem, 2)));
1061 deriv_wz *= (1. / delta_zbis);
1062 }
1063
1064 }
1065 else if ((ely0 == -1) || (ely1 == -1))
1066 {
1067
1068 if (ely0 == -1)
1069 {
1070 if (ely1 != -1)
1071 {
1072 delta_y_1 = domaine_VDF.dist_elem(ely1, num_elem, 1);
1073 delta_ybis = domaine_VDF.dim_elem(num_elem, 1);
1074
1075 deriv_uy = (val(elem_faces(ely1, 0)) - val(elem_faces(num_elem, 0)) + val(elem_faces(ely1, 3)) - val(elem_faces(num_elem, 3)));
1076 deriv_uy *= 0.5 / delta_y_1;
1077 deriv_wy = (val(elem_faces(ely1, 2)) - val(elem_faces(num_elem, 2)) + val(elem_faces(ely1, 5)) - val(elem_faces(num_elem, 5)));
1078 deriv_wy *= 0.5 / delta_y_1;
1079 deriv_vy = (val(elem_faces(num_elem, 3)) - val(elem_faces(num_elem, 1)));
1080 deriv_vy *= (1. / delta_ybis);
1081 }
1082 else
1083 {
1084 deriv_uy = 0;
1085 deriv_vy = 0;
1086 deriv_wy = 0;
1087 }
1088 }
1089 else // ely1 = -1 et ely0 != -1
1090 {
1091 delta_y_0 = domaine_VDF.dist_elem(num_elem, ely0, 1);
1092 delta_ybis = domaine_VDF.dim_elem(num_elem, 1);
1093
1094 deriv_uy = (val(elem_faces(num_elem, 0)) - val(elem_faces(ely0, 0)) + val(elem_faces(num_elem, 3)) - val(elem_faces(ely0, 3)));
1095 deriv_uy *= 0.5 / delta_y_0;
1096 deriv_wy = (val(elem_faces(num_elem, 2)) - val(elem_faces(ely0, 2)) + val(elem_faces(num_elem, 5)) - val(elem_faces(ely0, 5)));
1097 deriv_wy *= 0.5 / delta_y_0;
1098 deriv_vy = (val(elem_faces(num_elem, 3)) - val(elem_faces(num_elem, 1)));
1099 deriv_vy *= (1. / delta_ybis);
1100 }
1101
1102 delta_x_0 = domaine_VDF.dist_elem(num_elem, elx0, 0);
1103 delta_x_1 = domaine_VDF.dist_elem(elx1, num_elem, 0);
1104 delta_x = (delta_x_1 - delta_x_0) * (delta_x_1 + delta_x_0) / (delta_x_1 * delta_x_0);
1105 delta_xbis = domaine_VDF.dim_elem(num_elem, 0);
1106
1107 deriv_wx = (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 2)) + delta_x * val(elem_faces(num_elem, 2)) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 2)));
1108 deriv_wx += (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 5)) + delta_x * val(elem_faces(num_elem, 5)) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 5)));
1109 deriv_wx *= 0.5 / (delta_x_0 + delta_x_1);
1110
1111 deriv_vx = (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 1)) + delta_x * val(elem_faces(num_elem, 1)) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 1)));
1112 deriv_vx += (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 4)) + delta_x * val(elem_faces(num_elem, 4)) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 4)));
1113 deriv_vx *= 0.5 / (delta_x_0 + delta_x_1);
1114 deriv_ux = (val(elem_faces(num_elem, 3)) - val(elem_faces(num_elem, 0)));
1115 deriv_ux *= (1. / delta_xbis);
1116
1117 if ((elz0 == -1) || (elz1 == -1))
1118 {
1119 if (elz0 == -1)
1120 {
1121 if (elz1 != -1)
1122 {
1123 delta_z_1 = domaine_VDF.dist_elem(elz1, num_elem, 2);
1124 delta_zbis = domaine_VDF.dim_elem(num_elem, 2);
1125
1126 deriv_uz = (val(elem_faces(elz1, 0)) - val(elem_faces(num_elem, 0)) + val(elem_faces(elz1, 3)) - val(elem_faces(num_elem, 3)));
1127 deriv_uz *= 0.5 / delta_z_1;
1128 deriv_vz = (val(elem_faces(elz1, 1)) - val(elem_faces(num_elem, 1)) + val(elem_faces(elz1, 4)) - val(elem_faces(num_elem, 4)));
1129 deriv_vz *= 0.5 / delta_z_1;
1130 deriv_wz = (val(elem_faces(num_elem, 5)) - val(elem_faces(num_elem, 2)));
1131 deriv_wz *= (1. / delta_zbis);
1132 }
1133 else
1134 {
1135 deriv_uz = 0;
1136 deriv_vz = 0;
1137 deriv_wz = 0;
1138 }
1139 }
1140 else
1141 {
1142 delta_z_0 = domaine_VDF.dist_elem(num_elem, elz0, 2);
1143 delta_zbis = domaine_VDF.dim_elem(num_elem, 2);
1144
1145 deriv_uz = (val(elem_faces(num_elem, 0)) - val(elem_faces(elz0, 0)) + val(elem_faces(num_elem, 3)) - val(elem_faces(elz0, 3)));
1146 deriv_uz *= 0.5 / delta_z_0;
1147 deriv_vz = (val(elem_faces(num_elem, 1)) - val(elem_faces(elz0, 1)) + val(elem_faces(num_elem, 4)) - val(elem_faces(elz0, 4)));
1148 deriv_vz *= 0.5 / delta_z_0;
1149 deriv_wz = (val(elem_faces(num_elem, 5)) - val(elem_faces(num_elem, 2)));
1150 deriv_wz *= (1. / delta_zbis);
1151 }
1152 }
1153 else
1154 {
1155 delta_z_0 = domaine_VDF.dist_elem(num_elem, elz0, 2);
1156 delta_z_1 = domaine_VDF.dist_elem(elz1, num_elem, 2);
1157 delta_z = (delta_z_1 - delta_z_0) * (delta_z_1 + delta_z_0) / (delta_z_1 * delta_z_0);
1158 delta_zbis = domaine_VDF.dim_elem(num_elem, 2);
1159
1160 deriv_uz = (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 0)) + delta_z * val(elem_faces(num_elem, 0)) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 0)));
1161 deriv_uz += (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 3)) + delta_z * val(elem_faces(num_elem, 3)) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 3)));
1162 deriv_uz *= 0.5 / (delta_z_0 + delta_z_1);
1163
1164 deriv_vz = (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 1)) + delta_z * val(elem_faces(num_elem, 1)) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 1)));
1165 deriv_vz += (delta_z_0 / delta_z_1 * val(elem_faces(elz1, 4)) + delta_z * val(elem_faces(num_elem, 4)) - delta_z_1 / delta_z_0 * val(elem_faces(elz0, 4)));
1166 deriv_vz *= 0.5 / (delta_z_0 + delta_z_1);
1167 deriv_wz = (val(elem_faces(num_elem, 5)) - val(elem_faces(num_elem, 2)));
1168 deriv_wz *= (1. / delta_zbis);
1169 }
1170 }
1171 else if ((elz0 == -1) || (elz1 == -1))
1172 {
1173 if (elz0 == -1)
1174 {
1175 if (elz1 != -1)
1176 {
1177 delta_z_1 = domaine_VDF.dist_elem(elz1, num_elem, 2);
1178 delta_zbis = domaine_VDF.dim_elem(num_elem, 2);
1179
1180 deriv_uz = (val(elem_faces(elz1, 0)) - val(elem_faces(num_elem, 0)) + val(elem_faces(elz1, 3)) - val(elem_faces(num_elem, 3)));
1181 deriv_uz *= 0.5 / delta_z_1;
1182 deriv_vz = (val(elem_faces(elz1, 1)) - val(elem_faces(num_elem, 1)) + val(elem_faces(elz1, 4)) - val(elem_faces(num_elem, 4)));
1183 deriv_vz *= 0.5 / delta_z_1;
1184 deriv_wz = (val(elem_faces(num_elem, 5)) - val(elem_faces(num_elem, 2)));
1185 deriv_wz *= (1. / delta_zbis);
1186 }
1187 else
1188 {
1189 deriv_uz = 0;
1190 deriv_vz = 0;
1191 deriv_wz = 0;
1192 }
1193 }
1194 else // elz1 = -1 et elz0 != -1
1195 {
1196 delta_z_0 = domaine_VDF.dist_elem(num_elem, elz0, 2);
1197 delta_zbis = domaine_VDF.dim_elem(num_elem, 2);
1198
1199 deriv_uz = (val(elem_faces(num_elem, 0)) - val(elem_faces(elz0, 0)) + val(elem_faces(num_elem, 3)) - val(elem_faces(elz0, 3)));
1200 deriv_uz *= 0.5 / delta_z_0;
1201 deriv_vz = (val(elem_faces(num_elem, 1)) - val(elem_faces(elz0, 1)) + val(elem_faces(num_elem, 4)) - val(elem_faces(elz0, 4)));
1202 deriv_vz *= 0.5 / delta_z_0;
1203 deriv_wz = (val(elem_faces(num_elem, 5)) - val(elem_faces(num_elem, 2)));
1204 deriv_wz *= (1. / delta_zbis);
1205 }
1206
1207 delta_x_0 = domaine_VDF.dist_elem(num_elem, elx0, 0);
1208 delta_x_1 = domaine_VDF.dist_elem(elx1, num_elem, 0);
1209 delta_x = (delta_x_1 - delta_x_0) * (delta_x_1 + delta_x_0) / (delta_x_1 * delta_x_0);
1210 delta_xbis = domaine_VDF.dim_elem(num_elem, 0);
1211
1212 deriv_wx = (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 2)) + delta_x * val(elem_faces(num_elem, 2)) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 2)));
1213 deriv_wx += (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 5)) + delta_x * val(elem_faces(num_elem, 5)) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 5)));
1214 deriv_wx *= 0.5 / (delta_x_0 + delta_x_1);
1215
1216 deriv_vx = (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 1)) + delta_x * val(elem_faces(num_elem, 1)) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 1)));
1217 deriv_vx += (delta_x_0 / delta_x_1 * val(elem_faces(elx1, 4)) + delta_x * val(elem_faces(num_elem, 4)) - delta_x_1 / delta_x_0 * val(elem_faces(elx0, 4)));
1218 deriv_vx *= 0.5 / (delta_x_0 + delta_x_1);
1219 deriv_ux = (val(elem_faces(num_elem, 3)) - val(elem_faces(num_elem, 0)));
1220 deriv_ux *= (1. / delta_xbis);
1221
1222 delta_y_0 = domaine_VDF.dist_elem(num_elem, ely0, 1);
1223 delta_y_1 = domaine_VDF.dist_elem(ely1, num_elem, 1);
1224 delta_y = (delta_y_1 - delta_y_0) * (delta_y_1 + delta_y_0) / (delta_y_1 * delta_y_0);
1225 delta_ybis = domaine_VDF.dim_elem(num_elem, 1);
1226
1227 deriv_wy = (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 2)) + delta_y * val(elem_faces(num_elem, 2)) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 2)));
1228 deriv_wy += (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 5)) + delta_y * val(elem_faces(num_elem, 5)) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 5)));
1229 deriv_wy *= 0.5 / (delta_y_0 + delta_y_1);
1230
1231 deriv_uy = (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 0)) + delta_y * val(elem_faces(num_elem, 0)) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 0)));
1232 deriv_uy += (delta_y_0 / delta_y_1 * val(elem_faces(ely1, 3)) + delta_y * val(elem_faces(num_elem, 3)) - delta_y_1 / delta_y_0 * val(elem_faces(ely0, 3)));
1233 deriv_uy *= 0.5 / (delta_y_0 + delta_y_1);
1234 deriv_vy = (val(elem_faces(num_elem, 4)) - val(elem_faces(num_elem, 1)));
1235 deriv_ux *= (1. / delta_ybis);
1236 }
1237
1238 dscald[num_elem] = 4. * (deriv_ux * deriv_ux + deriv_vy * deriv_vy + deriv_wz * deriv_wz)
1239 + 2. * ((deriv_uy + deriv_vx) * (deriv_uy + deriv_vx) + (deriv_wx + deriv_uz) * (deriv_wx + deriv_uz) + (deriv_wy + deriv_vz) * (deriv_wy + deriv_vz));
1240}
1241
1242// Cas 2D
1243void caldscaldcentelemdim2(DoubleTab& dscald, const DoubleTab& val, const Domaine_VDF& domaine_VDF, int nb_elem, const IntTab& face_voisins, const IntTab& elem_faces)
1244{
1245 if (dscald.dimension(0) != nb_elem)
1246 dscald.resize(nb_elem);
1247 int elx0, elx1, ely0, ely1;
1248
1249 for (int num_elem = 0; num_elem < nb_elem; num_elem++)
1250 {
1251
1252 elx0 = face_voisins(elem_faces(num_elem, 0), 0);
1253 elx1 = face_voisins(elem_faces(num_elem, 2), 1);
1254 ely0 = face_voisins(elem_faces(num_elem, 1), 0);
1255 ely1 = face_voisins(elem_faces(num_elem, 3), 1);
1256
1257 if ((elx0 != -1) && (elx1 != -1) && (ely0 != -1) && (ely1 != -1))
1258 // Cas d'un element interne
1259
1260 calcul_dscald_interne2D(num_elem, elx0, elx1, ely0, ely1, domaine_VDF, val, dscald);
1261 else
1262 calcul_dscald_bord2D(num_elem, elx0, elx1, ely0, ely1, domaine_VDF, val, dscald);
1263 }
1264}
1265
1266// Cas 3D
1267void caldscaldcentelemdim3(DoubleTab& dscald, const DoubleTab& val, const Domaine_VDF& domaine_VDF, int nb_elem, const IntTab& face_voisins, const IntTab& elem_faces)
1268{
1269 if (dscald.dimension(0) != nb_elem)
1270 dscald.resize(nb_elem);
1271 int elx0, elx1, ely0, ely1, elz0, elz1;
1272
1273 for (int num_elem = 0; num_elem < nb_elem; num_elem++)
1274 {
1275
1276 elx0 = face_voisins(elem_faces(num_elem, 0), 0);
1277 elx1 = face_voisins(elem_faces(num_elem, 3), 1);
1278 ely0 = face_voisins(elem_faces(num_elem, 1), 0);
1279 ely1 = face_voisins(elem_faces(num_elem, 4), 1);
1280 elz0 = face_voisins(elem_faces(num_elem, 2), 0);
1281 elz1 = face_voisins(elem_faces(num_elem, 5), 1);
1282
1283 if ((elx0 != -1) && (elx1 != -1) && (ely0 != -1) && (ely1 != -1) && (elz0 != -1) && (elz1 != -1))
1284 // Cas d'un element interne
1285
1286 calcul_dscald_interne3D(num_elem, elx0, elx1, ely0, ely1, elz0, elz1, domaine_VDF, val, dscald);
1287 else
1288 calcul_dscald_bord3D(num_elem, elx0, elx1, ely0, ely1, elz0, elz1, domaine_VDF, val, dscald);
1289 }
1290}
Domaine_VDF
class Domaine_VDF
Definition Domaine_VDF.h:64
Domaine_VDF::dim_elem
double dim_elem(int, int) const
Definition Domaine_VDF.h:627
Domaine_VDF::dist_elem_period
double dist_elem_period(int, int, int) const
Definition Domaine_VDF.h:618
Domaine_VDF::dist_elem
double dist_elem(int, int, int) const
Definition Domaine_VDF.h:610
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
Objet_U::dimension
static int dimension
Definition Objet_U.h:99
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
TRUSTTab::resize
void resize(_SIZE_ n, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTTab.tpp:469
TRUSTTab::dimension
_SIZE_ dimension(int d) const
Definition TRUSTTab.tpp:133
TRUSTVect::line_size
int line_size() const
Definition TRUSTVect.tpp:67