Serac  0.1
Serac is an implicit thermal strucural mechanics simulation code.
tetrahedron_L2.inl
Go to the documentation of this file.
1 // Copyright (c) Lawrence Livermore National Security, LLC and
2 // other Serac Project Developers. See the top-level LICENSE file for
3 // details.
4 //
5 // SPDX-License-Identifier: (BSD-3-Clause)
6 
13 // this specialization defines shape functions (and their gradients) that
14 // interpolate at Gauss-Lobatto nodes for the appropriate polynomial order
15 //
16 // note: mfem assumes the parent element domain is the convex hull of {{0,0,0}, {1,0,0}, {0,1,0}, {0,0,1}}
17 // for additional information on the finite_element concept requirements, see finite_element.hpp
18 //
19 // for exact positions of nodes for different polynomial orders, see simplex_basis_function_unit_tests.cpp
21 template <int p, int c>
22 struct finite_element<mfem::Geometry::TETRAHEDRON, L2<p, c> > {
23  static constexpr auto geometry = mfem::Geometry::TETRAHEDRON;
24  static constexpr auto family = Family::L2;
25  static constexpr int components = c;
26  static constexpr int dim = 3;
27  static constexpr int n = (p + 1);
28  static constexpr int ndof = (p + 1) * (p + 2) * (p + 3) / 6;
29  static constexpr int nqpts(int q) { return num_quadrature_points(mfem::Geometry::TETRAHEDRON, q); }
30 
31  static constexpr int VALUE = 0, GRADIENT = 1;
32  static constexpr int SOURCE = 0, FLUX = 1;
33 
34  using residual_type =
35  typename std::conditional<components == 1, tensor<double, ndof>, tensor<double, ndof, components> >::type;
36 
37  using dof_type = tensor<double, c, ndof>;
38 
39  using value_type = typename std::conditional<components == 1, double, tensor<double, components> >::type;
40  using derivative_type =
41  typename std::conditional<components == 1, tensor<double, dim>, tensor<double, components, dim> >::type;
42  using qf_input_type = tuple<value_type, derivative_type>;
43 
44  SERAC_HOST_DEVICE static constexpr double shape_function([[maybe_unused]] tensor<double, dim> xi,
45  [[maybe_unused]] int i)
46  {
47  if constexpr (p == 0) {
48  return 1.0;
49  }
50 
51  if constexpr (p == 1) {
52  switch (i) {
53  case 0:
54  return 1 - xi[0] - xi[1] - xi[2];
55  case 1:
56  return xi[0];
57  case 2:
58  return xi[1];
59  case 3:
60  return xi[2];
61  }
62  }
63  if constexpr (p == 2) {
64  switch (i) {
65  case 0:
66  return (-1 + xi[0] + xi[1] + xi[2]) * (-1 + 2 * xi[0] + 2 * xi[1] + 2 * xi[2]);
67  case 1:
68  return -4 * xi[0] * (-1 + xi[0] + xi[1] + xi[2]);
69  case 2:
70  return xi[0] * (-1 + 2 * xi[0]);
71  case 3:
72  return -4 * xi[1] * (-1 + xi[0] + xi[1] + xi[2]);
73  case 4:
74  return 4 * xi[0] * xi[1];
75  case 5:
76  return xi[1] * (-1 + 2 * xi[1]);
77  case 6:
78  return -4 * xi[2] * (-1 + xi[0] + xi[1] + xi[2]);
79  case 7:
80  return 4 * xi[0] * xi[2];
81  case 8:
82  return 4 * xi[1] * xi[2];
83  case 9:
84  return xi[2] * (-1 + 2 * xi[2]);
85  }
86  }
87 
88  if constexpr (p == 3) {
89  constexpr double sqrt5 = 2.23606797749978981;
90  switch (i) {
91  case 0:
92  return -((-1 + xi[0] + xi[1] + xi[2]) *
93  (1 + 5 * xi[0] * xi[0] + 5 * xi[1] * xi[1] + 5 * (-1 + xi[2]) * xi[2] + xi[1] * (-5 + 11 * xi[2]) +
94  xi[0] * (-5 + 11 * xi[1] + 11 * xi[2])));
95  case 1:
96  return (5 * xi[0] * (-1 + xi[0] + xi[1] + xi[2]) *
97  (-1 - sqrt5 + 2 * sqrt5 * xi[0] + (3 + sqrt5) * xi[1] + (3 + sqrt5) * xi[2])) /
98  2.;
99  case 2:
100  return (-5 * xi[0] * (-1 + xi[0] + xi[1] + xi[2]) *
101  (1 - sqrt5 + 2 * sqrt5 * xi[0] + (-3 + sqrt5) * xi[1] + (-3 + sqrt5) * xi[2])) /
102  2.;
103  case 3:
104  return xi[0] * (1 + 5 * xi[0] * xi[0] + xi[1] - xi[1] * xi[1] + xi[2] - xi[1] * xi[2] - xi[2] * xi[2] -
105  xi[0] * (5 + xi[1] + xi[2]));
106  case 4:
107  return (5 * xi[1] * (-1 + xi[0] + xi[1] + xi[2]) *
108  (-1 - sqrt5 + (3 + sqrt5) * xi[0] + 2 * sqrt5 * xi[1] + (3 + sqrt5) * xi[2])) /
109  2.;
110  case 5:
111  return -27 * xi[0] * xi[1] * (-1 + xi[0] + xi[1] + xi[2]);
112  case 6:
113  return (5 * xi[0] * xi[1] * (-2 + (3 + sqrt5) * xi[0] - (-3 + sqrt5) * xi[1])) / 2.;
114  case 7:
115  return (-5 * xi[1] * (-1 + xi[0] + xi[1] + xi[2]) *
116  (1 - sqrt5 + (-3 + sqrt5) * xi[0] + 2 * sqrt5 * xi[1] + (-3 + sqrt5) * xi[2])) /
117  2.;
118  case 8:
119  return (-5 * xi[0] * xi[1] * (2 + (-3 + sqrt5) * xi[0] - (3 + sqrt5) * xi[1])) / 2.;
120  case 9:
121  return xi[1] * (1 - xi[0] * xi[0] + 5 * xi[1] * xi[1] + xi[2] - xi[2] * xi[2] - xi[1] * (5 + xi[2]) -
122  xi[0] * (-1 + xi[1] + xi[2]));
123  case 10:
124  return (5 * xi[2] * (-1 + xi[0] + xi[1] + xi[2]) *
125  (-439204 - 196418 * sqrt5 + (710647 + 317811 * sqrt5) * xi[0] + (710647 + 317811 * sqrt5) * xi[1] +
126  606965 * xi[2] + 271443 * sqrt5 * xi[2])) /
127  (271443 + 121393 * sqrt5);
128  case 11:
129  return -27 * xi[0] * xi[2] * (-1 + xi[0] + xi[1] + xi[2]);
130  case 12:
131  return (5 * xi[0] * xi[2] * (-5 - 3 * sqrt5 + (15 + 7 * sqrt5) * xi[0] + 2 * sqrt5 * xi[2])) /
132  (5 + 3 * sqrt5);
133  case 13:
134  return -27 * xi[1] * xi[2] * (-1 + xi[0] + xi[1] + xi[2]);
135  case 14:
136  return 27 * xi[0] * xi[1] * xi[2];
137  case 15:
138  return (5 * xi[1] * xi[2] * (-5 - 3 * sqrt5 + (15 + 7 * sqrt5) * xi[1] + 2 * sqrt5 * xi[2])) /
139  (5 + 3 * sqrt5);
140  case 16:
141  return (5 * xi[2] * (-1 + xi[0] + xi[1] + xi[2]) *
142  (88555 + 39603 * sqrt5 + (54730 + 24476 * sqrt5) * xi[0] + (54730 + 24476 * sqrt5) * xi[1] -
143  5 * (64079 + 28657 * sqrt5) * xi[2])) /
144  (143285 + 64079 * sqrt5);
145  case 17:
146  return (-5 * xi[0] * xi[2] * (2 + (-3 + sqrt5) * xi[0] - (3 + sqrt5) * xi[2])) / 2.;
147  case 18:
148  return (-5 * xi[1] * xi[2] * (2 + (-3 + sqrt5) * xi[1] - (3 + sqrt5) * xi[2])) / 2.;
149  case 19:
150  return -(xi[2] * (-1 + xi[0] * xi[0] + xi[1] * xi[1] + xi[1] * (-1 + xi[2]) - 5 * (-1 + xi[2]) * xi[2] +
151  xi[0] * (-1 + xi[1] + xi[2])));
152  }
153  }
154 
155  return 0.0;
156  }
157 
158  SERAC_HOST_DEVICE static constexpr tensor<double, dim> shape_function_gradient(tensor<double, dim> xi, int i)
159  {
160  if (p == 0) {
161  return {0.0, 0.0, 0.0};
162  }
163  if (p == 1) {
164  switch (i) {
165  case 0:
166  return {-1, -1, -1};
167  case 1:
168  return {1, 0, 0};
169  case 2:
170  return {0, 1, 0};
171  case 3:
172  return {0, 0, 1};
173  }
174  }
175  if (p == 2) {
176  switch (i) {
177  case 0:
178  return {-3 + 4 * xi[0] + 4 * xi[1] + 4 * xi[2], -3 + 4 * xi[0] + 4 * xi[1] + 4 * xi[2],
179  -3 + 4 * xi[0] + 4 * xi[1] + 4 * xi[2]};
180  case 1:
181  return {-4 * (-1 + 2 * xi[0] + xi[1] + xi[2]), -4 * xi[0], -4 * xi[0]};
182  case 2:
183  return {-1 + 4 * xi[0], 0, 0};
184  case 3:
185  return {-4 * xi[1], -4 * (-1 + xi[0] + 2 * xi[1] + xi[2]), -4 * xi[1]};
186  case 4:
187  return {4 * xi[1], 4 * xi[0], 0};
188  case 5:
189  return {0, -1 + 4 * xi[1], 0};
190  case 6:
191  return {-4 * xi[2], -4 * xi[2], -4 * (-1 + xi[0] + xi[1] + 2 * xi[2])};
192  case 7:
193  return {4 * xi[2], 0, 4 * xi[0]};
194  case 8:
195  return {0, 4 * xi[2], 4 * xi[1]};
196  case 9:
197  return {0, 0, -1 + 4 * xi[2]};
198  }
199  }
200 
201  if (p == 3) {
202  constexpr double sqrt5 = 2.23606797749978981;
203  switch (i) {
204  case 0:
205  return {-6 - 15 * xi[0] * xi[0] - 16 * xi[1] * xi[1] + xi[1] * (21 - 33 * xi[2]) + (21 - 16 * xi[2]) * xi[2] -
206  4 * xi[0] * (-5 + 8 * xi[1] + 8 * xi[2]),
207  -6 - 16 * xi[0] * xi[0] + 20 * xi[1] + xi[0] * (21 - 32 * xi[1] - 33 * xi[2]) + 21 * xi[2] -
208  (3 * xi[1] + 4 * xi[2]) * (5 * xi[1] + 4 * xi[2]),
209  -6 - 16 * xi[0] * xi[0] + 21 * xi[1] + xi[0] * (21 - 33 * xi[1] - 32 * xi[2]) + 20 * xi[2] -
210  (4 * xi[1] + 3 * xi[2]) * (4 * xi[1] + 5 * xi[2])};
211  case 1:
212  return {(5 * (6 * sqrt5 * xi[0] * xi[0] +
213  xi[0] * (-2 - 6 * sqrt5 + 6 * (1 + sqrt5) * xi[1] + 6 * (1 + sqrt5) * xi[2]) +
214  (-1 + xi[1] + xi[2]) * (-1 - sqrt5 + (3 + sqrt5) * xi[1] + (3 + sqrt5) * xi[2]))) /
215  2.,
216  (5 * xi[0] *
217  (-4 - 2 * sqrt5 + 3 * (1 + sqrt5) * xi[0] + 2 * (3 + sqrt5) * xi[1] + 2 * (3 + sqrt5) * xi[2])) /
218  2.,
219  (5 * xi[0] *
220  (-4 - 2 * sqrt5 + 3 * (1 + sqrt5) * xi[0] + 2 * (3 + sqrt5) * xi[1] + 2 * (3 + sqrt5) * xi[2])) /
221  2.};
222  case 2:
223  return {-15 * sqrt5 * xi[0] * xi[0] -
224  (5 * (-1 + xi[1] + xi[2]) * (1 - sqrt5 + (-3 + sqrt5) * xi[1] + (-3 + sqrt5) * xi[2])) / 2. -
225  5 * xi[0] * (1 - 3 * sqrt5 + 3 * (-1 + sqrt5) * xi[1] + 3 * (-1 + sqrt5) * xi[2]),
226  (-5 * xi[0] *
227  (4 - 2 * sqrt5 + 3 * (-1 + sqrt5) * xi[0] + 2 * (-3 + sqrt5) * xi[1] + 2 * (-3 + sqrt5) * xi[2])) /
228  2.,
229  (-5 * xi[0] *
230  (4 - 2 * sqrt5 + 3 * (-1 + sqrt5) * xi[0] + 2 * (-3 + sqrt5) * xi[1] + 2 * (-3 + sqrt5) * xi[2])) /
231  2.};
232  case 3:
233  return {1 + 15 * xi[0] * xi[0] + xi[1] - xi[1] * xi[1] + xi[2] - xi[1] * xi[2] - xi[2] * xi[2] -
234  2 * xi[0] * (5 + xi[1] + xi[2]),
235  -(xi[0] * (-1 + xi[0] + 2 * xi[1] + xi[2])), -(xi[0] * (-1 + xi[0] + xi[1] + 2 * xi[2]))};
236  case 4:
237  return {(5 * xi[1] *
238  (-2 * (2 + sqrt5) + 2 * (3 + sqrt5) * xi[0] + 3 * (1 + sqrt5) * xi[1] + 2 * (3 + sqrt5) * xi[2])) /
239  2.,
240  15 * sqrt5 * xi[1] * xi[1] +
241  5 * xi[1] * (-1 - 3 * sqrt5 + 3 * (1 + sqrt5) * xi[0] + 3 * (1 + sqrt5) * xi[2]) +
242  (5 * (-1 + xi[0] + xi[2]) * (-1 - sqrt5 + (3 + sqrt5) * xi[0] + (3 + sqrt5) * xi[2])) / 2.,
243  (5 * xi[1] *
244  (-2 * (2 + sqrt5) + 2 * (3 + sqrt5) * xi[0] + 3 * (1 + sqrt5) * xi[1] + 2 * (3 + sqrt5) * xi[2])) /
245  2.};
246  case 5:
247  return {-27 * xi[1] * (-1 + 2 * xi[0] + xi[1] + xi[2]), -27 * xi[0] * (-1 + xi[0] + 2 * xi[1] + xi[2]),
248  -27 * xi[0] * xi[1]};
249  case 6:
250  return {(-5 * xi[1] * (2 - 2 * (3 + sqrt5) * xi[0] + (-3 + sqrt5) * xi[1])) / 2.,
251  (5 * xi[0] * (-2 + (3 + sqrt5) * xi[0] - 2 * (-3 + sqrt5) * xi[1])) / 2., 0};
252  case 7:
253  return {(-5 * xi[1] *
254  (4 - 2 * sqrt5 + 2 * (-3 + sqrt5) * xi[0] + 3 * (-1 + sqrt5) * xi[1] + 2 * (-3 + sqrt5) * xi[2])) /
255  2.,
256  -15 * sqrt5 * xi[1] * xi[1] -
257  (5 * (-1 + xi[0] + xi[2]) * (1 - sqrt5 + (-3 + sqrt5) * xi[0] + (-3 + sqrt5) * xi[2])) / 2. -
258  5 * xi[1] * (1 - 3 * sqrt5 + 3 * (-1 + sqrt5) * xi[0] + 3 * (-1 + sqrt5) * xi[2]),
259  (-5 * xi[1] *
260  (4 - 2 * sqrt5 + 2 * (-3 + sqrt5) * xi[0] + 3 * (-1 + sqrt5) * xi[1] + 2 * (-3 + sqrt5) * xi[2])) /
261  2.};
262  case 8:
263  return {(5 * xi[1] * (-2 - 2 * (-3 + sqrt5) * xi[0] + (3 + sqrt5) * xi[1])) / 2.,
264  (-5 * xi[0] * (2 + (-3 + sqrt5) * xi[0] - 2 * (3 + sqrt5) * xi[1])) / 2., 0};
265  case 9:
266  return {-(xi[1] * (-1 + 2 * xi[0] + xi[1] + xi[2])),
267  1 - xi[0] * xi[0] + 15 * xi[1] * xi[1] + xi[2] - xi[2] * xi[2] - 2 * xi[1] * (5 + xi[2]) -
268  xi[0] * (-1 + 2 * xi[1] + xi[2]),
269  -(xi[1] * (-1 + xi[0] + xi[1] + 2 * xi[2]))};
270  case 10:
271  return {
272  (5 * xi[2] *
273  (-2 * (2 + sqrt5) + 2 * (3 + sqrt5) * xi[0] + 2 * (3 + sqrt5) * xi[1] + 3 * (1 + sqrt5) * xi[2])) /
274  2.,
275  (5 * xi[2] *
276  (-2 * (2 + sqrt5) + 2 * (3 + sqrt5) * xi[0] + 2 * (3 + sqrt5) * xi[1] + 3 * (1 + sqrt5) * xi[2])) /
277  2.,
278  (5 * (1 + sqrt5 + (3 + sqrt5) * xi[0] * xi[0] - 2 * (2 + sqrt5) * xi[1] + (3 + sqrt5) * xi[1] * xi[1] +
279  6 * (1 + sqrt5) * xi[1] * xi[2] + 2 * xi[2] * (-1 - 3 * sqrt5 + 3 * sqrt5 * xi[2]) +
280  2 * xi[0] * (-2 - sqrt5 + (3 + sqrt5) * xi[1] + 3 * (1 + sqrt5) * xi[2]))) /
281  2.};
282  case 11:
283  return {-27 * xi[2] * (-1 + 2 * xi[0] + xi[1] + xi[2]), -27 * xi[0] * xi[2],
284  -27 * xi[0] * (-1 + xi[0] + xi[1] + 2 * xi[2])};
285  case 12:
286  return {(-5 * xi[2] * (2 - 2 * (3 + sqrt5) * xi[0] + (-3 + sqrt5) * xi[2])) / 2., 0,
287  (5 * xi[0] * (-2 + (3 + sqrt5) * xi[0] - 2 * (-3 + sqrt5) * xi[2])) / 2.};
288  case 13:
289  return {-27 * xi[1] * xi[2], -27 * xi[2] * (-1 + xi[0] + 2 * xi[1] + xi[2]),
290  -27 * xi[1] * (-1 + xi[0] + xi[1] + 2 * xi[2])};
291  case 14:
292  return {27 * xi[1] * xi[2], 27 * xi[0] * xi[2], 27 * xi[0] * xi[1]};
293  case 15:
294  return {0, (-5 * xi[2] * (2 - 2 * (3 + sqrt5) * xi[1] + (-3 + sqrt5) * xi[2])) / 2.,
295  (5 * xi[1] * (-2 + (3 + sqrt5) * xi[1] - 2 * (-3 + sqrt5) * xi[2])) / 2.};
296  case 16:
297  return {
298  (-5 * xi[2] *
299  (4 - 2 * sqrt5 + 2 * (-3 + sqrt5) * xi[0] + 2 * (-3 + sqrt5) * xi[1] + 3 * (-1 + sqrt5) * xi[2])) /
300  2.,
301  (-5 * xi[2] *
302  (4 - 2 * sqrt5 + 2 * (-3 + sqrt5) * xi[0] + 2 * (-3 + sqrt5) * xi[1] + 3 * (-1 + sqrt5) * xi[2])) /
303  2.,
304  (-5 * (-3 + sqrt5) * xi[0] * xi[0]) / 2. -
305  5 * xi[0] * (2 - sqrt5 + (-3 + sqrt5) * xi[1] + 3 * (-1 + sqrt5) * xi[2]) -
306  (5 * (-1 + sqrt5 + (-3 + sqrt5) * xi[1] * xi[1] + 2 * xi[2] * (1 - 3 * sqrt5 + 3 * sqrt5 * xi[2]) +
307  xi[1] * (4 - 2 * sqrt5 + 6 * (-1 + sqrt5) * xi[2]))) /
308  2.};
309  case 17:
310  return {(5 * xi[2] * (-2 - 2 * (-3 + sqrt5) * xi[0] + (3 + sqrt5) * xi[2])) / 2., 0,
311  (-5 * xi[0] * (2 + (-3 + sqrt5) * xi[0] - 2 * (3 + sqrt5) * xi[2])) / 2.};
312  case 18:
313  return {0, (5 * xi[2] * (-2 - 2 * (-3 + sqrt5) * xi[1] + (3 + sqrt5) * xi[2])) / 2.,
314  (-5 * xi[1] * (2 + (-3 + sqrt5) * xi[1] - 2 * (3 + sqrt5) * xi[2])) / 2.};
315  case 19:
316  return {-(xi[2] * (-1 + 2 * xi[0] + xi[1] + xi[2])), -(xi[2] * (-1 + xi[0] + 2 * xi[1] + xi[2])),
317  1 + xi[0] - xi[0] * xi[0] + xi[1] - xi[0] * xi[1] - xi[1] * xi[1] - 2 * (5 + xi[0] + xi[1]) * xi[2] +
318  15 * xi[2] * xi[2]};
319  }
320  }
321 
322  return {};
323  }
324 
325  SERAC_HOST_DEVICE static constexpr tensor<double, ndof> shape_functions(tensor<double, dim> xi)
326  {
327  tensor<double, ndof> output{};
328  for (int i = 0; i < ndof; i++) {
329  output[i] = shape_function(xi, i);
330  }
331  return output;
332  }
333 
334  SERAC_HOST_DEVICE static constexpr tensor<double, ndof, dim> shape_function_gradients(tensor<double, dim> xi)
335  {
336  tensor<double, ndof, dim> output{};
337  for (int i = 0; i < ndof; i++) {
338  output[i] = shape_function_gradient(xi, i);
339  }
340  return output;
341  }
342 
343  template <typename in_t, int q>
344  static auto batch_apply_shape_fn(int j, tensor<in_t, nqpts(q)> input, const TensorProductQuadratureRule<q>&)
345  {
346  using source_t = decltype(get<0>(get<0>(in_t{})) + dot(get<1>(get<0>(in_t{})), tensor<double, dim>{}));
347  using flux_t = decltype(get<0>(get<1>(in_t{})) + dot(get<1>(get<1>(in_t{})), tensor<double, dim>{}));
348 
349  constexpr auto xi = GaussLegendreNodes<q, mfem::Geometry::TETRAHEDRON>();
350 
351  tensor<tuple<source_t, flux_t>, nqpts(q)> output;
352 
353  for (int i = 0; i < nqpts(q); i++) {
354  double phi_j = shape_function(xi[i], j);
355  tensor<double, dim> dphi_j_dxi = shape_function_gradient(xi[i], j);
356 
357  const auto& d00 = get<0>(get<0>(input(i)));
358  const auto& d01 = get<1>(get<0>(input(i)));
359  const auto& d10 = get<0>(get<1>(input(i)));
360  const auto& d11 = get<1>(get<1>(input(i)));
361 
362  output[i] = {d00 * phi_j + dot(d01, dphi_j_dxi), d10 * phi_j + dot(d11, dphi_j_dxi)};
363  }
364 
365  return output;
366  }
367 
368  template <int q>
369  SERAC_HOST_DEVICE static auto interpolate(const tensor<double, c, ndof>& X, const TensorProductQuadratureRule<q>&)
370  {
371  constexpr auto xi = GaussLegendreNodes<q, mfem::Geometry::TETRAHEDRON>();
372 
373  // transpose the quadrature data into a flat tensor of tuples
374  union {
375  tensor<tuple<tensor<double, c>, tensor<double, c, dim> >, nqpts(q)> unflattened;
376  tensor<qf_input_type, nqpts(q)> flattened;
377  } output{};
378 
379  for (int i = 0; i < c; i++) {
380  for (int j = 0; j < nqpts(q); j++) {
381  for (int k = 0; k < ndof; k++) {
382  get<VALUE>(output.unflattened[j])[i] += X(i, k) * shape_function(xi[j], k);
383  get<GRADIENT>(output.unflattened[j])[i] += X(i, k) * shape_function_gradient(xi[j], k);
384  }
385  }
386  }
387 
388  return output.flattened;
389  }
390 
391  template <typename source_type, typename flux_type, int q>
392  SERAC_HOST_DEVICE static void integrate(const tensor<tuple<source_type, flux_type>, nqpts(q)>& qf_output,
393  const TensorProductQuadratureRule<q>&,
394  tensor<double, c, ndof>* element_residual, int step = 1)
395  {
396  if constexpr (is_zero<source_type>{} && is_zero<flux_type>{}) {
397  return;
398  }
399 
400  using source_component_type = std::conditional_t<is_zero<source_type>{}, zero, double>;
401  using flux_component_type = std::conditional_t<is_zero<flux_type>{}, zero, tensor<double, dim> >;
402 
403  constexpr int ntrial = std::max(size(source_type{}), size(flux_type{}) / dim) / c;
404  constexpr auto integration_points = GaussLegendreNodes<q, mfem::Geometry::TETRAHEDRON>();
405  constexpr auto integration_weights = GaussLegendreWeights<q, mfem::Geometry::TETRAHEDRON>();
406 
407  for (int j = 0; j < ntrial; j++) {
408  for (int i = 0; i < c; i++) {
409  for (int Q = 0; Q < nqpts(q); Q++) {
410  tensor<double, dim> xi = integration_points[Q];
411  double wt = integration_weights[Q];
412 
413  source_component_type source;
414  if constexpr (!is_zero<source_type>{}) {
415  source = reinterpret_cast<const double*>(&get<SOURCE>(qf_output[Q]))[i * ntrial + j];
416  }
417 
418  flux_component_type flux;
419  if constexpr (!is_zero<flux_type>{}) {
420  for (int k = 0; k < dim; k++) {
421  flux[k] = reinterpret_cast<const double*>(&get<FLUX>(qf_output[Q]))[(i * dim + k) * ntrial + j];
422  }
423  }
424 
425  for (int k = 0; k < ndof; k++) {
426  element_residual[j * step](i, k) +=
427  (source * shape_function(xi, k) + dot(flux, shape_function_gradient(xi, k))) * wt;
428  }
429  }
430  }
431  }
432  }
433 };
SERAC_HOST_DEVICE
#define SERAC_HOST_DEVICE
Macro that evaluates to __host__ __device__ when compiling with nvcc and does nothing on a host compi...
Definition: accelerator.hpp:38
serac::num_quadrature_points
constexpr int num_quadrature_points(mfem::Geometry::Type g, int Q)
return the number of quadrature points in a Gauss-Legendre rule with parameter "Q"
Definition: geometry.hpp:31
serac::dot
constexpr SERAC_HOST_DEVICE auto dot(const isotropic_tensor< S, m, m > &I, const tensor< T, m, n... > &A)
dot product between an isotropic and (nonisotropic) tensor
Definition: isotropic_tensor.hpp:203
serac::tensor
tensor(const T(&data)[n1]) -> tensor< T, n1 >
class template argument deduction guide for type tensor.
serac::max
SERAC_HOST_DEVICE auto max(dual< gradient_type > a, double b)
Implementation of max for dual numbers.
Definition: dual.hpp:229
serac::size
constexpr SERAC_HOST_DEVICE int size(const tensor< T, n... > &)
returns the total number of stored values in a tensor
Definition: tensor.hpp:1934
serac::type
constexpr SERAC_HOST_DEVICE auto type(const tuple< T... > &values)
a function intended to be used for extracting the ith type from a tuple.
Definition: tuple.hpp:379
Generated by doxygen 1.9.1