Basix
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
12namespace basix::maps
13{
14
16enum class type
17{
18 identity = 0,
19 L2Piola = 1,
20 covariantPiola = 2,
21 contravariantPiola = 3,
22 doubleCovariantPiola = 4,
23 doubleContravariantPiola = 5,
24};
25
27template <typename O, typename P, typename Q, typename R>
28void 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
38template <typename O, typename P, typename Q, typename R>
39void 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
57template <typename O, typename P, typename Q, typename R>
58void 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
78template <typename O, typename P, typename Q, typename R>
79void 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
106template <typename O, typename P, typename Q, typename R>
107void 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