TRUST 1.9.8
HPC thermohydraulic platform
Loading...
Searching...
No Matches
Champ_implementation_P1.h
1/****************************************************************************
2* Copyright (c) 2024, 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 Champ_implementation_P1_included
17#define Champ_implementation_P1_included
18
19#include <Champ_implementation_sommet.h>
20#include <TRUSTTab.h>
21#include <kokkos++.h>
22#include <TRUSTArray_kokkos.tpp>
23
24class Champ_implementation_P1: public Champ_implementation_sommet
25{
26public:
27 ~Champ_implementation_P1() override { }
28 static void init_from_file(DoubleTab& val, const Domaine& dom, int nb_comp, double tolerance, Entree& input);
29
30protected:
31 void value_interpolation(const DoubleTab& positions, const ArrOfInt& cells, const DoubleTab& values, DoubleTab& resu, int ncomp = -1) const override;
32 Champ_base& le_champ() override =0;
33 const Champ_base& le_champ() const override =0;
34
35private:
36 virtual double form_function(const ArrOfDouble& position, const IntTab& les_elems, const DoubleTab& nodes, ArrOfInt& index, int cell, int ddl) const;
37};
38
39extern double coord_barycentrique_P1(const IntTab& polys, const DoubleTab& coord, double x, double y, int le_poly, int i);
40extern double coord_barycentrique_P1(const IntTab& polys, const DoubleTab& coord, double x, double y, double z, int le_poly, int i);
41
42// Defines 4 optimized inlined function called in Champ_implementation_P1.cpp by the 2 previous coord_barycentrique_P1 functions:
43// coord_barycentrique_P1_triangle
44// coord_barycentrique_P1_rectangle
45// coord_barycentrique_P1_tetraedre
46// coord_barycentrique_P1_hexaedre
47KOKKOS_INLINE_FUNCTION double coord_barycentrique_P1_triangle(CIntTabView polys, CDoubleTabView coord, double x, double y, int le_poly, int i)
48{
49 int som0, som1, som2;
50 switch(i)
51 {
52 case 0:
53 {
54 som0 = polys(le_poly, 0);
55 som1 = polys(le_poly, 1);
56 som2 = polys(le_poly, 2);
57 break;
58 }
59 case 1:
60 {
61 som0 = polys(le_poly, 1);
62 som1 = polys(le_poly, 2);
63 som2 = polys(le_poly, 0);
64 break;
65 }
66 case 2:
67 {
68 som0 = polys(le_poly, 2);
69 som1 = polys(le_poly, 0);
70 som2 = polys(le_poly, 1);
71 break;
72 }
73 default:
74 {
75 som0 = -1;
76 som1 = -1;
77 som2 = -1;
78#ifdef KOKKOS
79 Kokkos::abort("A triangle does not have correct nodes in Champ_P1::coord_barycentrique.");
80#else
81 Process::exit("A triangle does not have correct nodes in Champ_P1::coord_barycentrique.");
82#endif
83 }
84 }
85 double den = (coord(som2, 0) - coord(som1, 0)) * (coord(som0, 1) - coord(som1, 1)) - (coord(som2, 1) - coord(som1, 1)) * (coord(som0, 0) - coord(som1, 0));
86
87 double num = (coord(som2, 0) - coord(som1, 0)) * (y - coord(som1, 1)) - (coord(som2, 1) - coord(som1, 1)) * (x - coord(som1, 0));
88
89 double coord_bary = num / den;
90 return coord_bary;
91}
92KOKKOS_INLINE_FUNCTION double coord_barycentrique_P1_tetraedre(CIntTabView polys, CDoubleTabView coord, double x, double y, double z, int le_poly, int i)
93{
94 int som0, som1, som2, som3;
95 switch(i)
96 {
97 case 0:
98 {
99 som0 = polys(le_poly, 0);
100 som1 = polys(le_poly, 1);
101 som2 = polys(le_poly, 2);
102 som3 = polys(le_poly, 3);
103 break;
104 }
105 case 1:
106 {
107 som0 = polys(le_poly, 1);
108 som1 = polys(le_poly, 2);
109 som2 = polys(le_poly, 3);
110 som3 = polys(le_poly, 0);
111 break;
112 }
113 case 2:
114 {
115 som0 = polys(le_poly, 2);
116 som1 = polys(le_poly, 3);
117 som2 = polys(le_poly, 0);
118 som3 = polys(le_poly, 1);
119 break;
120 }
121 case 3:
122 {
123 som0 = polys(le_poly, 3);
124 som1 = polys(le_poly, 0);
125 som2 = polys(le_poly, 1);
126 som3 = polys(le_poly, 2);
127 break;
128 }
129 default:
130 {
131 som0 = -1;
132 som1 = -1;
133 som2 = -1;
134 som3 = -1;
135#ifdef KOKKOS
136 Kokkos::abort("Error in Champ_P1::coord_barycentrique : A tetrahedron does not have correct nodes ");
137#else
138 Process::exit("Error in Champ_P1::coord_barycentrique : A tetrahedron does not have correct nodes ");
139#endif
140 }
141 }
142
143 double xp = (coord(som2, 1) - coord(som1, 1)) * (coord(som0, 2) - coord(som1, 2)) - (coord(som2, 2) - coord(som1, 2)) * (coord(som0, 1) - coord(som1, 1));
144 double yp = (coord(som2, 2) - coord(som1, 2)) * (coord(som0, 0) - coord(som1, 0)) - (coord(som2, 0) - coord(som1, 0)) * (coord(som0, 2) - coord(som1, 2));
145 double zp = (coord(som2, 0) - coord(som1, 0)) * (coord(som0, 1) - coord(som1, 1)) - (coord(som2, 1) - coord(som1, 1)) * (coord(som0, 0) - coord(som1, 0));
146 double den = xp * (coord(som3, 0) - coord(som1, 0)) + yp * (coord(som3, 1) - coord(som1, 1)) + zp * (coord(som3, 2) - coord(som1, 2));
147
148 xp = (coord(som2, 1) - coord(som1, 1)) * (z - coord(som1, 2)) - (coord(som2, 2) - coord(som1, 2)) * (y - coord(som1, 1));
149 yp = (coord(som2, 2) - coord(som1, 2)) * (x - coord(som1, 0)) - (coord(som2, 0) - coord(som1, 0)) * (z - coord(som1, 2));
150 zp = (coord(som2, 0) - coord(som1, 0)) * (y - coord(som1, 1)) - (coord(som2, 1) - coord(som1, 1)) * (x - coord(som1, 0));
151 double num = xp * (coord(som3, 0) - coord(som1, 0)) + yp * (coord(som3, 1) - coord(som1, 1)) + zp * (coord(som3, 2) - coord(som1, 2));
152
153 double coord_bary = num / den;
154 return coord_bary;
155}
156inline double coord_barycentrique_P1_triangle(const IntTab& polys, const DoubleTab& coord, double x, double y, int le_poly, int i)
157{
158 assert(polys.dimension(1) == 3);
159 int som0, som1, som2;
160 switch(i)
161 {
162 case 0:
163 {
164 som0 = polys(le_poly, 0);
165 som1 = polys(le_poly, 1);
166 som2 = polys(le_poly, 2);
167 break;
168 }
169 case 1:
170 {
171 som0 = polys(le_poly, 1);
172 som1 = polys(le_poly, 2);
173 som2 = polys(le_poly, 0);
174 break;
175 }
176 case 2:
177 {
178 som0 = polys(le_poly, 2);
179 som1 = polys(le_poly, 0);
180 som2 = polys(le_poly, 1);
181 break;
182 }
183 default:
184 {
185 som0 = -1;
186 som1 = -1;
187 som2 = -1;
188 Cerr << "A triangle does not have " << i << " nodes in Champ_P1::coord_barycentrique." << finl;
189 Process::exit();
190 }
191 }
192 double den = (coord(som2, 0) - coord(som1, 0)) * (coord(som0, 1) - coord(som1, 1)) - (coord(som2, 1) - coord(som1, 1)) * (coord(som0, 0) - coord(som1, 0));
193
194 double num = (coord(som2, 0) - coord(som1, 0)) * (y - coord(som1, 1)) - (coord(som2, 1) - coord(som1, 1)) * (x - coord(som1, 0));
195
196 assert(den != 0.);
197 double coord_bary = num / den;
198#ifndef NDEBUG
199 if ((coord_bary < -Objet_U::precision_geom) || (coord_bary > 1. + Objet_U::precision_geom))
200 {
201 Cerr << "WARNING: The barycentric coordinate of point :" << finl;
202 Cerr << "x= " << x << " y=" << y << finl;
203 Cerr << "is not between 0 and 1 : " << coord_bary << finl;
204 Cerr << "On the element " << le_poly << " of the processor " << Process::me() << finl;
205 }
206#endif
207 return coord_bary;
208}
209
210inline double coord_barycentrique_P1_rectangle(const IntTab& polys, const DoubleTab& coord, double x, double y, int le_poly, int i)
211{
212 assert(polys.dimension(1) == 4);
213 double alpha_x, alpha_y;
214 switch(i)
215 {
216 case 0:
217 {
218 alpha_x = -1.;
219 alpha_y = -1;
220 break;
221 }
222 case 1:
223 {
224 alpha_x = 1.;
225 alpha_y = -1.;
226 break;
227 }
228 case 2:
229 {
230 alpha_x = -1.;
231 alpha_y = 1.;
232 break;
233 }
234 case 3:
235 {
236 alpha_x = 1.;
237 alpha_y = 1.;
238 break;
239 }
240 default:
241 {
242 alpha_x = 0;
243 alpha_y = 0;
244 Cerr << "Error in Champ_P1::coord_barycentrique : " << finl;
245 Process::exit();
246 }
247 }
248
249 int som0 = polys(le_poly, 0);
250 int som1 = polys(le_poly, 1);
251 int som2 = polys(le_poly, 2);
252 double delta_x = coord(som1, 0) - coord(som0, 0);
253 double delta_y = coord(som2, 1) - coord(som0, 1);
254 double x0 = coord(som0, 0) + delta_x / 2.;
255 double y0 = coord(som0, 1) + delta_y / 2.;
256 double coord_bary = 0.25 * (1. + 2. * alpha_x * (x - x0) / delta_x) * (1. + 2. * alpha_y * (y - y0) / delta_y);
257#ifndef NDEBUG
258 if ((coord_bary < -Objet_U::precision_geom) || (coord_bary > 1. + Objet_U::precision_geom))
259 {
260 Cerr << "WARNING: The barycentric coordinate of point :" << finl;
261 Cerr << "x= " << x << " y=" << y << finl;
262 Cerr << "is not between 0 and 1 : " << coord_bary << finl;
263 Cerr << "On the element " << le_poly << " of the processor " << Process::me() << finl;
264 }
265#endif
266 return coord_bary;
267}
268
269inline double coord_barycentrique_P1_tetraedre(const IntTab& polys, const DoubleTab& coord, double x, double y, double z, int le_poly, int i)
270{
271 int som0, som1, som2, som3;
272 assert(polys.dimension(1) == 4);
273 switch(i)
274 {
275 case 0:
276 {
277 som0 = polys(le_poly, 0);
278 som1 = polys(le_poly, 1);
279 som2 = polys(le_poly, 2);
280 som3 = polys(le_poly, 3);
281 break;
282 }
283 case 1:
284 {
285 som0 = polys(le_poly, 1);
286 som1 = polys(le_poly, 2);
287 som2 = polys(le_poly, 3);
288 som3 = polys(le_poly, 0);
289 break;
290 }
291 case 2:
292 {
293 som0 = polys(le_poly, 2);
294 som1 = polys(le_poly, 3);
295 som2 = polys(le_poly, 0);
296 som3 = polys(le_poly, 1);
297 break;
298 }
299 case 3:
300 {
301 som0 = polys(le_poly, 3);
302 som1 = polys(le_poly, 0);
303 som2 = polys(le_poly, 1);
304 som3 = polys(le_poly, 2);
305 break;
306 }
307 default:
308 {
309 som0 = -1;
310 som1 = -1;
311 som2 = -1;
312 som3 = -1;
313 Cerr << "Error in Champ_P1::coord_barycentrique : " << finl;
314 Cerr << "A tetrahedron does not have " << i << "nodes " << finl;
315 Process::exit();
316 }
317 }
318
319 double xp = (coord(som2, 1) - coord(som1, 1)) * (coord(som0, 2) - coord(som1, 2)) - (coord(som2, 2) - coord(som1, 2)) * (coord(som0, 1) - coord(som1, 1));
320 double yp = (coord(som2, 2) - coord(som1, 2)) * (coord(som0, 0) - coord(som1, 0)) - (coord(som2, 0) - coord(som1, 0)) * (coord(som0, 2) - coord(som1, 2));
321 double zp = (coord(som2, 0) - coord(som1, 0)) * (coord(som0, 1) - coord(som1, 1)) - (coord(som2, 1) - coord(som1, 1)) * (coord(som0, 0) - coord(som1, 0));
322 double den = xp * (coord(som3, 0) - coord(som1, 0)) + yp * (coord(som3, 1) - coord(som1, 1)) + zp * (coord(som3, 2) - coord(som1, 2));
323
324 xp = (coord(som2, 1) - coord(som1, 1)) * (z - coord(som1, 2)) - (coord(som2, 2) - coord(som1, 2)) * (y - coord(som1, 1));
325 yp = (coord(som2, 2) - coord(som1, 2)) * (x - coord(som1, 0)) - (coord(som2, 0) - coord(som1, 0)) * (z - coord(som1, 2));
326 zp = (coord(som2, 0) - coord(som1, 0)) * (y - coord(som1, 1)) - (coord(som2, 1) - coord(som1, 1)) * (x - coord(som1, 0));
327 double num = xp * (coord(som3, 0) - coord(som1, 0)) + yp * (coord(som3, 1) - coord(som1, 1)) + zp * (coord(som3, 2) - coord(som1, 2));
328
329 assert(den != 0.);
330 double coord_bary = num / den;
331#ifndef NDEBUG
332 if ((coord_bary < -Objet_U::precision_geom) || (coord_bary > 1. + Objet_U::precision_geom))
333 {
334 Cerr << "WARNING: The barycentric coordinate of point :" << finl;
335 Cerr << "x= " << x << " y=" << y << " z=" << z << finl;
336 Cerr << "is not between 0 and 1 : " << coord_bary << finl;
337 Cerr << "On the element " << le_poly << " of the processor " << Process::me() << finl;
338 }
339#endif
340 return coord_bary;
341}
342inline double coord_barycentrique_P1_hexaedre(const IntTab& polys, const DoubleTab& coord, double x, double y, double z, int le_poly, int i)
343{
344 assert(polys.dimension(1) == 8);
345 double alpha_x, alpha_y, alpha_z;
346 switch(i)
347 {
348 case 0:
349 {
350 alpha_x = -1.;
351 alpha_y = -1.;
352 alpha_z = -1.;
353 break;
354 }
355 case 1:
356 {
357 alpha_x = 1.;
358 alpha_y = -1.;
359 alpha_z = -1.;
360 break;
361 }
362 case 2:
363 {
364 alpha_x = -1.;
365 alpha_y = 1.;
366 alpha_z = -1.;
367 break;
368 }
369 case 3:
370 {
371 alpha_x = 1.;
372 alpha_y = 1.;
373 alpha_z = -1.;
374 break;
375 }
376 case 4:
377 {
378 alpha_x = -1.;
379 alpha_y = -1;
380 alpha_z = 1.;
381 break;
382 }
383 case 5:
384 {
385 alpha_x = 1.;
386 alpha_y = -1.;
387 alpha_z = 1.;
388 break;
389 }
390 case 6:
391 {
392 alpha_x = -1.;
393 alpha_y = 1.;
394 alpha_z = 1.;
395 break;
396 }
397 case 7:
398 {
399 alpha_x = 1.;
400 alpha_y = 1.;
401 alpha_z = 1.;
402 break;
403 }
404 default:
405 {
406 alpha_x = 0;
407 alpha_y = 0;
408 alpha_z = 0;
409 Cerr << "Error in Champ_P1::coord_barycentrique : " << finl;
410 Process::exit();
411 }
412 }
413
414 int som0 = polys(le_poly, 0);
415 int som1 = polys(le_poly, 1);
416 int som2 = polys(le_poly, 2);
417 int som3 = polys(le_poly, 4);
418
419 double delta_x = coord(som1, 0) - coord(som0, 0);
420 double delta_y = coord(som2, 1) - coord(som0, 1);
421 double delta_z = coord(som3, 2) - coord(som0, 2);
422
423 double x0 = coord(som0, 0) + delta_x / 2.;
424 double y0 = coord(som0, 1) + delta_y / 2.;
425 double z0 = coord(som0, 2) + delta_z / 2.;
426
427 double coord_bary = (1. / 8.) * (1. + 2. * alpha_x * (x - x0) / delta_x) * (1. + 2. * alpha_y * (y - y0) / delta_y) * (1. + 2. * alpha_z * (z - z0) / delta_z);
428#ifndef NDEBUG
429 if ((coord_bary < -Objet_U::precision_geom) || (coord_bary > 1. + Objet_U::precision_geom))
430 {
431 Cerr << "WARNING: The barycentric coordinate of point :" << finl;
432 Cerr << "x= " << x << " y=" << y << " z=" << z << finl;
433 Cerr << "is not between 0 and 1 : " << coord_bary << finl;
434 Cerr << "On the element " << le_poly << " of the processor " << Process::me() << finl;
435 }
436#endif
437 return coord_bary;
438}
439
440#endif /* Champ_implementation_P1_inclus */
Champ_base
classe Champ_base Cette classe est la base de la hierarchie des champs.
Definition Champ_base.h:43
Champ_implementation_P1
Definition Champ_implementation_P1.h:25
Champ_implementation_P1::le_champ
const Champ_base & le_champ() const override=0
Champ_implementation_P1::le_champ
Champ_base & le_champ() override=0
Champ_implementation_P1::init_from_file
static void init_from_file(DoubleTab &val, const Domaine &dom, int nb_comp, double tolerance, Entree &input)
Initialise le tableau de valeurs aux sommets du domaine dom a partir de valeurs lues dans "input".
Definition Champ_implementation_P1.cpp:177
Champ_implementation_P1::~Champ_implementation_P1
~Champ_implementation_P1() override
Definition Champ_implementation_P1.h:27
Champ_implementation_P1::value_interpolation
void value_interpolation(const DoubleTab &positions, const ArrOfInt &cells, const DoubleTab &values, DoubleTab &resu, int ncomp=-1) const override
Definition Champ_implementation_P1.cpp:98
Champ_implementation_sommet
: class Champ_implementation_sommet
Definition Champ_implementation_sommet.h:29
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Objet_U::precision_geom
static double precision_geom
Definition Objet_U.h:86
Process::me
static int me()
renvoie mon rang dans le groupe de communication courant.
Definition Process.cpp:125
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
TRUSTTab::dimension
_SIZE_ dimension(int d) const
Definition TRUSTTab.tpp:133