8#include "element-families.h"
12#include "precompute.h"
13#include "sobolev-spaces.h"
32template <
typename T, std::
size_t d>
33using mdspan_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
34 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, d>>;
35template <
typename T, std::
size_t d>
37 = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE::mdarray<
38 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, d>>;
43std::array<std::vector<mdspan_t<const T, 2>>, 4>
44to_mdspan(std::array<std::vector<mdarray_t<T, 2>>, 4>& x)
46 std::array<std::vector<mdspan_t<const T, 2>>, 4> x1;
47 for (std::size_t i = 0; i < x.size(); ++i)
48 for (std::size_t j = 0; j < x[i].size(); ++j)
49 x1[i].emplace_back(x[i][j].data(), x[i][j].extents());
57std::array<std::vector<mdspan_t<const T, 4>>, 4>
58to_mdspan(std::array<std::vector<mdarray_t<T, 4>>, 4>& M)
60 std::array<std::vector<mdspan_t<const T, 4>>, 4> M1;
61 for (std::size_t i = 0; i < M.size(); ++i)
62 for (std::size_t j = 0; j < M[i].size(); ++j)
63 M1[i].emplace_back(M[i][j].data(), M[i][j].extents());
71std::array<std::vector<mdspan_t<const T, 2>>, 4>
72to_mdspan(
const std::array<std::vector<std::vector<T>>, 4>& x,
73 const std::array<std::vector<std::array<std::size_t, 2>>, 4>& shape)
75 std::array<std::vector<mdspan_t<const T, 2>>, 4> x1;
76 for (std::size_t i = 0; i < x.size(); ++i)
77 for (std::size_t j = 0; j < x[i].size(); ++j)
78 x1[i].push_back(mdspan_t<const T, 2>(x[i][j].data(), shape[i][j]));
86std::array<std::vector<mdspan_t<const T, 4>>, 4>
87to_mdspan(
const std::array<std::vector<std::vector<T>>, 4>& M,
88 const std::array<std::vector<std::array<std::size_t, 4>>, 4>& shape)
90 std::array<std::vector<mdspan_t<const T, 4>>, 4> M1;
91 for (std::size_t i = 0; i < M.size(); ++i)
92 for (std::size_t j = 0; j < M[i].size(); ++j)
93 M1[i].push_back(mdspan_t<const T, 4>(M[i][j].data(), shape[i][j]));
103template <
typename T, std::
size_t d>
121template <std::
floating_po
int T>
122std::tuple<std::array<std::vector<std::vector<T>>, 4>,
123 std::array<std::vector<std::array<std::size_t, 2>>, 4>,
124 std::array<std::vector<std::vector<T>>, 4>,
125 std::array<std::vector<std::array<std::size_t, 4>>, 4>>
137template <std::
floating_po
int F>
140 template <
typename T, std::
size_t d>
141 using mdspan_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
142 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, d>>;
329 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
368 for (std::size_t
i = 1;
i <=
nd; ++
i)
370 for (std::size_t
i = 1;
i <=
nd; ++
i)
372 std::size_t
vs = std::accumulate(_value_shape.begin(), _value_shape.end(),
373 1, std::multiplies{});
374 std::size_t
ndofs = _coeffs.second[0];
399 std::pair<std::vector<F>, std::array<std::size_t, 4>>
400 tabulate(
int nd, impl::mdspan_t<const F, 2>
x)
const;
425 std::pair<std::vector<F>, std::array<std::size_t, 4>>
427 std::array<std::size_t, 2>
shape)
const;
454 void tabulate(
int nd, impl::mdspan_t<const F, 2>
x,
483 std::span<F>
basis)
const;
511 const std::vector<std::size_t>&
value_shape()
const {
return _value_shape; }
516 int dim()
const {
return _coeffs.second[0]; }
526 return _lagrange_variant;
550 return _dof_transformations_are_permutations;
557 return _dof_transformations_are_identity;
573 std::pair<std::vector<F>, std::array<std::size_t, 3>>
574 push_forward(impl::mdspan_t<const F, 3>
U, impl::mdspan_t<const F, 3>
J,
575 std::span<const F>
detJ, impl::mdspan_t<const F, 3>
K)
const;
584 std::pair<std::vector<F>, std::array<std::size_t, 3>>
585 pull_back(impl::mdspan_t<const F, 3>
u, impl::mdspan_t<const F, 3>
J,
586 std::span<const F>
detJ, impl::mdspan_t<const F, 3>
K)
const;
619 template <
typename O,
typename P,
typename Q,
typename R>
624 case maps::type::identity:
625 return [](
O&
u,
const P&
U,
const Q&,
F,
const R&)
629 for (std::size_t
i = 0;
i <
U.extent(0); ++
i)
630 for (std::size_t
j = 0;
j <
U.extent(1); ++
j)
633 case maps::type::covariantPiola:
634 return [](
O&
u,
const P&
U,
const Q&
J,
F detJ,
const R&
K)
636 case maps::type::contravariantPiola:
637 return [](
O&
u,
const P&
U,
const Q&
J,
F detJ,
const R&
K)
639 case maps::type::doubleCovariantPiola:
640 return [](
O&
u,
const P&
U,
const Q&
J,
F detJ,
const R&
K)
642 case maps::type::doubleContravariantPiola:
643 return [](
O&
u,
const P&
U,
const Q&
J,
F detJ,
const R&
K)
646 throw std::runtime_error(
"Map not implemented");
656 const std::vector<std::vector<std::vector<int>>>&
entity_dofs()
const
670 return _e_closure_dofs;
754 std::pair<std::vector<F>, std::array<std::size_t, 3>>
761 std::map<cell::type, std::pair<std::vector<F>, std::array<std::size_t, 3>>>
764 return _entity_transformations;
776 if (!_dof_transformations_are_permutations)
778 throw std::runtime_error(
779 "The DOF transformations for this element are not permutations");
782 if (_dof_transformations_are_identity)
798 if (!_dof_transformations_are_permutations)
800 throw std::runtime_error(
801 "The DOF transformations for this element are not permutations");
803 if (_dof_transformations_are_identity)
817 template <
typename T>
829 template <
typename T>
841 template <
typename T>
854 template <
typename T>
866 template <
typename T>
878 template <
typename T>
891 template <
typename T>
903 template <
typename T>
911 const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
points()
const
968 const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
979 const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
1019 const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
wcoeffs()
const
1029 std::vector<std::pair<std::vector<F>, std::array<std::size_t, 2>>>, 4>&
1072 std::vector<std::pair<std::vector<F>, std::array<std::size_t, 4>>>, 4>&
1083 const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
1098 return _tensor_factors.size() > 0;
1112 std::vector<std::tuple<std::vector<FiniteElement<F>>, std::vector<int>>>
1116 throw std::runtime_error(
"Element has no tensor product representation.");
1117 return _tensor_factors;
1137 template <
typename T,
bool post>
1140 const std::map<
cell::type, std::vector<std::vector<std::size_t>>>&
eperm)
1143 using array2_t = std::pair<std::vector<F>, std::array<std::size_t, 2>>;
1144 using array3_t = std::pair<std::vector<F>, std::array<std::size_t, 3>>;
1146 = std::vector<std::pair<std::vector<std::size_t>, array2_t>>;
1154 template <
typename T,
bool post,
typename OP>
1157 const std::map<cell::type, trans_data_t>&
etrans,
OP op)
const;
1166 std::size_t _cell_tdim;
1169 std::vector<std::vector<cell::type>> _cell_subentity_types;
1184 int _interpolation_nderivs;
1187 int _highest_degree;
1190 int _highest_complete_degree;
1193 std::vector<std::size_t> _value_shape;
1206 std::pair<std::vector<F>, std::array<std::size_t, 2>> _coeffs;
1209 std::vector<std::vector<std::vector<int>>> _edofs;
1212 std::vector<std::vector<std::vector<int>>> _e_closure_dofs;
1215 std::map<cell::type, array3_t> _entity_transformations;
1222 std::pair<std::vector<F>, std::array<std::size_t, 2>> _points;
1226 std::array<std::vector<std::pair<std::vector<F>, std::array<std::size_t, 2>>>,
1231 std::pair<std::vector<F>, std::array<std::size_t, 2>> _matM;
1235 bool _dof_transformations_are_permutations;
1238 bool _dof_transformations_are_identity;
1243 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm;
1248 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_rev;
1251 std::map<cell::type, trans_data_t> _etrans;
1254 std::map<cell::type, trans_data_t> _etransT;
1257 std::map<cell::type, trans_data_t> _etrans_inv;
1260 std::map<cell::type, trans_data_t> _etrans_invT;
1264 bool _discontinuous;
1267 std::pair<std::vector<F>, std::array<std::size_t, 2>> _dual_matrix;
1274 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
1283 std::vector<int> _dof_ordering;
1286 bool _interpolation_is_identity;
1290 std::pair<std::vector<F>, std::array<std::size_t, 2>> _wcoeffs;
1294 = std::vector<std::pair<std::vector<F>, std::array<std::size_t, 4>>>;
1295 std::array<array4_t, 4> _M;
1326template <std::
floating_po
int T>
1328 cell::type cell_type,
const std::vector<std::size_t>& value_shape,
1329 impl::mdspan_t<const T, 2> wcoeffs,
1330 const std::array<std::vector<impl::mdspan_t<const T, 2>>, 4>& x,
1331 const std::array<std::vector<impl::mdspan_t<const T, 4>>, 4>& M,
1332 int interpolation_nderivs,
maps::type map_type,
1334 int highest_complete_degree,
int highest_degree,
polyset::type poly_type);
1347template <std::
floating_po
int T>
1352 std::vector<int> dof_ordering = {});
1359template <std::
floating_po
int F>
1360template <
typename T,
bool post>
1361void FiniteElement<F>::permute_data(
1362 std::span<T> data,
int block_size, std::uint32_t cell_info,
1363 const std::map<
cell::type, std::vector<std::vector<std::size_t>>>& eperm)
1366 if (_cell_tdim >= 2)
1370 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1374 auto& trans = eperm.at(cell::type::interval)[0];
1375 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1378 if (cell_info >> (face_start + e) & 1)
1384 if (_cell_tdim == 3)
1387 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1389 auto& trans = eperm.at(_cell_subentity_types[2][f]);
1392 if (!post and cell_info >> (3 * f) & 1)
1397 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1402 if (post and cell_info >> (3 * f) & 1)
1410template <std::
floating_po
int F>
1411template <
typename T,
bool post,
typename OP>
1412void FiniteElement<F>::transform_data(
1413 std::span<T> data,
int block_size, std::uint32_t cell_info,
1414 const std::map<cell::type, trans_data_t>& etrans, OP op)
const
1416 if (_cell_tdim >= 2)
1420 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1422 for (
auto& edofs0 : _edofs[0])
1423 dofstart += edofs0.size();
1427 auto& [v_size_t, matrix] = etrans.at(cell::type::interval)[0];
1428 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1431 if (cell_info >> (face_start + e) & 1)
1433 op(std::span(v_size_t),
1434 mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
1435 dofstart, block_size);
1437 dofstart += _edofs[1][e].size();
1441 if (_cell_tdim == 3)
1444 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1446 auto& trans = etrans.at(_cell_subentity_types[2][f]);
1449 if (!post and cell_info >> (3 * f) & 1)
1451 const auto& m = trans[1];
1452 const auto& v_size_t = std::get<0>(m);
1453 const auto& matrix = std::get<1>(m);
1454 op(std::span(v_size_t),
1455 mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
1456 dofstart, block_size);
1460 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1462 const auto& m = trans[0];
1463 const auto& v_size_t = std::get<0>(m);
1464 const auto& matrix = std::get<1>(m);
1465 op(std::span(v_size_t),
1466 mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
1467 dofstart, block_size);
1471 if (post and cell_info >> (3 * f) & 1)
1473 const auto& m = trans[1];
1474 const auto& v_size_t = std::get<0>(m);
1475 const auto& matrix = std::get<1>(m);
1476 op(std::span(v_size_t),
1477 mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
1478 dofstart, block_size);
1481 dofstart += _edofs[2][f].size();
1487template <std::
floating_po
int F>
1488template <
typename T>
1493 if (_dof_transformations_are_identity)
1496 if (_dof_transformations_are_permutations)
1500 precompute::apply_matrix<F, T>);
1503template <std::
floating_po
int F>
1504template <
typename T>
1508 if (_dof_transformations_are_identity)
1511 if (_dof_transformations_are_permutations)
1515 precompute::apply_matrix<F, T>);
1518template <std::
floating_po
int F>
1519template <
typename T>
1523 if (_dof_transformations_are_identity)
1526 if (_dof_transformations_are_permutations)
1530 precompute::apply_matrix<F, T>);
1533template <std::
floating_po
int F>
1534template <
typename T>
1538 if (_dof_transformations_are_identity)
1541 if (_dof_transformations_are_permutations)
1545 precompute::apply_matrix<F, T>);
1548template <std::
floating_po
int F>
1549template <
typename T>
1553 if (_dof_transformations_are_identity)
1556 if (_dof_transformations_are_permutations)
1568 precompute::apply_matrix_to_transpose<F, T>);
1571template <std::
floating_po
int F>
1572template <
typename T>
1576 if (_dof_transformations_are_identity)
1579 if (_dof_transformations_are_permutations)
1591 precompute::apply_matrix_to_transpose<F, T>);
1594template <std::
floating_po
int F>
1595template <
typename T>
1599 if (_dof_transformations_are_identity)
1602 if (_dof_transformations_are_permutations)
1614 precompute::apply_matrix_to_transpose<F, T>);
1617template <std::
floating_po
int F>
1618template <
typename T>
1622 if (_dof_transformations_are_identity)
1625 if (_dof_transformations_are_permutations)
1638 precompute::apply_matrix_to_transpose<F, T>);
A finite element.
Definition finite-element.h:139
const std::array< std::vector< std::pair< std::vector< F >, std::array< std::size_t, 4 > > >, 4 > & M() const
Definition finite-element.h:1073
void apply_inverse_transpose_dof_transformation_to_transpose(std::span< T > data, int block_size, std::uint32_t cell_info) const
Apply inverse transpose DOF transformations to some transposed data.
Definition finite-element.h:1573
void apply_inverse_transpose_dof_transformation(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1520
std::pair< std::vector< F >, std::array< std::size_t, 3 > > base_transformations() const
Get the base transformations.
Definition finite-element.cpp:1088
int highest_complete_degree() const
Definition finite-element.h:506
bool dof_transformations_are_identity() const
Definition finite-element.h:555
bool operator==(const FiniteElement &e) const
Definition finite-element.cpp:959
std::map< cell::type, std::pair< std::vector< F >, std::array< std::size_t, 3 > > > entity_transformations() const
Definition finite-element.h:762
const std::vector< int > & dof_ordering() const
Get dof layout.
Definition finite-element.h:1128
const std::vector< std::vector< std::vector< int > > > & entity_dofs() const
Definition finite-element.h:656
FiniteElement(FiniteElement &&element)=default
Move constructor.
bool has_tensor_product_factorisation() const
Definition finite-element.h:1096
int interpolation_nderivs() const
The number of derivatives needed when interpolating.
Definition finite-element.h:1125
std::pair< std::vector< F >, std::array< std::size_t, 4 > > tabulate(int nd, impl::mdspan_t< const F, 2 > x) const
Compute basis values and derivatives at set of points.
Definition finite-element.cpp:997
int dim() const
Definition finite-element.h:516
void apply_inverse_dof_transformation(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1535
const std::vector< std::size_t > & value_shape() const
Definition finite-element.h:511
element::lagrange_variant lagrange_variant() const
Definition finite-element.h:524
const std::array< std::vector< std::pair< std::vector< F >, std::array< std::size_t, 2 > > >, 4 > & x() const
Definition finite-element.h:1030
int degree() const
Definition finite-element.h:495
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & points() const
Definition finite-element.h:911
sobolev::space sobolev_space() const
Definition finite-element.h:539
const std::vector< std::vector< std::vector< int > > > & entity_closure_dofs() const
Definition finite-element.h:668
FiniteElement(const FiniteElement &element)=default
Copy constructor.
std::vector< std::tuple< std::vector< FiniteElement< F > >, std::vector< int > > > get_tensor_product_representation() const
Definition finite-element.h:1113
std::array< std::size_t, 4 > tabulate_shape(std::size_t nd, std::size_t num_points) const
Definition finite-element.h:364
void apply_dof_transformation_to_transpose(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1550
FiniteElement & operator=(FiniteElement &&element)=default
Move assignment operator.
void apply_inverse_dof_transformation_to_transpose(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1619
void apply_dof_transformation(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1489
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & interpolation_matrix() const
Return a matrix of weights interpolation,.
Definition finite-element.h:969
polyset::type polyset_type() const
Definition finite-element.h:491
void permute_dofs(std::span< std::int32_t > dofs, std::uint32_t cell_info) const
Definition finite-element.h:774
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & wcoeffs() const
Definition finite-element.h:1019
std::pair< std::vector< F >, std::array< std::size_t, 3 > > pull_back(impl::mdspan_t< const F, 3 > u, impl::mdspan_t< const F, 3 > J, std::span< const F > detJ, impl::mdspan_t< const F, 3 > K) const
Definition finite-element.cpp:1197
bool dof_transformations_are_permutations() const
Definition finite-element.h:548
std::pair< std::vector< F >, std::array< std::size_t, 3 > > push_forward(impl::mdspan_t< const F, 3 > U, impl::mdspan_t< const F, 3 > J, std::span< const F > detJ, impl::mdspan_t< const F, 3 > K) const
Definition finite-element.cpp:1157
cell::type cell_type() const
Definition finite-element.h:487
bool interpolation_is_identity() const
Definition finite-element.h:1122
bool discontinuous() const
Definition finite-element.h:544
FiniteElement(element::family family, cell::type cell_type, polyset::type poly_type, int degree, const std::vector< std::size_t > &value_shape, mdspan_t< const F, 2 > wcoeffs, const std::array< std::vector< mdspan_t< const F, 2 > >, 4 > &x, const std::array< std::vector< mdspan_t< const F, 4 > >, 4 > &M, int interpolation_nderivs, maps::type map_type, sobolev::space sobolev_space, bool discontinuous, int highest_complete_degree, int highest_degree, element::lagrange_variant lvariant, element::dpc_variant dvariant, std::vector< std::tuple< std::vector< FiniteElement >, std::vector< int > > > tensor_factors={}, std::vector< int > dof_ordering={})
Construct a finite element.
maps::type map_type() const
Definition finite-element.h:535
void apply_transpose_dof_transformation_to_transpose(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1596
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & dual_matrix() const
Definition finite-element.h:980
element::dpc_variant dpc_variant() const
Definition finite-element.h:531
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & coefficient_matrix() const
Definition finite-element.h:1084
FiniteElement & operator=(const FiniteElement &element)=default
Assignment operator.
int highest_degree() const
Definition finite-element.h:501
void unpermute_dofs(std::span< std::int32_t > dofs, std::uint32_t cell_info) const
Definition finite-element.h:795
element::family family() const
Definition finite-element.h:520
void apply_transpose_dof_transformation(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1505
~FiniteElement()=default
Destructor.
std::function< void(O &, const P &, const Q &, F, const R &)> map_fn() const
Definition finite-element.h:620
type
Cell type.
Definition cell.h:21
std::tuple< std::array< std::vector< std::vector< T > >, 4 >, std::array< std::vector< std::array< std::size_t, 2 > >, 4 >, std::array< std::vector< std::vector< T > >, 4 >, std::array< std::vector< std::array< std::size_t, 4 > >, 4 > > make_discontinuous(const std::array< std::vector< mdspan_t< const T, 2 > >, 4 > &x, const std::array< std::vector< mdspan_t< const T, 4 > >, 4 > &M, std::size_t tdim, std::size_t value_size)
Definition finite-element.cpp:329
lagrange_variant
Variants of a Lagrange space that can be created.
Definition element-families.h:12
impl::mdspan_t< T, d > mdspan_t
Typedef for mdspan.
Definition finite-element.h:104
dpc_variant
Definition element-families.h:33
family
Available element families.
Definition element-families.h:46
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
type
Cell type.
Definition polyset.h:136
void apply_permutation_mapped(std::span< const std::size_t > perm, std::span< E > data, std::span< const int > emap, std::size_t block_size=1)
Permutation of mapped data.
Definition precompute.h:156
space
Sobolev space type.
Definition sobolev-spaces.h:13
Basix: FEniCS runtime basis evaluation library.
Definition cell.h:17
FiniteElement< T > create_element(element::family family, cell::type cell, int degree, element::lagrange_variant lvariant, element::dpc_variant dvariant, bool discontinuous, std::vector< int > dof_ordering={})
Definition finite-element.cpp:189
std::string version()
Definition finite-element.cpp:1235
FiniteElement< T > create_custom_element(cell::type cell_type, const std::vector< std::size_t > &value_shape, impl::mdspan_t< const T, 2 > wcoeffs, const std::array< std::vector< impl::mdspan_t< const T, 2 > >, 4 > &x, const std::array< std::vector< impl::mdspan_t< const T, 4 > >, 4 > &M, int interpolation_nderivs, maps::type map_type, sobolev::space sobolev_space, bool discontinuous, int highest_complete_degree, int highest_degree, polyset::type poly_type)
Definition finite-element.cpp:418