12 #include <xtensor/xadapt.hpp>
13 #include <xtensor/xtensor.hpp>
14 #include <xtensor/xview.hpp>
16 #include <xtensor/xio.hpp>
29 doubleContravariantPiola,
46 template <
typename O,
typename Mat0,
typename Mat1,
typename Mat2>
47 void dot22(O& r,
const Mat0& A,
const Mat1& B,
const Mat2& C)
49 assert(A.shape(1) == B.shape(0));
51 for (std::size_t i = 0; i < r.shape(0); ++i)
52 for (std::size_t j = 0; j < r.shape(1); ++j)
53 for (std::size_t k = 0; k < A.shape(1); ++k)
54 for (std::size_t l = 0; l < B.shape(1); ++l)
55 r(i, j) += A(i, k) * B(k, l) * C(l, j);
58 template <
typename Vec,
typename Mat0,
typename Mat1>
59 void dot21(Vec& r,
const Mat0& A,
const Mat1& B)
63 for (std::size_t i = 0; i < r.shape(0); ++i)
64 for (std::size_t k = 0; k < A.shape(1); ++k)
65 r[i] += A(i, k) * B[k];
68 template <
typename Vec0,
typename Vec1,
typename Mat0,
typename Mat1>
69 void identity(Vec0& r,
const Vec1& U,
const Mat0& ,
double ,
75 template <
typename O,
typename P,
typename Q,
typename R>
76 void covariant_piola(O&& r,
const P& U,
const Q& ,
double ,
79 dot21(r, xt::transpose(K), U);
82 template <
typename O,
typename P,
typename Q,
typename R>
83 void contravariant_piola(O&& r,
const P& U,
const Q& J,
double detJ,
90 template <
typename O,
typename P,
typename Q,
typename R>
91 void double_covariant_piola(O& r,
const P& U,
const Q& J,
double ,
94 auto _U = xt::reshape_view(U, {J.shape(1), J.shape(1)});
95 auto _r = xt::reshape_view(r, {K.shape(1), K.shape(1)});
96 dot22(_r, xt::transpose(K), _U, K);
99 template <
typename O,
typename P,
typename Q,
typename R>
100 void double_contravariant_piola(O& r,
const P& U,
const Q& J,
double detJ,
103 auto _U = xt::reshape_view(U, {J.shape(1), J.shape(1)});
104 auto _r = xt::reshape_view(r, {J.shape(0), J.shape(0)});
105 dot22(_r, J, _U, xt::transpose(J));
120 template <
typename O,
typename P,
typename Mat0,
typename Mat1>
121 void apply_map(O&& u,
const P& U,
const Mat0& J,
double detJ,
const Mat1& K,
126 case maps::type::identity:
127 return impl::identity(u, U, J, detJ, K);
128 case maps::type::covariantPiola:
129 return impl::covariant_piola(u, U, J, detJ, K);
130 case maps::type::contravariantPiola:
131 return impl::contravariant_piola(u, U, J, detJ, K);
132 case maps::type::doubleCovariantPiola:
133 return impl::double_covariant_piola(u, U, J, detJ, K);
134 case maps::type::doubleContravariantPiola:
135 return impl::double_contravariant_piola(u, U, J, detJ, K);
137 throw std::runtime_error(
"Mapping not yet implemented");
Information about finite element maps.
Definition: maps.h:20
void apply_map(O &&u, const P &U, const Mat0 &J, double detJ, const Mat1 &K, maps::type map_type)
Definition: maps.h:121
type
Cell type.
Definition: maps.h:24
const std::string & type_to_str(maps::type type)
Convert mapping type enum to string.
Definition: maps.cpp:10