Basix
maps.h
1 // Copyright (c) 2021 Matthew Scroggs
2 // FEniCS Project
3 // SPDX-License-Identifier: MIT
4 
5 #pragma once
6 
7 #include "span.hpp"
8 #include <map>
9 #include <stdexcept>
10 #include <string>
11 #include <vector>
12 #include <xtensor/xadapt.hpp>
13 #include <xtensor/xtensor.hpp>
14 #include <xtensor/xview.hpp>
15 
16 #include <xtensor/xio.hpp>
17 
19 namespace basix::maps
20 {
21 
23 enum class type
24 {
25  identity,
26  covariantPiola,
27  contravariantPiola,
28  doubleCovariantPiola,
29  doubleContravariantPiola,
30 };
31 
32 // Get the function that maps data from the reference to the physical
33 // cell
34 // @param mapping_type Mapping type
35 // @return The mapping function
36 // std::function<std::vector<double>(const tcb::span<const double>&,
37 // const xt::xtensor<double, 2>&, double,
38 // const xt::xtensor<double, 2>&)>
39 // get_forward_map(maps::type mapping_type);
40 
42 const std::string& type_to_str(maps::type type);
43 
44 namespace impl
45 {
46 template <typename O, typename Mat0, typename Mat1, typename Mat2>
47 void dot22(O& r, const Mat0& A, const Mat1& B, const Mat2& C)
48 {
49  assert(A.shape(1) == B.shape(0));
50  r = 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);
56 }
57 
58 template <typename Vec, typename Mat0, typename Mat1>
59 void dot21(Vec& r, const Mat0& A, const Mat1& B)
60 {
61  // assert(A.shape(1) == B.shape(0));
62  r = 0;
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];
66 }
67 
68 template <typename Vec0, typename Vec1, typename Mat0, typename Mat1>
69 void identity(Vec0& r, const Vec1& U, const Mat0& /*J*/, double /*detJ*/,
70  const Mat1& /*K*/)
71 {
72  r = U;
73 }
74 
75 template <typename O, typename P, typename Q, typename R>
76 void covariant_piola(O&& r, const P& U, const Q& /*J*/, double /*detJ*/,
77  const R& K)
78 {
79  dot21(r, xt::transpose(K), U);
80 }
81 
82 template <typename O, typename P, typename Q, typename R>
83 void contravariant_piola(O&& r, const P& U, const Q& J, double detJ,
84  const R& /*K*/)
85 {
86  dot21(r, J, U);
87  r /= detJ;
88 }
89 
90 template <typename O, typename P, typename Q, typename R>
91 void double_covariant_piola(O& r, const P& U, const Q& J, double /*detJ*/,
92  const R& K)
93 {
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);
97 }
98 
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,
101  const R& /*K*/)
102 {
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));
106  _r /= (detJ * detJ);
107 }
108 } // namespace impl
109 
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,
122  maps::type map_type)
123 {
124  switch (map_type)
125  {
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);
136  default:
137  throw std::runtime_error("Mapping not yet implemented");
138  }
139 }
140 
141 } // namespace basix::maps
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