Basix 0.6.0

Home     Installation     Demos     C++ docs     Python docs

maps.h
1 // Copyright (c) 2021-2022 Matthew Scroggs and Garth N. Wells
2 // FEniCS Project
3 // SPDX-License-Identifier: MIT
4 
5 #pragma once
6 
7 #include "mdspan.hpp"
8 #include <stdexcept>
9 #include <type_traits>
10 
12 namespace basix::maps
13 {
14 
16 enum class type
17 {
18  identity = 0,
19  L2Piola = 1,
20  covariantPiola = 2,
21  contravariantPiola = 3,
22  doubleCovariantPiola = 4,
23  doubleContravariantPiola = 5,
24 };
25 
27 template <typename O, typename P, typename Q, typename R>
28 void l2_piola(O&& r, const P& U, const Q& /*J*/, double detJ, const R& /*K*/)
29 {
30  assert(U.extent(0) == r.extent(0));
31  assert(U.extent(1) == r.extent(1));
32  for (std::size_t i = 0; i < U.extent(0); ++i)
33  for (std::size_t j = 0; j < U.extent(1); ++j)
34  r(i, j) = U(i, j) / detJ;
35 }
36 
38 template <typename O, typename P, typename Q, typename R>
39 void covariant_piola(O&& r, const P& U, const Q& /*J*/, double /*detJ*/,
40  const R& K)
41 {
42  using T = typename std::decay_t<O>::value_type;
43  for (std::size_t p = 0; p < U.extent(0); ++p)
44  {
45  // r_p = K^T U_p, where p indicates the p-th row
46  for (std::size_t i = 0; i < r.extent(1); ++i)
47  {
48  T acc = 0;
49  for (std::size_t k = 0; k < K.extent(0); ++k)
50  acc += K(k, i) * U(p, k);
51  r(p, i) = acc;
52  }
53  }
54 }
55 
57 template <typename O, typename P, typename Q, typename R>
58 void contravariant_piola(O&& r, const P& U, const Q& J, double detJ,
59  const R& /*K*/)
60 {
61  using T = typename std::decay_t<O>::value_type;
62  for (std::size_t p = 0; p < U.extent(0); ++p)
63  {
64  for (std::size_t i = 0; i < r.extent(1); ++i)
65  {
66  T acc = 0;
67  for (std::size_t k = 0; k < J.extent(1); ++k)
68  acc += J(i, k) * U(p, k);
69  r(p, i) = acc;
70  }
71  }
72 
73  std::transform(r.data_handle(), r.data_handle() + r.size(), r.data_handle(),
74  [detJ](auto ri) { return ri / detJ; });
75 }
76 
78 template <typename O, typename P, typename Q, typename R>
79 void double_covariant_piola(O&& r, const P& U, const Q& J, double /*detJ*/,
80  const R& K)
81 {
82  namespace stdex = std::experimental;
83  using T = typename std::decay_t<O>::value_type;
84  for (std::size_t p = 0; p < U.extent(0); ++p)
85  {
86  stdex::mdspan<const T, stdex::dextents<std::size_t, 2>> _U(
87  U.data_handle() + p * U.extent(1), J.extent(1), J.extent(1));
88  stdex::mdspan<T, stdex::dextents<std::size_t, 2>> _r(
89  r.data_handle() + p * r.extent(1), K.extent(1), K.extent(1));
90  // _r = K^T _U K
91  for (std::size_t i = 0; i < _r.extent(0); ++i)
92  {
93  for (std::size_t j = 0; j < _r.extent(1); ++j)
94  {
95  T acc = 0;
96  for (std::size_t k = 0; k < K.extent(0); ++k)
97  for (std::size_t l = 0; l < _U.extent(1); ++l)
98  acc += K(k, i) * _U(k, l) * K(l, j);
99  _r(i, j) = acc;
100  }
101  }
102  }
103 }
104 
106 template <typename O, typename P, typename Q, typename R>
107 void double_contravariant_piola(O&& r, const P& U, const Q& J, double detJ,
108  const R& /*K*/)
109 {
110  namespace stdex = std::experimental;
111  using T = typename std::decay_t<O>::value_type;
112 
113  for (std::size_t p = 0; p < U.extent(0); ++p)
114  {
115  stdex::mdspan<const T, stdex::dextents<std::size_t, 2>> _U(
116  U.data_handle() + p * U.extent(1), J.extent(1), J.extent(1));
117  stdex::mdspan<T, stdex::dextents<std::size_t, 2>> _r(
118  r.data_handle() + p * r.extent(1), J.extent(0), J.extent(0));
119 
120  // _r = J U J^T
121  for (std::size_t i = 0; i < _r.extent(0); ++i)
122  {
123  for (std::size_t j = 0; j < _r.extent(1); ++j)
124  {
125  T acc = 0;
126  for (std::size_t k = 0; k < J.extent(1); ++k)
127  for (std::size_t l = 0; l < _U.extent(1); ++l)
128  acc += J(i, k) * _U(k, l) * J(j, l);
129  _r(i, j) = acc;
130  }
131  }
132  }
133 
134  std::transform(r.data_handle(), r.data_handle() + r.size(), r.data_handle(),
135  [detJ](auto ri) { return ri / (detJ * detJ); });
136 }
137 
138 } // namespace basix::maps
Information about finite element maps.
Definition: maps.h:13
void l2_piola(O &&r, const P &U, const Q &, double detJ, const R &)
L2 Piola map.
Definition: maps.h:28
void covariant_piola(O &&r, const P &U, const Q &, double, const R &K)
Covariant Piola map.
Definition: maps.h:39
void contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Contravariant Piola map.
Definition: maps.h:58
void double_contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Double contravariant Piola map.
Definition: maps.h:107
void double_covariant_piola(O &&r, const P &U, const Q &J, double, const R &K)
Double covariant Piola map.
Definition: maps.h:79
type
Map type.
Definition: maps.h:17