Note: this is documentation for an old release. View the latest documentation at docs.fenicsproject.org/basix/v0.8.0/cpp/maps_8h_source.html

Basix 0.4.1

Home     Installation     Demos     C++ docs     Python docs

maps.h
1 // Copyright (c) 2021 Matthew Scroggs
2 // FEniCS Project
3 // SPDX-License-Identifier: MIT
4 
5 #pragma once
6 
7 #include <stdexcept>
8 #include <type_traits>
9 #include <xtensor/xadapt.hpp>
10 #include <xtensor/xtensor.hpp>
11 #include <xtensor/xview.hpp>
12 #include <xtl/xspan.hpp>
13 
15 namespace basix::maps
16 {
17 
19 enum class type
20 {
21  identity = 0,
22  L2Piola = 1,
23  covariantPiola = 2,
24  contravariantPiola = 3,
25  doubleCovariantPiola = 4,
26  doubleContravariantPiola = 5,
27 };
28 
29 namespace impl
30 {
31 template <typename O, typename Mat0, typename Mat1, typename Mat2>
32 void dot22(O&& r, const Mat0& A, const Mat1& B, const Mat2& C)
33 {
34  assert(A.shape(1) == B.shape(0));
35  using T = typename std::decay_t<O>::value_type;
36  for (std::size_t i = 0; i < r.shape(0); ++i)
37  for (std::size_t j = 0; j < r.shape(1); ++j)
38  {
39  T acc = 0;
40  for (std::size_t k = 0; k < A.shape(1); ++k)
41  for (std::size_t l = 0; l < B.shape(1); ++l)
42  acc += A(i, k) * B(k, l) * C(l, j);
43  r(i, j) = acc;
44  }
45 }
46 
47 template <typename Vec, typename Mat0, typename Mat1>
48 void dot21(Vec&& r, const Mat0& A, const Mat1& B)
49 {
50  using T = typename std::decay_t<Vec>::value_type;
51  for (std::size_t i = 0; i < r.shape(0); ++i)
52  {
53  T acc = 0;
54  for (std::size_t k = 0; k < A.shape(1); ++k)
55  acc += A(i, k) * B[k];
56  r[i] = acc;
57  }
58 }
59 } // namespace impl
60 
62 template <typename O, typename P, typename Q, typename R>
63 void l2_piola(O&& r, const P& U, const Q& /*J*/, double detJ, const R& /*K*/)
64 {
65  r.assign(U);
66  std::for_each(r.begin(), r.end(), [detJ](auto& ri) { ri /= detJ; });
67 }
68 
70 template <typename O, typename P, typename Q, typename R>
71 void covariant_piola(O&& r, const P& U, const Q& /*J*/, double /*detJ*/,
72  const R& K)
73 {
74  auto Kt = xt::transpose(K);
75  for (std::size_t p = 0; p < U.shape(0); ++p)
76  {
77  auto r_p = xt::row(r, p);
78  auto U_p = xt::row(U, p);
79  impl::dot21(r_p, Kt, U_p);
80  }
81 }
82 
84 template <typename O, typename P, typename Q, typename R>
85 void contravariant_piola(O&& r, const P& U, const Q& J, double detJ,
86  const R& /*K*/)
87 {
88  for (std::size_t p = 0; p < U.shape(0); ++p)
89  {
90  auto r_p = xt::row(r, p);
91  auto U_p = xt::row(U, p);
92  impl::dot21(r_p, J, U_p);
93  }
94  double inv_detJ = 1 / detJ;
95  std::for_each(r.begin(), r.end(), [inv_detJ](auto& ri) { ri *= inv_detJ; });
96 }
97 
99 template <typename O, typename P, typename Q, typename R>
100 void double_covariant_piola(O&& r, const P& U, const Q& J, double /*detJ*/,
101  const R& K)
102 {
103  for (std::size_t p = 0; p < U.shape(0); ++p)
104  {
105  auto r_p = xt::row(r, p);
106  auto U_p = xt::row(U, p);
107  auto _U = xt::reshape_view(U_p, {J.shape(1), J.shape(1)});
108  auto _r = xt::reshape_view(r_p, {K.shape(1), K.shape(1)});
109  impl::dot22(_r, xt::transpose(K), _U, K);
110  }
111 }
112 
114 template <typename O, typename P, typename Q, typename R>
115 void double_contravariant_piola(O&& r, const P& U, const Q& J, double detJ,
116  const R& /*K*/)
117 {
118  auto Jt = xt::transpose(J);
119  for (std::size_t p = 0; p < U.shape(0); ++p)
120  {
121  auto r_p = xt::row(r, p);
122  auto U_p = xt::row(U, p);
123  auto _U = xt::reshape_view(U_p, {J.shape(1), J.shape(1)});
124  auto _r = xt::reshape_view(r_p, {J.shape(0), J.shape(0)});
125  impl::dot22(_r, J, _U, Jt);
126  }
127  double inv_detJ = 1 / (detJ * detJ);
128  std::for_each(r.begin(), r.end(), [inv_detJ](auto& ri) { ri *= inv_detJ; });
129 }
130 
131 } // namespace basix::maps
basix::maps
Information about finite element maps.
Definition: maps.h:15
basix::maps::l2_piola
void l2_piola(O &&r, const P &U, const Q &, double detJ, const R &)
L2 Piola map.
Definition: maps.h:63
basix::maps::double_contravariant_piola
void double_contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Double contravariant Piola map.
Definition: maps.h:115
basix::maps::type
type
Map type.
Definition: maps.h:19
basix::maps::double_covariant_piola
void double_covariant_piola(O &&r, const P &U, const Q &J, double, const R &K)
Double covariant Piola map.
Definition: maps.h:100
basix::maps::contravariant_piola
void contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Contravariant Piola map.
Definition: maps.h:85
basix::maps::covariant_piola
void covariant_piola(O &&r, const P &U, const Q &, double, const R &K)
Covariant Piola map.
Definition: maps.h:71