Basix 0.8.0

HomeInstallationDemosC++ docsPython docs

Searching...
No Matches
maps.h
1// Copyright (c) 2021-2022 Matthew Scroggs and Garth N. Wells
2// FEniCS Project
4
5#pragma once
6
7#include "mdspan.hpp"
8#include <stdexcept>
9#include <type_traits>
10
12namespace basix::maps
13{
14
15namespace impl
16{
19template <typename T, typename = void>
20struct scalar_value_type
21{
23 typedef T value_type;
24};
26template <typename T>
27struct scalar_value_type<T, std::void_t<typename T::value_type>>
28{
29 typedef typename T::value_type value_type;
30};
32template <typename T>
33using scalar_value_type_t = typename scalar_value_type<T>::value_type;
34} // namespace impl
35
37enum class type
38{
39 identity = 0,
40 L2Piola = 1,
41 covariantPiola = 2,
42 contravariantPiola = 3,
43 doubleCovariantPiola = 4,
44 doubleContravariantPiola = 5,
45};
46
48template <typename O, typename P, typename Q, typename R>
49void l2_piola(O&& r, const P& U, const Q& /*J*/, double detJ, const R& /*K*/)
50{
51 assert(U.extent(0) == r.extent(0));
52 assert(U.extent(1) == r.extent(1));
53 for (std::size_t i = 0; i < U.extent(0); ++i)
54 for (std::size_t j = 0; j < U.extent(1); ++j)
55 r(i, j) = U(i, j) / detJ;
56}
57
59template <typename O, typename P, typename Q, typename R>
60void covariant_piola(O&& r, const P& U, const Q& /*J*/, double /*detJ*/,
61 const R& K)
62{
63 using T = typename std::decay_t<O>::value_type;
64 using Z = typename impl::scalar_value_type_t<T>;
65 for (std::size_t p = 0; p < U.extent(0); ++p)
66 {
67 // r_p = K^T U_p, where p indicates the p-th row
68 for (std::size_t i = 0; i < r.extent(1); ++i)
69 {
70 T acc = 0;
71 for (std::size_t k = 0; k < K.extent(0); ++k)
72 acc += static_cast<Z>(K(k, i)) * U(p, k);
73 r(p, i) = acc;
74 }
75 }
76}
77
79template <typename O, typename P, typename Q, typename R>
80void contravariant_piola(O&& r, const P& U, const Q& J, double detJ,
81 const R& /*K*/)
82{
83 using T = typename std::decay_t<O>::value_type;
84 using Z = typename impl::scalar_value_type_t<T>;
85 for (std::size_t p = 0; p < U.extent(0); ++p)
86 {
87 for (std::size_t i = 0; i < r.extent(1); ++i)
88 {
89 T acc = 0;
90 for (std::size_t k = 0; k < J.extent(1); ++k)
91 acc += static_cast<Z>(J(i, k)) * U(p, k);
92 r(p, i) = acc;
93 }
94 }
95
96 std::transform(r.data_handle(), r.data_handle() + r.size(), r.data_handle(),
97 [detJ](auto ri) { return ri / static_cast<Z>(detJ); });
98}
99
101template <typename O, typename P, typename Q, typename R>
102void double_covariant_piola(O&& r, const P& U, const Q& J, double /*detJ*/,
103 const R& K)
104{
105 namespace stdex
106 = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE;
107 using T = typename std::decay_t<O>::value_type;
108 using Z = typename impl::scalar_value_type_t<T>;
109 for (std::size_t p = 0; p < U.extent(0); ++p)
110 {
111 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
112 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
113 _U(U.data_handle() + p * U.extent(1), J.extent(1), J.extent(1));
114 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
115 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
116 _r(r.data_handle() + p * r.extent(1), K.extent(1), K.extent(1));
117 // _r = K^T _U K
118 for (std::size_t i = 0; i < _r.extent(0); ++i)
119 {
120 for (std::size_t j = 0; j < _r.extent(1); ++j)
121 {
122 T acc = 0;
123 for (std::size_t k = 0; k < K.extent(0); ++k)
124 for (std::size_t l = 0; l < _U.extent(1); ++l)
125 acc += static_cast<Z>(K(k, i)) * _U(k, l) * static_cast<Z>(K(l, j));
126 _r(i, j) = acc;
127 }
128 }
129 }
130}
131
133template <typename O, typename P, typename Q, typename R>
134void double_contravariant_piola(O&& r, const P& U, const Q& J, double detJ,
135 const R& /*K*/)
136{
137 namespace stdex
138 = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE;
139 using T = typename std::decay_t<O>::value_type;
140 using Z = typename impl::scalar_value_type_t<T>;
141 for (std::size_t p = 0; p < U.extent(0); ++p)
142 {
143 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
144 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
145 _U(U.data_handle() + p * U.extent(1), J.extent(1), J.extent(1));
146 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
147 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
148 _r(r.data_handle() + p * r.extent(1), J.extent(0), J.extent(0));
149
150 // _r = J U J^T
151 for (std::size_t i = 0; i < _r.extent(0); ++i)
152 {
153 for (std::size_t j = 0; j < _r.extent(1); ++j)
154 {
155 T acc = 0;
156 for (std::size_t k = 0; k < J.extent(1); ++k)
157 for (std::size_t l = 0; l < _U.extent(1); ++l)
158 acc += static_cast<Z>(J(i, k)) * _U(k, l) * static_cast<Z>(J(j, l));
159 _r(i, j) = acc;
160 }
161 }
162 }
163
164 std::transform(r.data_handle(), r.data_handle() + r.size(), r.data_handle(),
165 [detJ](auto ri) { return ri / static_cast<Z>(detJ * detJ); });
166}
167
168} // namespace basix::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:49
void covariant_piola(O &&r, const P &U, const Q &, double, const R &K)
Covariant Piola map.
Definition: maps.h:60
void contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Contravariant Piola map.
Definition: maps.h:80
void double_contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Double contravariant Piola map.
Definition: maps.h:134
void double_covariant_piola(O &&r, const P &U, const Q &J, double, const R &K)
Double covariant Piola map.
Definition: maps.h:102
type
Map type.
Definition: maps.h:38