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