ExaDG
Loading...
Searching...
No Matches
evaluate_functions.h
1/* ______________________________________________________________________
2 *
3 * ExaDG - High-Order Discontinuous Galerkin for the Exa-Scale
4 *
5 * Copyright (C) 2021 by the ExaDG authors
6 *
7 * This program is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program. If not, see <https://www.gnu.org/licenses/>.
19 * ______________________________________________________________________
20 */
21
22#ifndef INCLUDE_EXADG_FUNCTIONS_AND_BOUNDARY_CONDITIONS_EVALUATE_FUNCTIONS_H_
23#define INCLUDE_EXADG_FUNCTIONS_AND_BOUNDARY_CONDITIONS_EVALUATE_FUNCTIONS_H_
24
25// deal.II
26#include <deal.II/base/function.h>
27#include <deal.II/base/point.h>
28#include <deal.II/base/tensor.h>
29#include <deal.II/base/vectorization.h>
30
31#include <exadg/functions_and_boundary_conditions/container_interface_data.h>
32#include <exadg/functions_and_boundary_conditions/function_with_normal.h>
33
34#include <memory>
35
36namespace ExaDG
37{
38template<int rank, int dim, typename Number>
39struct FunctionEvaluator
40{
41 static inline DEAL_II_ALWAYS_INLINE //
42 dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>>
43 value(dealii::Function<dim> const & function,
44 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points)
45 {
46 (void)function;
47 (void)q_points;
48
49 AssertThrow(false, dealii::ExcMessage("should not arrive here."));
50
51 return dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>>();
52 }
53
54 static inline DEAL_II_ALWAYS_INLINE //
55 dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>>
56 value(dealii::Function<dim> & function,
57 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
58 double const & time)
59 {
60 function.set_time(time);
61 return value(function, q_points);
62 }
63
64 static inline DEAL_II_ALWAYS_INLINE //
65 dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>>
66 value(ContainerInterfaceData<rank, dim, double> const & function,
67 unsigned int const face,
68 unsigned int const q,
69 unsigned int const quad_index)
70 {
71 (void)function;
72 (void)face;
73 (void)q;
74 (void)quad_index;
75
76 AssertThrow(false, dealii::ExcMessage("should not arrive here."));
77
78 return dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>>();
79 }
80
81 static inline DEAL_II_ALWAYS_INLINE //
82 dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>>
83 value(FunctionWithNormal<dim> const & function_with_normal,
84 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
85 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & normals)
86 {
87 (void)function_with_normal;
88 (void)q_points;
89 (void)normals;
90
91 AssertThrow(false, dealii::ExcMessage("not implemented."));
92
93 return dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>>();
94 }
95
96 static inline DEAL_II_ALWAYS_INLINE //
97 dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>>
98 value(FunctionWithNormal<dim> & function_with_normal,
99 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
100 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & normals,
101 double const & time)
102 {
103 function_with_normal.set_time(time);
104 return value(function_with_normal, q_points, normals);
105 }
106
107 static inline DEAL_II_ALWAYS_INLINE //
108 dealii::SymmetricTensor<rank, dim, dealii::VectorizedArray<Number>>
109 value_symmetric(dealii::Function<dim> const & function,
110 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points)
111 {
112 (void)function;
113 (void)q_points;
114
115 AssertThrow(false, dealii::ExcMessage("should not arrive here."));
116
117 return dealii::SymmetricTensor<rank, dim, dealii::VectorizedArray<Number>>();
118 }
119
120 static inline DEAL_II_ALWAYS_INLINE //
121 dealii::SymmetricTensor<rank, dim, dealii::VectorizedArray<Number>>
122 value_symmetric(dealii::Function<dim> & function,
123 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
124 double const & time)
125 {
126 function.set_time(time);
127 return value_symmetric(function, q_points);
128 }
129
130 static inline DEAL_II_ALWAYS_INLINE //
131 dealii::SymmetricTensor<rank, dim, dealii::VectorizedArray<Number>>
132 value_symmetric(ContainerInterfaceData<rank, dim, double> const & function,
133 unsigned int const face,
134 unsigned int const q,
135 unsigned int const quad_index)
136 {
137 (void)function;
138 (void)face;
139 (void)q;
140 (void)quad_index;
141
142 AssertThrow(false, dealii::ExcMessage("should not arrive here."));
143
144 return dealii::SymmetricTensor<rank, dim, dealii::VectorizedArray<Number>>();
145 }
146};
147
148template<int dim, typename Number>
149struct FunctionEvaluator<0, dim, Number>
150{
151 static inline DEAL_II_ALWAYS_INLINE //
152 dealii::Tensor<0, dim, dealii::VectorizedArray<Number>>
153 value(dealii::Function<dim> const & function,
154 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points)
155 {
156 dealii::VectorizedArray<Number> value = dealii::make_vectorized_array<Number>(0.0);
157
158 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
159 {
160 dealii::Point<dim> q_point;
161 for(unsigned int d = 0; d < dim; ++d)
162 q_point[d] = q_points[d][v];
163
164 value[v] = function.value(q_point);
165 }
166
167 return value;
168 }
169
170 static inline DEAL_II_ALWAYS_INLINE //
171 dealii::Tensor<0, dim, dealii::VectorizedArray<Number>>
172 value(dealii::Function<dim> & function,
173 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
174 double const & time)
175 {
176 function.set_time(time);
177 return value(function, q_points);
178 }
179
180
181 static inline DEAL_II_ALWAYS_INLINE //
182 dealii::Tensor<0, dim, dealii::VectorizedArray<Number>>
183 value(ContainerInterfaceData<0, dim, double> const & function,
184 unsigned int const face,
185 unsigned int const q,
186 unsigned int const quad_index)
187 {
188 dealii::VectorizedArray<Number> value = dealii::make_vectorized_array<Number>(0.0);
189
190 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
191 value[v] = function.get_data(quad_index, face, q, v);
192
193 return value;
194 }
195};
196
197template<int dim, typename Number>
198struct FunctionEvaluator<1, dim, Number>
199{
200 static inline DEAL_II_ALWAYS_INLINE //
201 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
202 value(dealii::Function<dim> const & function,
203 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points)
204 {
205 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> value;
206
207 for(unsigned int d = 0; d < dim; ++d)
208 {
209 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
210 {
211 dealii::Point<dim> q_point;
212 for(unsigned int i = 0; i < dim; ++i)
213 q_point[i] = q_points[i][v];
214
215 value[d][v] = function.value(q_point, d);
216 }
217 }
218
219 return value;
220 }
221
222 static inline DEAL_II_ALWAYS_INLINE //
223 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
224 value(dealii::Function<dim> & function,
225 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
226 double const & time)
227 {
228 function.set_time(time);
229 return value(function, q_points);
230 }
231
232 static inline DEAL_II_ALWAYS_INLINE //
233 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
234 value(ContainerInterfaceData<1, dim, double> const & function,
235 unsigned int const face,
236 unsigned int const q,
237 unsigned int const quad_index)
238 {
239 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> value;
240
241 std::array<dealii::Tensor<1, dim, Number>, dealii::VectorizedArray<Number>::size()>
242 tensor_array;
243
244 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
245 tensor_array[v] = function.get_data(quad_index, face, q, v);
246
247 for(unsigned int d = 0; d < dim; ++d)
248 {
249 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
250 value[d][v] = tensor_array[v][d];
251 }
252
253 return value;
254 }
255
256 static inline DEAL_II_ALWAYS_INLINE //
257 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
258 value(FunctionWithNormal<dim> & function_with_normal,
259 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
260 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & normals)
261 {
262 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> value;
263
264 for(unsigned int d = 0; d < dim; ++d)
265 {
266 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
267 {
268 dealii::Point<dim> q_point;
269 dealii::Tensor<1, dim> normal;
270 for(unsigned int i = 0; i < dim; ++i)
271 {
272 q_point[i] = q_points[i][v];
273 normal[i] = normals[i][v];
274 }
275 function_with_normal.set_normal_vector(normal);
276 value[d][v] = function_with_normal.value(q_point, d);
277 }
278 }
279
280 return value;
281 }
282
283 static inline DEAL_II_ALWAYS_INLINE //
284 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
285 value(FunctionWithNormal<dim> & function_with_normal,
286 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
287 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & normals,
288 double const & time)
289 {
290 function_with_normal.set_time(time);
291 return value(function_with_normal, q_points, normals);
292 }
293};
294
295template<int dim, typename Number>
296struct FunctionEvaluator<2, dim, Number>
297{
298 static inline DEAL_II_ALWAYS_INLINE //
299 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>>
300 value(dealii::Function<dim> const & function,
301 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points)
302 {
303 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> value;
304
305 for(unsigned int d1 = 0; d1 < dim; ++d1)
306 {
307 for(unsigned int d2 = 0; d2 < dim; ++d2)
308 {
309 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
310 {
311 dealii::Point<dim> q_point;
312
313 for(unsigned int i = 0; i < dim; ++i)
314 q_point[i] = q_points[i][v];
315
316 auto const unrolled_index =
317 dealii::Tensor<2, dim>::component_to_unrolled_index(dealii::TableIndices<2>(d1, d2));
318
319 value[d1][d2][v] = function.value(q_point, unrolled_index);
320 }
321 }
322 }
323
324 return value;
325 }
326
327 static inline DEAL_II_ALWAYS_INLINE //
328 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>>
329 value(dealii::Function<dim> & function,
330 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
331 double const & time)
332 {
333 function.set_time(time);
334 return value(function, q_points);
335 }
336
337 static inline DEAL_II_ALWAYS_INLINE //
338 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>>
339 value(ContainerInterfaceData<2, dim, double> const & function,
340 unsigned int const face,
341 unsigned int const q,
342 unsigned int const quad_index)
343 {
344 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> value;
345
346 std::array<dealii::Tensor<2, dim, Number>, dealii::VectorizedArray<Number>::size()>
347 tensor_array;
348
349 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
350 tensor_array[v] = function.get_data(quad_index, face, q, v);
351
352 for(unsigned int d1 = 0; d1 < dim; ++d1)
353 {
354 for(unsigned int d2 = 0; d2 < dim; ++d2)
355 {
356 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
357 value[d1][d2][v] = tensor_array[v][d1][d2];
358 }
359 }
360
361 return value;
362 }
363
364 static inline DEAL_II_ALWAYS_INLINE //
365 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>>
366 value_symmetric(dealii::Function<dim> const & function,
367 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points)
368 {
369 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>> value;
370
371 for(unsigned int d1 = 0; d1 < dim; ++d1)
372 {
373 for(unsigned int d2 = d1; d2 < dim; ++d2)
374 {
375 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
376 {
377 dealii::Point<dim> q_point;
378
379 for(unsigned int i = 0; i < dim; ++i)
380 q_point[i] = q_points[i][v];
381
382 auto const unrolled_index = dealii::SymmetricTensor<2, dim>::component_to_unrolled_index(
383 dealii::TableIndices<2>(d1, d2));
384
385 value[d1][d2][v] = function.value(q_point, unrolled_index);
386 }
387 }
388 }
389
390 return value;
391 }
392
393 static inline DEAL_II_ALWAYS_INLINE //
394 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>>
395 value_symmetric(dealii::Function<dim> & function,
396 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
397 double const & time)
398 {
399 function.set_time(time);
400 return value(q_points, time);
401 }
402
403 static inline DEAL_II_ALWAYS_INLINE //
404 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>>
405 value_symmetric(ContainerInterfaceData<2, dim, double> const & function,
406 unsigned int const face,
407 unsigned int const q,
408 unsigned int const quad_index)
409 {
410 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>> value;
411
412 std::array<dealii::SymmetricTensor<2, dim, Number>, dealii::VectorizedArray<Number>::size()>
413 tensor_array;
414
415 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
416 tensor_array[v] = function.get_data(quad_index, face, q, v);
417
418 for(unsigned int d1 = 0; d1 < dim; ++d1)
419 {
420 for(unsigned int d2 = d1; d2 < dim; ++d2)
421 {
422 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
423 value[d1][d2][v] = tensor_array[v][d1][d2];
424 }
425 }
426
427 return value;
428 }
429};
430
431} // namespace ExaDG
432
433#endif /* INCLUDE_EXADG_FUNCTIONS_AND_BOUNDARY_CONDITIONS_EVALUATE_FUNCTIONS_H_ */
ExaDG::ContainerInterfaceData
Definition container_interface_data.h:45
ExaDG::FunctionWithNormal
Definition function_with_normal.h:33
ExaDG
Definition driver.cpp:33
ExaDG::FunctionEvaluator
Definition evaluate_functions.h:40
Generated by doxygen 1.13.2