Smith  0.1
Smith is an implicit thermal structural mechanics simulation code.
triangle_L2.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}, {1,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::TRIANGLE, L2<p, c> > {
23  static constexpr auto geometry = mfem::Geometry::TRIANGLE;
24  static constexpr auto family = Family::L2;
25  static constexpr int components = c;
26  static constexpr int dim = 2;
27  static constexpr int n = (p + 1);
28  static constexpr int ndof = (p + 1) * (p + 2) / 2;
29 
30  static constexpr int VALUE = 0, GRADIENT = 1;
31  static constexpr int SOURCE = 0, FLUX = 1;
32 
33  using residual_type =
34  typename std::conditional<components == 1, tensor<double, ndof>, tensor<double, ndof, components> >::type;
35 
36  using dof_type = tensor<double, c, ndof>;
37  using dof_type_if = tensor<double, c, 2, 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  /*
45 
46  interpolation nodes and their associated numbering:
47 
48  linear
49  2
50  | .
51  | .
52  | .
53  | .
54  | .
55  0-----------1
56 
57  quadratic
58  5
59  | .
60  | .
61  3 4
62  | .
63  | .
64  0-----1-----2
65 
66 
67  cubic
68  9
69  | .
70  7 8
71  | .
72  4 5 6
73  | .
74  0---1---2---3
75 
76  */
77 
78  SMITH_HOST_DEVICE static constexpr double shape_function([[maybe_unused]] tensor<double, dim> xi,
79  [[maybe_unused]] int i)
80  {
81  // constant
82  if constexpr (n == 1) {
83  return 1;
84  }
85 
86  // linear
87  if constexpr (n == 2) {
88  switch (i) {
89  case 0:
90  return 1 - xi[0] - xi[1];
91  case 1:
92  return xi[0];
93  case 2:
94  return xi[1];
95  }
96  }
97 
98  // quadratic
99  if constexpr (n == 3) {
100  switch (i) {
101  case 0:
102  return (-1 + xi[0] + xi[1]) * (-1 + 2 * xi[0] + 2 * xi[1]);
103  case 1:
104  return -4 * xi[0] * (-1 + xi[0] + xi[1]);
105  case 2:
106  return xi[0] * (-1 + 2 * xi[0]);
107  case 3:
108  return -4 * xi[1] * (-1 + xi[0] + xi[1]);
109  case 4:
110  return 4 * xi[0] * xi[1];
111  case 5:
112  return xi[1] * (-1 + 2 * xi[1]);
113  }
114  }
115 
116  // cubic
117  if constexpr (n == 4) {
118  constexpr double sqrt5 = 2.23606797749978981;
119  switch (i) {
120  case 0:
121  return -((-1 + xi[0] + xi[1]) *
122  (1 + 5 * xi[0] * xi[0] + 5 * (-1 + xi[1]) * xi[1] + xi[0] * (-5 + 11 * xi[1])));
123  case 1:
124  return (5 * xi[0] * (-1 + xi[0] + xi[1]) * (-1 - sqrt5 + 2 * sqrt5 * xi[0] + (3 + sqrt5) * xi[1])) / 2.;
125  case 2:
126  return (-5 * xi[0] * (-1 + xi[0] + xi[1]) * (1 - sqrt5 + 2 * sqrt5 * xi[0] + (-3 + sqrt5) * xi[1])) / 2.;
127  case 3:
128  return xi[0] * (1 + 5 * xi[0] * xi[0] + xi[1] - xi[1] * xi[1] - xi[0] * (5 + xi[1]));
129  case 4:
130  return (5 * xi[1] * (-1 + xi[0] + xi[1]) * (-1 - sqrt5 + (3 + sqrt5) * xi[0] + 2 * sqrt5 * xi[1])) / 2.;
131  case 5:
132  return -27 * xi[0] * xi[1] * (-1 + xi[0] + xi[1]);
133  case 6:
134  return (5 * xi[0] * xi[1] * (-2 + (3 + sqrt5) * xi[0] - (-3 + sqrt5) * xi[1])) / 2.;
135  case 7:
136  return (5 * xi[1] * (-1 + xi[0] + xi[1]) *
137  (5 - 3 * sqrt5 + 2 * (-5 + 2 * sqrt5) * xi[0] + 5 * (-1 + sqrt5) * xi[1])) /
138  (-5 + sqrt5);
139  case 8:
140  return (-5 * xi[0] * xi[1] * (2 + (-3 + sqrt5) * xi[0] - (3 + sqrt5) * xi[1])) / 2.;
141  case 9:
142  return xi[1] * (1 + xi[0] - xi[0] * xi[0] - xi[0] * xi[1] + 5 * (-1 + xi[1]) * xi[1]);
143  }
144  }
145 
146  return 0.0;
147  }
148 
149  SMITH_HOST_DEVICE static constexpr tensor<double, dim> shape_function_gradient(
150  [[maybe_unused]] tensor<double, dim> xi, [[maybe_unused]] int i)
151  {
152  // constant
153  if constexpr (n == 1) {
154  return {0.0, 0.0};
155  }
156 
157  // linear
158  if constexpr (n == 2) {
159  switch (i) {
160  case 0:
161  return {-1.0, -1.0};
162  case 1:
163  return {1.0, 0.0};
164  case 2:
165  return {0.0, 1.0};
166  }
167  }
168 
169  // quadratic
170  if constexpr (n == 3) {
171  switch (i) {
172  case 0:
173  return {-3 + 4 * xi[0] + 4 * xi[1], -3 + 4 * xi[0] + 4 * xi[1]};
174  case 1:
175  return {-4 * (-1 + 2 * xi[0] + xi[1]), -4 * xi[0]};
176  case 2:
177  return {-1 + 4 * xi[0], 0};
178  case 3:
179  return {-4 * xi[1], -4 * (-1 + xi[0] + 2 * xi[1])};
180  case 4:
181  return {4 * xi[1], 4 * xi[0]};
182  case 5:
183  return {0, -1 + 4 * xi[1]};
184  }
185  }
186 
187  // cubic
188  if constexpr (n == 4) {
189  constexpr double sqrt5 = 2.23606797749978981;
190  switch (i) {
191  case 0:
192  return {-6 - 15 * xi[0] * xi[0] + 4 * xi[0] * (5 - 8 * xi[1]) + (21 - 16 * xi[1]) * xi[1],
193  -6 - 16 * xi[0] * xi[0] + xi[0] * (21 - 32 * xi[1]) + 5 * (4 - 3 * xi[1]) * xi[1]};
194  case 1:
195  return {(5 * (6 * sqrt5 * xi[0] * xi[0] + xi[0] * (-2 - 6 * sqrt5 + 6 * (1 + sqrt5) * xi[1]) +
196  (-1 + xi[1]) * (-1 - sqrt5 + (3 + sqrt5) * xi[1]))) /
197  2.,
198  (5 * xi[0] * (-2 * (2 + sqrt5) + 3 * (1 + sqrt5) * xi[0] + 2 * (3 + sqrt5) * xi[1])) / 2.};
199  case 2:
200  return {(-5 * (6 * sqrt5 * xi[0] * xi[0] + (-1 + xi[1]) * (1 - sqrt5 + (-3 + sqrt5) * xi[1]) +
201  xi[0] * (2 - 6 * sqrt5 + 6 * (-1 + sqrt5) * xi[1]))) /
202  2.,
203  (-5 * xi[0] * (4 - 2 * sqrt5 + 3 * (-1 + sqrt5) * xi[0] + 2 * (-3 + sqrt5) * xi[1])) / 2.};
204  case 3:
205  return {1 + 15 * xi[0] * xi[0] + xi[1] - xi[1] * xi[1] - 2 * xi[0] * (5 + xi[1]),
206  -(xi[0] * (-1 + xi[0] + 2 * xi[1]))};
207  case 4:
208  return {(5 * xi[1] * (-2 * (2 + sqrt5) + 2 * (3 + sqrt5) * xi[0] + 3 * (1 + sqrt5) * xi[1])) / 2.,
209  (5 * (1 + sqrt5 - 2 * (2 + sqrt5) * xi[0] + (3 + sqrt5) * xi[0] * xi[0] +
210  6 * (1 + sqrt5) * xi[0] * xi[1] + 2 * xi[1] * (-1 - 3 * sqrt5 + 3 * sqrt5 * xi[1]))) /
211  2.};
212  case 5:
213  return {-27 * xi[1] * (-1 + 2 * xi[0] + xi[1]), -27 * xi[0] * (-1 + xi[0] + 2 * xi[1])};
214  case 6:
215  return {(-5 * xi[1] * (2 - 2 * (3 + sqrt5) * xi[0] + (-3 + sqrt5) * xi[1])) / 2.,
216  (5 * xi[0] * (-2 + (3 + sqrt5) * xi[0] - 2 * (-3 + sqrt5) * xi[1])) / 2.};
217  case 7:
218  return {(-5 * xi[1] * (4 - 2 * sqrt5 + 2 * (-3 + sqrt5) * xi[0] + 3 * (-1 + sqrt5) * xi[1])) / 2.,
219  (-5 * (-1 + sqrt5 + (-3 + sqrt5) * xi[0] * xi[0] + 2 * xi[1] * (1 - 3 * sqrt5 + 3 * sqrt5 * xi[1]) +
220  xi[0] * (4 - 2 * sqrt5 + 6 * (-1 + sqrt5) * xi[1]))) /
221  2.};
222  case 8:
223  return {(5 * xi[1] * (-2 - 2 * (-3 + sqrt5) * xi[0] + (3 + sqrt5) * xi[1])) / 2.,
224  (-5 * xi[0] * (2 + (-3 + sqrt5) * xi[0] - 2 * (3 + sqrt5) * xi[1])) / 2.};
225  case 9:
226  return {-(xi[1] * (-1 + 2 * xi[0] + xi[1])),
227  1 + xi[0] - xi[0] * xi[0] - 2 * (5 + xi[0]) * xi[1] + 15 * xi[1] * xi[1]};
228  }
229  }
230 
231  return {};
232  }
233 
234  SMITH_HOST_DEVICE static constexpr tensor<double, ndof> shape_functions(tensor<double, dim> xi)
235  {
236  tensor<double, ndof> output{};
237  for (int i = 0; i < ndof; i++) {
238  output[i] = shape_function(xi, i);
239  }
240  return output;
241  }
242 
243  SMITH_HOST_DEVICE static constexpr tensor<double, ndof, dim> shape_function_gradients(tensor<double, dim> xi)
244  {
245  tensor<double, ndof, dim> output{};
246  for (int i = 0; i < ndof; i++) {
247  output[i] = shape_function_gradient(xi, i);
248  }
249  return output;
250  }
251 
252  template <typename in_t, int q>
253  static auto batch_apply_shape_fn(int j, const tensor<in_t, q*(q + 1) / 2>& input,
254  const TensorProductQuadratureRule<q>&)
255  {
256  using source_t = decltype(get<0>(get<0>(in_t{})) + dot(get<1>(get<0>(in_t{})), tensor<double, 2>{}));
257  using flux_t = decltype(get<0>(get<1>(in_t{})) + dot(get<1>(get<1>(in_t{})), tensor<double, 2>{}));
258 
259  constexpr auto xi = GaussLegendreNodes<q, mfem::Geometry::TRIANGLE>();
260 
261  static constexpr int Q = q * (q + 1) / 2;
262  tensor<tuple<source_t, flux_t>, Q> output;
263 
264  for (int i = 0; i < Q; i++) {
265  double phi_j = shape_function(xi[i], j);
266  tensor<double, dim> dphi_j_dxi = shape_function_gradient(xi[i], j);
267 
268  const auto& d00 = get<0>(get<0>(input(i)));
269  const auto& d01 = get<1>(get<0>(input(i)));
270  const auto& d10 = get<0>(get<1>(input(i)));
271  const auto& d11 = get<1>(get<1>(input(i)));
272 
273  output[i] = {d00 * phi_j + dot(d01, dphi_j_dxi), d10 * phi_j + dot(d11, dphi_j_dxi)};
274  }
275 
276  return output;
277  }
278 
279  template <typename T, int q>
280  static auto batch_apply_shape_fn_interior_face(int jx, const tensor<T, q*(q + 1) / 2>& input,
281  const TensorProductQuadratureRule<q>&)
282  {
283  constexpr auto xi = GaussLegendreNodes<q, mfem::Geometry::TRIANGLE>();
284 
285  using source0_t = decltype(get<0>(get<0>(T{})) + get<1>(get<0>(T{})));
286  using source1_t = decltype(get<0>(get<1>(T{})) + get<1>(get<1>(T{})));
287 
288  static constexpr int Q = q * (q + 1) / 2;
289  tensor<tuple<source0_t, source1_t>, Q> output;
290 
291  for (int i = 0; i < Q; i++) {
292  int j = jx % ndof;
293  int s = jx / ndof;
294 
295  double phi0_j = shape_function(xi[i], j) * (s == 0);
296  double phi1_j = shape_function(xi[i], j) * (s == 1);
297 
298  const auto& d00 = get<0>(get<0>(input(i)));
299  const auto& d01 = get<1>(get<0>(input(i)));
300  const auto& d10 = get<0>(get<1>(input(i)));
301  const auto& d11 = get<1>(get<1>(input(i)));
302 
303  output[i] = {d00 * phi0_j + d01 * phi1_j, d10 * phi0_j + d11 * phi1_j};
304  }
305 
306  return output;
307  }
308 
309  template <int q>
310  SMITH_HOST_DEVICE static auto interpolate(const tensor<double, c, ndof>& X, const TensorProductQuadratureRule<q>&)
311  {
312  constexpr auto xi = GaussLegendreNodes<q, mfem::Geometry::TRIANGLE>();
313  static constexpr int num_quadrature_points = q * (q + 1) / 2;
314 
315  // transpose the quadrature data into a flat tensor of tuples
316  union {
317  tensor<tuple<tensor<double, c>, tensor<double, c, dim> >, num_quadrature_points> unflattened;
318  tensor<qf_input_type, num_quadrature_points> flattened;
319  } output{};
320 
321  for (int i = 0; i < c; i++) {
322  for (int j = 0; j < num_quadrature_points; j++) {
323  for (int k = 0; k < ndof; k++) {
324  get<VALUE>(output.unflattened[j])[i] += X(i, k) * shape_function(xi[j], k);
325  get<GRADIENT>(output.unflattened[j])[i] += X(i, k) * shape_function_gradient(xi[j], k);
326  }
327  }
328  }
329 
330  return output.flattened;
331  }
332 
333  // overload for two-sided interior face kernels
334  template <int q>
335  SMITH_HOST_DEVICE static auto interpolate(const dof_type_if& X, const TensorProductQuadratureRule<q>&)
336  {
337  constexpr auto xi = GaussLegendreNodes<q, mfem::Geometry::TRIANGLE>();
338  static constexpr int num_quadrature_points = q * (q + 1) / 2;
339 
340  tensor<tuple<value_type, value_type>, num_quadrature_points> output{};
341 
342  // apply the shape functions
343  for (int i = 0; i < c; i++) {
344  for (int j = 0; j < num_quadrature_points; j++) {
345  for (int k = 0; k < ndof; k++) {
346  if constexpr (c == 1) {
347  get<0>(output[j]) += X[i][0][k] * shape_function(xi[j], k);
348  get<1>(output[j]) += X[i][1][k] * shape_function(xi[j], k);
349  } else {
350  get<0>(output[j])[i] += X[i][0][k] * shape_function(xi[j], k);
351  get<1>(output[j])[i] += X[i][1][k] * shape_function(xi[j], k);
352  }
353  }
354  }
355  }
356 
357  return output;
358  }
359 
360  template <typename source_type, typename flux_type, int q>
361  SMITH_HOST_DEVICE static void integrate(const tensor<tuple<source_type, flux_type>, q*(q + 1) / 2>& qf_output,
362  const TensorProductQuadratureRule<q>&,
363  tensor<double, c, ndof>* element_residual, int step = 1)
364  {
365  if constexpr (is_zero<source_type>{} && is_zero<flux_type>{}) {
366  return;
367  }
368 
369  using source_component_type = std::conditional_t<is_zero<source_type>{}, zero, double>;
370  using flux_component_type = std::conditional_t<is_zero<flux_type>{}, zero, tensor<double, dim> >;
371 
372  constexpr int num_quadrature_points = q * (q + 1) / 2;
373  constexpr int ntrial = std::max(size(source_type{}), size(flux_type{}) / dim) / c;
374  constexpr auto integration_points = GaussLegendreNodes<q, mfem::Geometry::TRIANGLE>();
375  constexpr auto integration_weights = GaussLegendreWeights<q, mfem::Geometry::TRIANGLE>();
376 
377  for (int j = 0; j < ntrial; j++) {
378  for (int i = 0; i < c; i++) {
379  for (int Q = 0; Q < num_quadrature_points; Q++) {
380  tensor<double, 2> xi = integration_points[Q];
381  double wt = integration_weights[Q];
382 
383  source_component_type source;
384  if constexpr (!is_zero<source_type>{}) {
385  source = reinterpret_cast<const double*>(&get<SOURCE>(qf_output[Q]))[i * ntrial + j];
386  }
387 
388  flux_component_type flux;
389  if constexpr (!is_zero<flux_type>{}) {
390  for (int k = 0; k < dim; k++) {
391  flux[k] = reinterpret_cast<const double*>(&get<FLUX>(qf_output[Q]))[(i * dim + k) * ntrial + j];
392  }
393  }
394 
395  for (int k = 0; k < ndof; k++) {
396  element_residual[j * step](i, k) +=
397  (source * shape_function(xi, k) + dot(flux, shape_function_gradient(xi, k))) * wt;
398  }
399  }
400  }
401  }
402  }
403 
404  template <typename T, int q>
405  SMITH_HOST_DEVICE static void integrate(const tensor<tuple<T, T>, (q * (q + 1)) / 2>& qf_output,
406  const TensorProductQuadratureRule<q>&, dof_type_if* element_residual,
407  [[maybe_unused]] int step = 1)
408  {
409  constexpr int ntrial = size(T{}) / c;
410  constexpr int num_quadrature_points = (q * (q + 1)) / 2;
411  constexpr auto integration_points = GaussLegendreNodes<q, mfem::Geometry::TRIANGLE>();
412  constexpr auto integration_weights = GaussLegendreWeights<q, mfem::Geometry::TRIANGLE>();
413 
414  for (int j = 0; j < ntrial; j++) {
415  for (int i = 0; i < c; i++) {
416  for (int Q = 0; Q < num_quadrature_points; Q++) {
417  tensor<double, 2> xi = integration_points[Q];
418  double wt = integration_weights[Q];
419 
420  double source_0 = reinterpret_cast<const double*>(&get<0>(qf_output[Q]))[i * ntrial + j];
421  double source_1 = reinterpret_cast<const double*>(&get<1>(qf_output[Q]))[i * ntrial + j];
422 
423  for (int k = 0; k < ndof; k++) {
424  element_residual[j * step](i, 0, k) += (source_0 * shape_function(xi, k)) * wt;
425  element_residual[j * step](i, 1, k) += (source_1 * shape_function(xi, k)) * wt;
426  }
427  }
428  }
429  }
430  }
431 };
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