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