TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
distances_VEF.cpp
1/****************************************************************************
2* Copyright (c) 2025, 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#include <distances_VEF.h>
16#include <Domaine.h>
17#include <Motcle.h>
18
19// See norm_vit1 in distances_VEF.h for Kokkos functions
20double norm_2D_vit1(const DoubleTab& vit,int num1,int num2,int fac,const Domaine_VEF& domaine,double& u)
21{
22 const DoubleTab& face_normale = domaine.face_normales();
23 const Domaine& domaine_geom = domaine.domaine();
24 int nfac = domaine_geom.nb_faces_elem();
25
26 // fac numero de la face a paroi fixe
27 double r0,r1;
28 calcule_r0r1(face_normale,fac,r0,r1);
29
30 double v1=(vit(num1,0)+vit(num2,0))/nfac;
31 double v2=(vit(num1,1)+vit(num2,1))/nfac;
32 double v = vitesse_tangentielle(v1,v2,r0,r1);
33 if (v==0)
34 u = 0;
35 else
36 u = 1;
37 return v;
38}
39
40// See norm_vit1 in distances_VEF.h for Kokkos functions
41double norm_2D_vit1(const DoubleTab& vit,int num1,int num2,int fac,const Domaine_VEF& domaine,double& val1, double& val2)
42{
43 const DoubleTab& face_normale = domaine.face_normales();
44 const Domaine& domaine_geom = domaine.domaine();
45 int nfac = domaine_geom.nb_faces_elem();
46
47 // fac numero de la face a paroi fixe
48 double r0,r1;
49 calcule_r0r1(face_normale,fac,r0,r1);
50
51 double v1=(vit(num1,0)+vit(num2,0))/nfac;
52 double v2=(vit(num1,1)+vit(num2,1))/nfac;
53 double v = vitesse_tangentielle(v1,v2,r0,r1);
54 double psc = r0*v1+r1*v2;
55
56 // val1 et val2 sont les vitesses tangentielles
57 val1=(v1-psc*r0)/(v+DMINFLOAT);
58 val2=(v2-psc*r1)/(v+DMINFLOAT);
59
60 return v;
61}
62
63
64double norm_2D_vit1_lp(const DoubleTab& vit,int fac,int num1,int num2,const Domaine_VEF& domaine,double& val1, double& val2)
65{
66 // PQ : 03/03 : Definition de la vitesse tangentielle a une paroi
67 // obtenue par "moyennage" des vitesses aux faces associees a des elements standards
68 // (dont on suppose le bon comportement suivant la loi log)
69
70 const DoubleTab& face_normale = domaine.face_normales();
71 double c1,c2;
72 c1=0.;
73 c2=0.;
74
75 // int rang1 = rang_elem_non_std(num1);
76 // if (rang1<0) c1=1.;
77 // int rang2 = rang_elem_non_std(num2);
78 // if (rang2<0) c2=1.;
79
80 c1=1.;
81 c2=1.;
82
83 double r0,r1;
84 calcule_r0r1(face_normale,fac,r0,r1);
85
86 double v0=(c1*vit(num1,0)+c2*vit(num2,0))/(c1+c2);
87 double v1=(c1*vit(num1,1)+c2*vit(num2,1))/(c1+c2);
88
89 double psc = r0 * v0 + r1 * v1;
90 double norm_vit = vitesse_tangentielle(v0,v1,r0,r1);
91 // composantes du vecteur tangent normalise
92 val1 = (v0-psc*r0)/(norm_vit+DMINFLOAT);
93 val2 = (v1-psc*r1)/(norm_vit+DMINFLOAT);
94
95 return norm_vit;
96}
97
98// See norm_vit1 in distances_VEF.h for Kokkos functions
99double norm_2D_vit1(const DoubleTab& vit,int num1,int num2,int num3,int fac,const Domaine_VEF& domaine,double& u)
100{
101 const DoubleTab& face_normale = domaine.face_normales();
102 const Domaine& domaine_geom = domaine.domaine();
103 int nfac = domaine_geom.nb_faces_elem();
104
105 // fac numero de la face a paroi fixe
106 double r0,r1;
107 calcule_r0r1(face_normale,fac,r0,r1);
108 double v1=(vit(num1,0)+vit(num2,0)+vit(num3,0))/nfac;
109 double v2=(vit(num1,1)+vit(num2,1)+vit(num3,1))/nfac;
110 double v = vitesse_tangentielle(v1,v2,r0,r1);
111 if (v==0)
112 u = 0;
113 else
114 u = 1;
115 return v;
116}
117
118// See norm_vit1 in distances_VEF.h for Kokkos functions
119double norm_2D_vit1(const DoubleTab& vit,int num1,int num2,int num3,int num4,int fac,const Domaine_VEF& domaine,double& u)
120{
121 const DoubleTab& face_normale = domaine.face_normales();
122 const Domaine& domaine_geom = domaine.domaine();
123 int nfac = domaine_geom.nb_faces_elem();
124
125 // fac numero de la face a paroi fixe
126 double r0,r1;
127 calcule_r0r1(face_normale,fac,r0,r1);
128 double v1=(vit(num1,0)+vit(num2,0)+vit(num3,0)+vit(num4,0))/nfac;
129 double v2=(vit(num1,1)+vit(num2,1)+vit(num3,1)+vit(num4,1))/nfac;
130 double v = vitesse_tangentielle(v1,v2,r0,r1);
131 if (v==0)
132 u = 0;
133 else
134 u = 1;
135 return v;
136}
137
138
139double norm_2D_vit2(const DoubleTab& vit,int num1,int num2,int fac,const Domaine_VEF& domaine, double& u)
140{
141 const DoubleTab& face_normale = domaine.face_normales();
142 const Domaine& domaine_geom = domaine.domaine();
143 int nfac = domaine_geom.nb_faces_elem();
144 // fac numero de la face a paroi defilante
145 double r0,r1;
146 calcule_r0r1(face_normale,fac,r0,r1);
147 double v1=(vit(num1,0)+vit(num2,0)-2.*vit(fac,0))/nfac;
148 double v2=(vit(num1,1)+vit(num2,1)-2.*vit(fac,1))/nfac;
149 double v = vitesse_tangentielle(v1,v2,r0,r1);
150 if (v==0)
151 u = 0;
152 else
153 u = 1;
154 return v;
155}
156
157double norm_2D_vit2(const DoubleTab& vit,int num1,int num2,int num3,int fac,const Domaine_VEF& domaine, double& u)
158{
159 const DoubleTab& face_normale = domaine.face_normales();
160 const Domaine& domaine_geom = domaine.domaine();
161 int nfac = domaine_geom.nb_faces_elem();
162 // fac numero de la face a paroi defilante
163 double r0,r1;
164 calcule_r0r1(face_normale,fac,r0,r1);
165 double v1=(vit(num1,0)+vit(num2,0)+vit(num3,0)-3.*vit(fac,0))/nfac;
166 double v2=(vit(num1,1)+vit(num2,1)+vit(num3,1)-3.*vit(fac,1))/nfac;
167 double v = vitesse_tangentielle(v1,v2,r0,r1);
168 if (v==0)
169 u = 0;
170 else
171 u = 1;
172 return v;
173}
174
175double norm_2D_vit2(const DoubleTab& vit,int num1,int num2,int num3,int num4,int fac,const Domaine_VEF& domaine,double& u)
176{
177 const DoubleTab& face_normale = domaine.face_normales();
178 const Domaine& domaine_geom = domaine.domaine();
179 int nfac = domaine_geom.nb_faces_elem();
180
181 // fac numero de la face a paroi fixe
182 double r0,r1;
183 calcule_r0r1(face_normale,fac,r0,r1);
184 double v1=(vit(num1,0)+vit(num2,0)+vit(num3,0)+vit(num4,0)-4.*vit(fac,0))/nfac;
185 double v2=(vit(num1,1)+vit(num2,1)+vit(num3,1)+vit(num4,1)-4.*vit(fac,1))/nfac;
186 double v = vitesse_tangentielle(v1,v2,r0,r1);
187 if (v==0)
188 u = 0;
189 else
190 u = 1;
191 return v;
192}
193
194// See distance in distances_VEF.h for Kokkos functions
195double distance_2D(int fac,int elem,const Domaine_VEF& domaine)
196{
197 const DoubleTab& xp = domaine.xp(); // centre de gravite des elements
198 const DoubleTab& xv = domaine.xv(); // centre de gravite des faces
199
200 const DoubleTab& face_normale = domaine.face_normales();
201 double r0,r1;
202 calcule_r0r1(face_normale,fac,r0,r1);
203
204 double x0=xv(fac,0);
205 double y0=xv(fac,1);
206 double x1=xp(elem,0);
207 double y1=xp(elem,1);
208
209 return std::fabs(r0*(x1-x0)+r1*(y1-y0));
210}
211
212
213double norm_2D_vit1_k(const DoubleTab& vit,int fac,int num1,
214 const Domaine_VEF& domaine,
215 double& val1,double& val2)
216{
217 const DoubleTab& face_normale = domaine.face_normales();
218 // fac numero de la face a paroi fixe
219 double r0,r1;
220 calcule_r0r1(face_normale,fac,r0,r1);
221 double v1=vit(num1,0);
222 double v2=vit(num1,1);
223 double norm_vit = vitesse_tangentielle(v1,v2,r0,r1);
224
225 double psc = r0*v1+r1*v2;
226
227 // val1 et val2 sont les vitesses tangetentielles
228 val1=(v1-psc*r0)/(norm_vit+DMINFLOAT);
229 val2=(v2-psc*r1)/(norm_vit+DMINFLOAT);
230
231
232 return norm_vit;
233}
234
235double norm_3D_vit1_k(const DoubleTab& vit,int fac,int num1,
236 const Domaine_VEF& domaine,
237 double& val1,double& val2,double& val3)
238{
239 const DoubleTab& face_normale = domaine.face_normales();
240 // fac numero de la face a paroi fixe
241 double r0,r1,r2;
242 calcule_r0r1r2(face_normale,fac,r0,r1,r2);
243 double v1=vit(num1,0);
244 double v2=vit(num1,1);
245 double v3=vit(num1,2);
246 double norm_vit = vitesse_tangentielle(v1,v2,v3,r0,r1,r2);
247
248 double psc = r0*v1+r1*v2+r2*v3;
249 // val1,val2 val3 sont les vitesses tangentielles
250 val1=(v1-psc*r0)/(norm_vit+DMINFLOAT);
251 val2=(v2-psc*r1)/(norm_vit+DMINFLOAT);
252 val3=(v3-psc*r2)/(norm_vit+DMINFLOAT);
253
254 return norm_vit;
255}
256
257// See norm_vit1 in distances_VEF.h for Kokkos functions
258double norm_3D_vit1(const DoubleTab& vit,int fac,int num1,int num2,int num3,
259 const Domaine_VEF& domaine,
260 double& val1,double& val2,double& val3)
261{
262 const DoubleTab& face_normale = domaine.face_normales();
263 const Domaine& domaine_geom = domaine.domaine();
264 int nfac = domaine_geom.nb_faces_elem();
265 // fac numero de la face a paroi fixe
266 double r0,r1,r2;
267 calcule_r0r1r2(face_normale,fac,r0,r1,r2);
268 double v1=(vit(num1,0)+vit(num2,0)+vit(num3,0))/nfac;
269 double v2=(vit(num1,1)+vit(num2,1)+vit(num3,1))/nfac;
270 double v3=(vit(num1,2)+vit(num2,2)+vit(num3,2))/nfac;
271 double norm_vit = vitesse_tangentielle(v1,v2,v3,r0,r1,r2);
272 // val1,val2 val3 sont les vitesses tangentielles
273 double psc = r0*v1+r1*v2+r2*v3;
274 val1=(v1-psc*r0)/(norm_vit+DMINFLOAT);
275 val2=(v2-psc*r1)/(norm_vit+DMINFLOAT);
276 val3=(v3-psc*r2)/(norm_vit+DMINFLOAT);
277
278 return norm_vit;
279}
280
281double norm_3D_vit1(const DoubleTab& vit,int fac,int num1,int num2,int num3,int num4,
282 const Domaine_VEF& domaine,
283 double& val1,double& val2,double& val3)
284{
285 const DoubleTab& face_normale = domaine.face_normales();
286 const Domaine& domaine_geom = domaine.domaine();
287 int nfac = domaine_geom.nb_faces_elem();
288 // fac numero de la face a paroi fixe
289 double r0,r1,r2;
290 calcule_r0r1r2(face_normale,fac,r0,r1,r2);
291 double v1=(vit(num1,0)+vit(num2,0)+vit(num3,0)+vit(num4,0))/nfac;
292 double v2=(vit(num1,1)+vit(num2,1)+vit(num3,1)+vit(num4,1))/nfac;
293 double v3=(vit(num1,2)+vit(num2,2)+vit(num3,2)+vit(num4,2))/nfac;
294 double norm_vit = vitesse_tangentielle(v1,v2,v3,r0,r1,r2);
295 // val1,val2 val3 sont les vitesses tangentielles
296 double psc = r0*v1+r1*v2+r2*v3;
297 val1=(v1-psc*r0)/(norm_vit+DMINFLOAT);
298 val2=(v2-psc*r1)/(norm_vit+DMINFLOAT);
299 val3=(v3-psc*r2)/(norm_vit+DMINFLOAT);
300 return norm_vit;
301}
302double norm_3D_vit1(const DoubleTab& vit,int fac,int num1,int num2,int num3,int num4,
303 int num5,const Domaine_VEF& domaine,
304 double& val1,double& val2,double& val3)
305{
306 const DoubleTab& face_normale = domaine.face_normales();
307 const Domaine& domaine_geom = domaine.domaine();
308 int nfac = domaine_geom.nb_faces_elem();
309 // fac numero de la face a paroi fixe
310 double r0,r1,r2;
311 calcule_r0r1r2(face_normale,fac,r0,r1,r2);
312 double v1=(vit(num1,0)+vit(num2,0)+vit(num3,0)+vit(num4,0)+vit(num5,0))/nfac;
313 double v2=(vit(num1,1)+vit(num2,1)+vit(num3,1)+vit(num4,1)+vit(num5,1))/nfac;
314 double v3=(vit(num1,2)+vit(num2,2)+vit(num3,2)+vit(num4,2)+vit(num5,2))/nfac;
315 double norm_vit = vitesse_tangentielle(v1,v2,v3,r0,r1,r2);
316 // val1,val2 val3 sont les vitesses tangentielles
317 double psc = r0*v1+r1*v2+r2*v3;
318 val1=(v1-psc*r0)/(norm_vit+DMINFLOAT);
319 val2=(v2-psc*r1)/(norm_vit+DMINFLOAT);
320 val3=(v3-psc*r2)/(norm_vit+DMINFLOAT);
321
322 return norm_vit;
323}
324
325double norm_2D_vit2_k(const DoubleTab& vit,int fac,int num1,
326 const Domaine_VEF& domaine,
327 double& val1,double& val2)
328{
329 const DoubleTab& face_normale = domaine.face_normales();
330 // fac numero de la face a paroi defilante
331 double r0,r1;
332 calcule_r0r1(face_normale,fac,r0,r1);
333
334 double v1=vit(num1,0)-vit(fac,0);
335 double v2=vit(num1,1)-vit(fac,1);
336 double norm_vit = vitesse_tangentielle(v1,v2,r0,r1);
337 // val1 et val2 sont les vitesses tangetentielles
338 double psc = r0*v1+r1*v2;
339 val1=(v1-psc*r0)/(norm_vit+DMINFLOAT);
340 val2=(v2-psc*r1)/(norm_vit+DMINFLOAT);
341
342 return norm_vit;
343}
344double norm_3D_vit2_k(const DoubleTab& vit,int fac,int num1,
345 const Domaine_VEF& domaine,
346 double& val1,double& val2,double& val3)
347{
348 const DoubleTab& face_normale = domaine.face_normales();
349 // fac numero de la face a paroi defilante
350 double r0,r1,r2;
351 calcule_r0r1r2(face_normale,fac,r0,r1,r2);
352 double v1=vit(num1,0)-vit(fac,0);
353 double v2=vit(num1,1)-vit(fac,1);
354 double v3=vit(num1,2)-vit(fac,2);
355 double norm_vit = vitesse_tangentielle(v1,v2,v3,r0,r1,r2);
356 // val1,val2 val3 sont les vitesses tangentielles
357 double psc = r0*v1+r1*v2+r2*v3;
358 val1=(v1-psc*r0)/(norm_vit+DMINFLOAT);
359 val2=(v2-psc*r1)/(norm_vit+DMINFLOAT);
360 val3=(v3-psc*r2)/(norm_vit+DMINFLOAT);
361
362 return norm_vit;
363}
364
365double norm_3D_vit2(const DoubleTab& vit,int fac,int num1,int num2,int num3,
366 const Domaine_VEF& domaine,
367 double& val1,double& val2,double& val3)
368{
369 const DoubleTab& face_normale = domaine.face_normales();
370 const Domaine& domaine_geom = domaine.domaine();
371 int nfac = domaine_geom.nb_faces_elem();
372 // fac numero de la face a paroi defilante
373 double r0,r1,r2;
374 calcule_r0r1r2(face_normale,fac,r0,r1,r2);
375 double v1=(vit(num1,0)+vit(num2,0)+vit(num3,0)-(nfac-1)*vit(fac,0))/nfac;
376 double v2=(vit(num1,1)+vit(num2,1)+vit(num3,1)-(nfac-1)*vit(fac,1))/nfac;
377 double v3=(vit(num1,2)+vit(num2,2)+vit(num3,2)-(nfac-1)*vit(fac,2))/nfac;
378 double norm_vit = vitesse_tangentielle(v1,v2,v3,r0,r1,r2);
379 // val1,val2 val3 sont les vitesses tangentielles
380 double psc = r0*v1+r1*v2+r2*v3;
381 val1=(v1-psc*r0)/(norm_vit+DMINFLOAT);
382 val2=(v2-psc*r1)/(norm_vit+DMINFLOAT);
383 val3=(v3-psc*r2)/(norm_vit+DMINFLOAT);
384 return norm_vit;
385}
386
387double norm_3D_vit1_lp(const DoubleTab& vit,int fac,int num1,int num2,int num3,
388 const Domaine_VEF& domaine,
389 double& val1,double& val2,double& val3)
390{
391 // PQ : 03/03 : Definition de la vitesse tangentielle a une paroi
392 // obtenue par "moyennage" des vitesses aux faces associees a des elements standards
393 // (dont on suppose le bon comportement suivant la loi log)
394
395 const DoubleTab& face_normale = domaine.face_normales();
396 // const IntVect& rang_elem_non_std = domaine.rang_elem_non_std();
397 double c1,c2,c3;
398 c1=0.;
399 c2=0.;
400 c3=0.;
401
402 // int rang1 = rang_elem_non_std(num1);
403 // if (rang1<0) c1=1.;
404 // int rang2 = rang_elem_non_std(num2);
405 // if (rang2<0) c2=1.;
406 // int rang3 = rang_elem_non_std(num3);
407 // if (rang3<0) c3=1.;
408
409 c1=1.;
410 c2=1.;
411 c3=1.;
412
413 double r0,r1,r2;
414 calcule_r0r1r2(face_normale,fac,r0,r1,r2);
415 double v0=(c1*vit(num1,0)+c2*vit(num2,0)+c3*vit(num3,0))/(c1+c2+c3);
416 double v1=(c1*vit(num1,1)+c2*vit(num2,1)+c3*vit(num3,1))/(c1+c2+c3);
417 double v2=(c1*vit(num1,2)+c2*vit(num2,2)+c3*vit(num3,2))/(c1+c2+c3);
418
419 double psc = r0 * v0 + r1 * v1 + r2 * v2;
420 double norm_vit = vitesse_tangentielle(v0,v1,v2,r0,r1,r2);
421 // composantes du vecteur tangent normalise
422 val1 = (v0-psc*r0)/(norm_vit+DMINFLOAT);
423 val2 = (v1-psc*r1)/(norm_vit+DMINFLOAT);
424 val3 = (v2-psc*r2)/(norm_vit+DMINFLOAT);
425
426 return norm_vit;
427}
428
429double norm_3D_vit2(const DoubleTab& vit,int fac,int num1,int num2,int num3,int num4,
430 const Domaine_VEF& domaine,
431 double& val1,double& val2,double& val3)
432{
433 const DoubleTab& face_normale = domaine.face_normales();
434 const Domaine& domaine_geom = domaine.domaine();
435 int nfac = domaine_geom.nb_faces_elem();
436 // fac numero de la face a paroi defilante
437 double r0,r1,r2;
438 calcule_r0r1r2(face_normale,fac,r0,r1,r2);
439
440 double v1=(vit(num1,0)+vit(num2,0)+vit(num3,0)+vit(num4,0)-nfac*vit(fac,0))/nfac;
441 double v2=(vit(num1,1)+vit(num2,1)+vit(num3,1)+vit(num4,1)-nfac*vit(fac,1))/nfac;
442 double v3=(vit(num1,2)+vit(num2,2)+vit(num3,2)+vit(num4,2)-nfac*vit(fac,2))/nfac;
443 double norm_vit = vitesse_tangentielle(v1,v2,v3,r0,r1,r2);
444 // val1,val2 val3 sont les vitesses tangentielles
445 double psc = r0*v1+r1*v2+r2*v3;
446 val1=(v1-psc*r0)/(norm_vit+DMINFLOAT);
447 val2=(v2-psc*r1)/(norm_vit+DMINFLOAT);
448 val3=(v3-psc*r2)/(norm_vit+DMINFLOAT);
449
450 return norm_vit;
451}
452double norm_3D_vit2(const DoubleTab& vit,int fac,int num1,int num2,int num3,int num4,
453 int num5, const Domaine_VEF& domaine,
454 double& val1,double& val2,double& val3)
455{
456 const DoubleTab& face_normale = domaine.face_normales();
457 const Domaine& domaine_geom = domaine.domaine();
458 int nfac = domaine_geom.nb_faces_elem();
459 // fac numero de la face a paroi defilante
460 double r0,r1,r2;
461 calcule_r0r1r2(face_normale,fac,r0,r1,r2);
462
463 double v1=(vit(num1,0)+vit(num2,0)+vit(num3,0)+vit(num4,0)+vit(num5,0)-nfac*vit(fac,0))/nfac;
464 double v2=(vit(num1,1)+vit(num2,1)+vit(num3,1)+vit(num4,1)+vit(num5,1)-nfac*vit(fac,1))/nfac;
465 double v3=(vit(num1,2)+vit(num2,2)+vit(num3,2)+vit(num4,2)+vit(num5,2)-nfac*vit(fac,2))/nfac;
466 double norm_vit = vitesse_tangentielle(v1,v2,v3,r0,r1,r2);
467 // val1,val2 val3 sont les vitesses tangentielles
468 double psc = r0*v1+r1*v2+r2*v3;
469 val1=(v1-psc*r0)/(norm_vit+DMINFLOAT);
470 val2=(v2-psc*r1)/(norm_vit+DMINFLOAT);
471 val3=(v3-psc*r2)/(norm_vit+DMINFLOAT);
472
473 return norm_vit;
474}
475
476// See distance in distances_VEF.h for Kokkos functions
477double distance_3D(int fac,int elem,const Domaine_VEF& domaine)
478{
479 const DoubleTab& xp = domaine.xp(); // centre de gravite des elements
480 const DoubleTab& xv = domaine.xv(); // centre de gravite des faces
481 const DoubleTab& face_normale = domaine.face_normales();
482 double r0,r1,r2;
483 calcule_r0r1r2(face_normale,fac,r0,r1,r2);
484 double x0=xv(fac,0);
485 double y0=xv(fac,1);
486 double z0=xv(fac,2);
487 double x1=xp(elem,0);
488 double y1=xp(elem,1);
489 double z1=xp(elem,2);
490
491 return std::fabs(r0*(x1-x0)+r1*(y1-y0)+r2*(z1-z0));
492}
493
494
495
496DoubleVect& calcul_longueur_filtre(DoubleVect& longueur_filtre, const Motcle& methode, const Domaine_VEF& domaine)
497{
498 int nbr_element=domaine.nb_elem_tot();
499 int element;
500 int dim=Objet_U::dimension;
501 const Domaine& domaine_geom = domaine.domaine();
502
503 if (longueur_filtre.size() != nbr_element)
504 {
505 Cerr << "erreur dans la taille du DoubleVect valeurs de la longueur du filtre" << finl;
506 Process::exit();
507 }
508
509 if (methode == Motcle("volume") || methode == Motcle("volume_sans_lissage")) // racine cubique du volume
510 {
511 longueur_filtre=-1.;
512
513 const DoubleVect& volume = domaine.volumes();
514 for (element=0; element<nbr_element; element ++)
515 {
516 longueur_filtre(element) = exp(log(volume[element])/double(dim));
517 }
518 }
519 else if (methode == Motcle("arete")) // recherche de la plus longue arete d'un element
520 {
521 longueur_filtre=-1.;
522
523 const IntTab& les_sommets = domaine_geom.les_elems();
524 int som1,som2,som_1,som_2;
525 double distance;
526
527 for (element=0; element<nbr_element; element ++)
528 for (som1=0; som1<dim; som1++)
529 for (som2=som1; som2<dim+1; som2++)
530 {
531 som_1 = les_sommets(element, som1);
532 som_2 = les_sommets(element, som2);
533
534 distance = distance_sommets(som_1, som_2, domaine);
535 distance/=sqrt(2.);
536 longueur_filtre(element) = std::max(longueur_filtre(element), distance);
537 }
538 }
539 else if (methode == Motcle("Scotti")) // application de Scotti a un pseudo-cube
540 {
541 longueur_filtre=-1.;
542
543 const IntTab& les_sommets = domaine_geom.les_elems();
544 int som0,som1,som2,som3;
545 int som_0,som_1,som_2,som_3;
546
547 if(Objet_U::dimension==2) // On revient a la racine carre du volume
548 {
549 int nbr_elementb=domaine.nb_elem_tot();
550 int elementb;
551 //int dim=Objet_U::dimension;
552 const DoubleVect& volume = domaine.volumes();
553
554 for (elementb=0; elementb<nbr_elementb; elementb ++)
555 {
556 longueur_filtre(elementb) = exp(log(volume[elementb])/double(dim));
557 }
558 }
559 else
560 {
561 ArrOfDouble psc(4);
562 ArrOfDouble dist(3);
563 double dist_tot,dist_min,dist_max,dist_moy;
564 double a1,a2,f_scotti;
565
566 som_0=0;
567 som_1=0;
568 som_2=0;
569 som_3=0;
570
571 for (element=0; element<nbr_element; element ++)
572 {
573 som0 = les_sommets(element, 0);
574 som1 = les_sommets(element, 1);
575 som2 = les_sommets(element, 2);
576 som3 = les_sommets(element, 3);
577
578 psc[0] = som_pscal(som0,som1,som2,som3,domaine);
579 psc[1] = som_pscal(som1,som0,som2,som3,domaine);
580 psc[2] = som_pscal(som2,som0,som1,som3,domaine);
581 psc[3] = som_pscal(som3,som0,som1,som2,domaine);
582
583 const int indice_min = imin_array(psc);
584 if(indice_min==0)
585 {
586 som_0=som0;
587 som_1=som1;
588 som_2=som2;
589 som_3=som3;
590 }
591 if(indice_min==1)
592 {
593 som_0=som1;
594 som_1=som0;
595 som_2=som2;
596 som_3=som3;
597 }
598 if(indice_min==2)
599 {
600 som_0=som2;
601 som_1=som0;
602 som_2=som1;
603 som_3=som3;
604 }
605 if(indice_min==3)
606 {
607 som_0=som3;
608 som_1=som0;
609 som_2=som1;
610 som_3=som2;
611 }
612
613 dist[0]=distance_sommets(som_0,som_1,domaine);
614 dist[1]=distance_sommets(som_0,som_2,domaine);
615 dist[2]=distance_sommets(som_0,som_3,domaine);
616
617 dist_tot=dist[0]+dist[1]+dist[2];
618
619 dist_min=min_array(dist);
620 dist_max=max_array(dist);
621 dist_moy=dist_tot-dist_min-dist_max;
622
623 a1=dist_min/dist_max;
624 a2=dist_moy/dist_max;
625
626 f_scotti=cosh(sqrt((4./27.)*( log(a1)*log(a1) - log(a1)*log(a2) + log(a2)*log(a2) )));
627
628 longueur_filtre(element) = f_scotti * exp(log(dist_min*dist_max*dist_moy)/3.);
629 }
630 }//3D
631 }
632 else
633 {
634 Cerr << "calcul_longueur_filtre.cpp n'a pas reconnu l'argument : " << methode << finl;
635 Cerr << "les arguments possibles sont : \"volume\", \"volume_sans_lissage\", \"Scotti\", \"arete\"." << finl;
636 Process::exit();
637
638 }
639
640
641 if ( ! (methode == Motcle("volume_sans_lissage")) ) // processus de "regularisation"
642 {
643 const Domaine& dom=domaine.domaine();
644 const IntTab& les_sommets = domaine_geom.les_elems();
645 int nb_sommet = domaine.nb_som_tot();
646 ArrOfDouble longueur_filtre_sommet(nb_sommet);
647 int som1;
648 int som_0,som_1;
649
650 longueur_filtre_sommet=-1.;
651
652 for (element=0; element<nbr_element; element ++)
653 for (som1=0; som1<dim+1; som1++)
654 {
655 som_0 = les_sommets(element, som1);
656 som_1 = dom.get_renum_som_perio(som_0);
657 longueur_filtre_sommet[som_1] = std::max(longueur_filtre(element), longueur_filtre_sommet[som_1]);
658 }
659
660 longueur_filtre=-1.;
661
662 for (element=0; element<nbr_element; element ++)
663 for (som1=0; som1<dim+1; som1++)
664 {
665 som_0 = les_sommets(element, som1);
666 som_1 = dom.get_renum_som_perio(som_0);
667 longueur_filtre(element) = std::max (longueur_filtre(element), longueur_filtre_sommet[som_1]);
668 }
669 }
670
671 assert(nbr_element == 0 || min_array(longueur_filtre)>0.);
672
673 return longueur_filtre;
674}
675
676double distance_sommets(const int sommet1, const int sommet2, const Domaine_VEF& domaine)
677{
678 const Domaine& domaine_geom = domaine.domaine();
679 const DoubleTab& xs = domaine_geom.coord_sommets();
680 double x1 = xs(sommet1,0);
681 double y1 = xs(sommet1,1);
682 double x2 = xs(sommet2,0);
683 double y2 = xs(sommet2,1);
684
685 if (Objet_U::dimension == 3)
686 {
687 double z1 = xs(sommet1,2);
688 double z2 = xs(sommet2,2);
689 return sqrt( carre(x2-x1) + carre(y2-y1) + carre(z2-z1) );
690 }
691 else // en dimension 2
692 {
693 return sqrt( carre(x2-x1) + carre(y2-y1) );
694 }
695}
696
697double som_pscal(const int som0, const int som1, const int som2, const int som3, const Domaine_VEF& domaine)
698{
699 const Domaine& domaine_geom = domaine.domaine();
700 const DoubleTab& xs = domaine_geom.coord_sommets();
701
702 ArrOfDouble v1(3),v2(3),v3(3);
703 double n1,n2,n3;
704
705 for(int i=0 ; i<Objet_U::dimension ; i++)
706 {
707 v1[i]=xs(som1,i)-xs(som0,i);
708 v2[i]=xs(som2,i)-xs(som0,i);
709 v3[i]=xs(som3,i)-xs(som0,i);
710 }
711
712 n1=norme_array(v1);
713 n2=norme_array(v2);
714 n3=norme_array(v3);
715
716 v1/=n1;
717 v2/=n2;
718 v3/=n3;
719
720 double som_psc=std::fabs(dotproduct_array(v1,v2))+std::fabs(dotproduct_array(v2,v3))+std::fabs(dotproduct_array(v3,v1));
721 return som_psc;
722}
723
724double norm_vit_lp_k(const DoubleTab& vit,int face,int face_b,const Domaine_VEF& domaine,ArrOfDouble& val,int is_defilante)
725{
726 int dimension=Objet_U::dimension;
727 //assert(face_b<domaine.premiere_face_int()); assert a ameliorer ou faces_virt...
728 if (is_defilante==0)
729 if (dimension==3)
730 return norm_3D_vit1_k(vit,face_b,face,domaine,val[0],val[1],val[2]);
731 else
732 return norm_2D_vit1_k(vit,face_b,face,domaine,val[0],val[1]);
733
734 else
735 {
736 if (dimension==3)
737 return norm_3D_vit2_k(vit,face_b,face,domaine,val[0],val[1],val[2]);
738 else
739 return norm_2D_vit2_k(vit,face_b,face,domaine,val[0],val[1]);
740 }
741}
742
743double distance_face_elem(int fac,int elem,const Domaine_VEF& domaine)
744{
745 int dimension=Objet_U::dimension;
746 if (dimension==3)
747 return distance_3D(fac,elem,domaine);
748 else
749 return distance_2D(fac,elem,domaine);
750}
751
Domaine_32_64::get_renum_som_perio
int_t get_renum_som_perio(int_t i) const
Definition Domaine.h:281
Domaine_32_64::les_elems
IntTab_t & les_elems()
Definition Domaine.h:129
Domaine_32_64::nb_faces_elem
int nb_faces_elem(int=0) const
Renvoie le nombre de face de type i des elements geometriques constituants le domaine.
Definition Domaine.h:484
Domaine_32_64::coord_sommets
const DoubleTab_t & coord_sommets() const
Definition Domaine.h:112
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:54
Motcle
Une chaine de caractere (Nom) en majuscules.
Definition Motcle.h:26
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
TRUSTVect::size
_SIZE_ size() const
Definition TRUSTVect.tpp:45