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> & function,
44 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
45 double const & time)
46 {
47 (void)function;
48 (void)q_points;
49 (void)time;
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(ContainerInterfaceData<rank, dim, double> const & function,
59 unsigned int const face,
60 unsigned int const q,
61 unsigned int const quad_index)
62 {
63 (void)function;
64 (void)face;
65 (void)q;
66 (void)quad_index;
67
68 AssertThrow(false, dealii::ExcMessage("should not arrive here."));
69
70 return dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>>();
71 }
72
73 static inline DEAL_II_ALWAYS_INLINE //
74 dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>>
75 value(FunctionWithNormal<dim> & function_with_normal,
76 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
77 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & normals,
78 double const & time)
79 {
80 (void)function_with_normal;
81 (void)q_points;
82 (void)normals;
83 (void)time;
84
85 AssertThrow(false, dealii::ExcMessage("not implemented."));
86
87 return dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>>();
88 }
89
90 static inline DEAL_II_ALWAYS_INLINE //
91 dealii::SymmetricTensor<rank, dim, dealii::VectorizedArray<Number>>
92 value_symmetric(dealii::Function<dim> & function,
93 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
94 double const & time)
95 {
96 (void)function;
97 (void)q_points;
98 (void)time;
99
100 AssertThrow(false, dealii::ExcMessage("should not arrive here."));
101
102 return dealii::SymmetricTensor<rank, dim, dealii::VectorizedArray<Number>>();
103 }
104
105 static inline DEAL_II_ALWAYS_INLINE //
106 dealii::SymmetricTensor<rank, dim, dealii::VectorizedArray<Number>>
107 value_symmetric(ContainerInterfaceData<rank, dim, double> const & function,
108 unsigned int const face,
109 unsigned int const q,
110 unsigned int const quad_index)
111 {
112 (void)function;
113 (void)face;
114 (void)q;
115 (void)quad_index;
116
117 AssertThrow(false, dealii::ExcMessage("should not arrive here."));
118
119 return dealii::SymmetricTensor<rank, dim, dealii::VectorizedArray<Number>>();
120 }
121};
122
123template<int dim, typename Number>
124struct FunctionEvaluator<0, dim, Number>
125{
126 static inline DEAL_II_ALWAYS_INLINE //
127 dealii::Tensor<0, dim, dealii::VectorizedArray<Number>>
128 value(dealii::Function<dim> & function,
129 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
130 double const & time)
131 {
132 dealii::VectorizedArray<Number> value = dealii::make_vectorized_array<Number>(0.0);
133
134 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
135 {
136 dealii::Point<dim> q_point;
137 for(unsigned int d = 0; d < dim; ++d)
138 q_point[d] = q_points[d][v];
139
140 function.set_time(time);
141 value[v] = function.value(q_point);
142 }
143
144 return value;
145 }
146
147 static inline DEAL_II_ALWAYS_INLINE //
148 dealii::Tensor<0, dim, dealii::VectorizedArray<Number>>
149 value(ContainerInterfaceData<0, dim, double> const & function,
150 unsigned int const face,
151 unsigned int const q,
152 unsigned int const quad_index)
153 {
154 dealii::VectorizedArray<Number> value = dealii::make_vectorized_array<Number>(0.0);
155
156 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
157 value[v] = function.get_data(quad_index, face, q, v);
158
159 return value;
160 }
161};
162
163template<int dim, typename Number>
164struct FunctionEvaluator<1, dim, Number>
165{
166 static inline DEAL_II_ALWAYS_INLINE //
167 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
168 value(dealii::Function<dim> & function,
169 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
170 double const & time)
171 {
172 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> value;
173
174 for(unsigned int d = 0; d < dim; ++d)
175 {
176 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
177 {
178 dealii::Point<dim> q_point;
179 for(unsigned int i = 0; i < dim; ++i)
180 q_point[i] = q_points[i][v];
181
182 function.set_time(time);
183 value[d][v] = function.value(q_point, d);
184 }
185 }
186
187 return value;
188 }
189
190 static inline DEAL_II_ALWAYS_INLINE //
191 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
192 value(ContainerInterfaceData<1, dim, double> const & function,
193 unsigned int const face,
194 unsigned int const q,
195 unsigned int const quad_index)
196 {
197 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> value;
198
199 std::array<dealii::Tensor<1, dim, Number>, dealii::VectorizedArray<Number>::size()>
200 tensor_array;
201
202 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
203 tensor_array[v] = function.get_data(quad_index, face, q, v);
204
205 for(unsigned int d = 0; d < dim; ++d)
206 {
207 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
208 value[d][v] = tensor_array[v][d];
209 }
210
211 return value;
212 }
213
214 static inline DEAL_II_ALWAYS_INLINE //
215 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
216 value(FunctionWithNormal<dim> & function_with_normal,
217 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
218 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & normals,
219 double const & time)
220 {
221 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> value;
222
223 for(unsigned int d = 0; d < dim; ++d)
224 {
225 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
226 {
227 dealii::Point<dim> q_point;
228 dealii::Tensor<1, dim> normal;
229 for(unsigned int i = 0; i < dim; ++i)
230 {
231 q_point[i] = q_points[i][v];
232 normal[i] = normals[i][v];
233 }
234 function_with_normal.set_time(time);
235 function_with_normal.set_normal_vector(normal);
236 value[d][v] = function_with_normal.value(q_point, d);
237 }
238 }
239
240 return value;
241 }
242};
243
244template<int dim, typename Number>
245struct FunctionEvaluator<2, dim, Number>
246{
247 static inline DEAL_II_ALWAYS_INLINE //
248 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>>
249 value(dealii::Function<dim> & function,
250 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
251 double const & time)
252 {
253 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> value;
254
255 for(unsigned int d1 = 0; d1 < dim; ++d1)
256 {
257 for(unsigned int d2 = 0; d2 < dim; ++d2)
258 {
259 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
260 {
261 dealii::Point<dim> q_point;
262
263 for(unsigned int i = 0; i < dim; ++i)
264 q_point[i] = q_points[i][v];
265
266 function.set_time(time);
267
268 auto const unrolled_index =
269 dealii::Tensor<2, dim>::component_to_unrolled_index(dealii::TableIndices<2>(d1, d2));
270
271 value[d1][d2][v] = function.value(q_point, unrolled_index);
272 }
273 }
274 }
275
276 return value;
277 }
278
279 static inline DEAL_II_ALWAYS_INLINE //
280 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>>
281 value(ContainerInterfaceData<2, dim, double> const & function,
282 unsigned int const face,
283 unsigned int const q,
284 unsigned int const quad_index)
285 {
286 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> value;
287
288 std::array<dealii::Tensor<2, dim, Number>, dealii::VectorizedArray<Number>::size()>
289 tensor_array;
290
291 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
292 tensor_array[v] = function.get_data(quad_index, face, q, v);
293
294 for(unsigned int d1 = 0; d1 < dim; ++d1)
295 {
296 for(unsigned int d2 = 0; d2 < dim; ++d2)
297 {
298 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
299 value[d1][d2][v] = tensor_array[v][d1][d2];
300 }
301 }
302
303 return value;
304 }
305
306 static inline DEAL_II_ALWAYS_INLINE //
307 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>>
308 value_symmetric(dealii::Function<dim> & function,
309 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_points,
310 double const & time)
311 {
312 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>> value;
313
314 for(unsigned int d1 = 0; d1 < dim; ++d1)
315 {
316 for(unsigned int d2 = d1; d2 < dim; ++d2)
317 {
318 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
319 {
320 dealii::Point<dim> q_point;
321
322 for(unsigned int i = 0; i < dim; ++i)
323 q_point[i] = q_points[i][v];
324
325 function.set_time(time);
326
327 auto const unrolled_index = dealii::SymmetricTensor<2, dim>::component_to_unrolled_index(
328 dealii::TableIndices<2>(d1, d2));
329
330 value[d1][d2][v] = function.value(q_point, unrolled_index);
331 }
332 }
333 }
334
335 return value;
336 }
337
338 static inline DEAL_II_ALWAYS_INLINE //
339 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>>
340 value_symmetric(ContainerInterfaceData<2, dim, double> const & function,
341 unsigned int const face,
342 unsigned int const q,
343 unsigned int const quad_index)
344 {
345 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>> value;
346
347 std::array<dealii::SymmetricTensor<2, dim, Number>, dealii::VectorizedArray<Number>::size()>
348 tensor_array;
349
350 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
351 tensor_array[v] = function.get_data(quad_index, face, q, v);
352
353 for(unsigned int d1 = 0; d1 < dim; ++d1)
354 {
355 for(unsigned int d2 = d1; d2 < dim; ++d2)
356 {
357 for(unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
358 value[d1][d2][v] = tensor_array[v][d1][d2];
359 }
360 }
361
362 return value;
363 }
364};
365
366} // namespace ExaDG
367
368#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.10.0