8 #include "element-families.h"
10 #include "precompute.h"
16 #include <xtensor/xtensor.hpp>
17 #include <xtensor/xview.hpp>
18 #include <xtl/xspan.hpp>
39 std::tuple<std::array<std::vector<xt::xtensor<double, 2>>, 4>,
40 std::array<std::vector<xt::xtensor<double, 3>>, 4>>
42 const std::array<std::vector<xt::xtensor<double, 3>>, 4>& M,
43 int tdim,
int value_size);
222 const xt::xtensor<double, 2>&
wcoeffs,
223 const std::array<std::vector<xt::xtensor<double, 2>>, 4>&
x,
224 const std::array<std::vector<xt::xtensor<double, 3>>, 4>&
M,
226 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
263 std::size_t num_points)
const;
286 xt::xtensor<double, 4>
tabulate(
int nd,
const xt::xarray<double>&
x)
const;
313 void tabulate(
int nd,
const xt::xarray<double>&
x,
314 xt::xtensor<double, 4>& basis)
const;
376 xt::xtensor<double, 3>
push_forward(
const xt::xtensor<double, 3>& U,
377 const xt::xtensor<double, 3>& J,
378 const xtl::span<const double>& detJ,
379 const xt::xtensor<double, 3>& K)
const;
387 xt::xtensor<double, 3>
pull_back(
const xt::xtensor<double, 3>& u,
388 const xt::xtensor<double, 3>& J,
389 const xtl::span<const double>& detJ,
390 const xt::xtensor<double, 3>& K)
const;
423 template <
typename O,
typename P,
typename Q,
typename R>
424 std::function<void(O&,
const P&,
const Q&,
double,
const R&)>
map_fn()
const
428 case maps::type::identity:
429 return [](O& u,
const P& U,
const Q&, double,
const R&) { u.assign(U); };
430 case maps::type::L2Piola:
431 return [](O& u,
const P& U,
const Q& J,
double detJ,
const R& K) {
434 case maps::type::covariantPiola:
435 return [](O& u,
const P& U,
const Q& J,
double detJ,
const R& K) {
438 case maps::type::contravariantPiola:
439 return [](O& u,
const P& U,
const Q& J,
double detJ,
const R& K) {
442 case maps::type::doubleCovariantPiola:
443 return [](O& u,
const P& U,
const Q& J,
double detJ,
const R& K) {
446 case maps::type::doubleContravariantPiola:
447 return [](O& u,
const P& U,
const Q& J,
double detJ,
const R& K) {
451 throw std::runtime_error(
"Map not implemented");
482 const std::vector<std::vector<std::vector<int>>>&
entity_dofs()
const;
585 std::uint32_t cell_info)
const;
595 std::uint32_t cell_info)
const;
605 template <
typename T>
607 std::uint32_t cell_info)
const;
617 template <
typename T>
620 std::uint32_t cell_info)
const;
630 template <
typename T>
632 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const;
642 template <
typename T>
645 std::uint32_t cell_info)
const;
655 template <
typename T>
658 std::uint32_t cell_info)
const;
668 template <
typename T>
670 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const;
680 template <
typename T>
682 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const;
692 template <
typename T>
694 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const;
700 const xt::xtensor<double, 2>&
points()
const;
791 const xt::xtensor<double, 2>&
wcoeffs()
const;
796 const std::array<std::vector<xt::xtensor<double, 2>>, 4>&
x()
const;
830 const std::array<std::vector<xt::xtensor<double, 3>>, 4>&
M()
const;
861 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
873 std::size_t _cell_tdim;
876 std::vector<std::vector<cell::type>> _cell_subentity_types;
891 std::vector<int> _value_shape;
901 xt::xtensor<double, 2> _coeffs;
909 std::vector<std::vector<int>> _num_edofs;
913 std::vector<std::vector<int>> _num_e_closure_dofs;
916 std::vector<std::vector<std::vector<int>>> _edofs;
919 std::vector<std::vector<std::vector<int>>> _e_closure_dofs;
922 std::map<cell::type, xt::xtensor<double, 3>> _entity_transformations;
929 xt::xtensor<double, 2> _points;
933 std::array<std::vector<xt::xtensor<double, 2>>, 4> _x;
936 xt::xtensor<double, 2> _matM;
940 bool _dof_transformations_are_permutations;
943 bool _dof_transformations_are_identity;
948 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm;
953 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_rev;
957 std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
958 xt::xtensor<double, 2>>>>
963 std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
964 xt::xtensor<double, 2>>>>
969 std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
970 xt::xtensor<double, 2>>>>
975 std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
976 xt::xtensor<double, 2>>>>
984 xt::xtensor<double, 2> _dual_matrix;
990 std::array<int, 2> _degree_bounds;
997 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
1001 bool _interpolation_is_identity;
1005 xt::xtensor<double, 2> _wcoeffs;
1008 std::array<std::vector<xt::xtensor<double, 3>>, 4> _M;
1032 const std::vector<std::size_t>& value_shape,
1033 const xt::xtensor<double, 2>& wcoeffs,
1034 const std::array<std::vector<xt::xtensor<double, 2>>, 4>& x,
1035 const std::array<std::vector<xt::xtensor<double, 3>>, 4>& M,
1036 maps::type map_type,
bool discontinuous,
int highest_complete_degree);
1049 bool discontinuous);
1076 bool discontinuous);
1089 int degree,
bool discontinuous);
1138 template <
typename T>
1141 std::uint32_t cell_info)
const
1143 if (_dof_transformations_are_identity)
1146 if (_cell_tdim >= 2)
1150 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1152 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1155 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1158 if (cell_info >> (face_start + e) & 1)
1160 dofstart, block_size);
1161 dofstart += _num_edofs[1][e];
1164 if (_cell_tdim == 3)
1167 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1170 if (cell_info >> (3 * f) & 1)
1172 data, dofstart, block_size);
1175 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1177 data, dofstart, block_size);
1178 dofstart += _num_edofs[2][f];
1184 template <
typename T>
1186 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1188 if (_dof_transformations_are_identity)
1191 if (_cell_tdim >= 2)
1195 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1197 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1200 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1203 if (cell_info >> (face_start + e) & 1)
1205 dofstart, block_size);
1206 dofstart += _num_edofs[1][e];
1209 if (_cell_tdim == 3)
1212 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1215 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1217 data, dofstart, block_size);
1219 if (cell_info >> (3 * f) & 1)
1221 data, dofstart, block_size);
1222 dofstart += _num_edofs[2][f];
1228 template <
typename T>
1230 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1232 if (_dof_transformations_are_identity)
1235 if (_cell_tdim >= 2)
1239 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1241 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1244 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1247 if (cell_info >> (face_start + e) & 1)
1249 dofstart, block_size);
1250 dofstart += _num_edofs[1][e];
1253 if (_cell_tdim == 3)
1256 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1259 if (cell_info >> (3 * f) & 1)
1261 _etrans_invT.at(_cell_subentity_types[2][f])[1], data, dofstart,
1265 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1267 _etrans_invT.at(_cell_subentity_types[2][f])[0], data, dofstart,
1269 dofstart += _num_edofs[2][f];
1275 template <
typename T>
1277 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1279 if (_dof_transformations_are_identity)
1282 if (_cell_tdim >= 2)
1286 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1288 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1291 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1294 if (cell_info >> (face_start + e) & 1)
1296 dofstart, block_size);
1297 dofstart += _num_edofs[1][e];
1300 if (_cell_tdim == 3)
1303 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1306 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1308 _etrans_inv.at(_cell_subentity_types[2][f])[0], data, dofstart,
1311 if (cell_info >> (3 * f) & 1)
1313 _etrans_inv.at(_cell_subentity_types[2][f])[1], data, dofstart,
1315 dofstart += _num_edofs[2][f];
1321 template <
typename T>
1323 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1325 if (_dof_transformations_are_identity)
1328 if (_cell_tdim >= 2)
1332 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1334 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1337 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1340 if (cell_info >> (face_start + e) & 1)
1342 _etrans.at(cell::type::interval)[0], data, dofstart, block_size);
1343 dofstart += _num_edofs[1][e];
1346 if (_cell_tdim == 3)
1349 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1352 if (cell_info >> (3 * f) & 1)
1354 _etrans.at(_cell_subentity_types[2][f])[1], data, dofstart,
1358 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1360 _etrans.at(_cell_subentity_types[2][f])[0], data, dofstart,
1362 dofstart += _num_edofs[2][f];
1368 template <
typename T>
1370 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1372 if (_dof_transformations_are_identity)
1375 if (_cell_tdim >= 2)
1379 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1381 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1384 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1387 if (cell_info >> (face_start + e) & 1)
1389 _etrans_invT.at(cell::type::interval)[0], data, dofstart,
1391 dofstart += _num_edofs[1][e];
1394 if (_cell_tdim == 3)
1397 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1400 if (cell_info >> (3 * f) & 1)
1402 _etrans_invT.at(_cell_subentity_types[2][f])[1], data, dofstart,
1406 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1408 _etrans_invT.at(_cell_subentity_types[2][f])[0], data, dofstart,
1410 dofstart += _num_edofs[2][f];
1416 template <
typename T>
1418 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1420 if (_dof_transformations_are_identity)
1423 if (_cell_tdim >= 2)
1427 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1429 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1432 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1435 if (cell_info >> (face_start + e) & 1)
1437 _etransT.at(cell::type::interval)[0], data, dofstart, block_size);
1438 dofstart += _num_edofs[1][e];
1441 if (_cell_tdim == 3)
1444 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1447 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1449 _etransT.at(_cell_subentity_types[2][f])[0], data, dofstart,
1453 if (cell_info >> (3 * f) & 1)
1455 _etransT.at(_cell_subentity_types[2][f])[1], data, dofstart,
1458 dofstart += _num_edofs[2][f];
1464 template <
typename T>
1466 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1468 if (_dof_transformations_are_identity)
1471 if (_cell_tdim >= 2)
1475 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1477 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1480 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1483 if (cell_info >> (face_start + e) & 1)
1485 _etrans_inv.at(cell::type::interval)[0], data, dofstart,
1487 dofstart += _num_edofs[1][e];
1490 if (_cell_tdim == 3)
1493 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1496 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1498 _etrans_inv.at(_cell_subentity_types[2][f])[0], data, dofstart,
1502 if (cell_info >> (3 * f) & 1)
1504 _etrans_inv.at(_cell_subentity_types[2][f])[1], data, dofstart,
1507 dofstart += _num_edofs[2][f];
A finite element.
Definition: finite-element.h:53
FiniteElement & operator=(const FiniteElement &element)=default
Assignment operator.
cell::type cell_type() const
Definition: finite-element.cpp:805
const std::array< std::vector< xt::xtensor< double, 2 > >, 4 > & x() const
Definition: finite-element.cpp:1110
FiniteElement & operator=(FiniteElement &&element)=default
Move assignment operator.
xt::xtensor< double, 4 > tabulate(int nd, const xt::xarray< double > &x) const
Definition: finite-element.cpp:761
std::array< int, 2 > degree_bounds() const
Definition: finite-element.cpp:1128
xt::xtensor< double, 3 > push_forward(const xt::xtensor< double, 3 > &U, const xt::xtensor< double, 3 > &J, const xtl::span< const double > &detJ, const xt::xtensor< double, 3 > &K) const
Definition: finite-element.cpp:923
xt::xtensor< double, 3 > pull_back(const xt::xtensor< double, 3 > &u, const xt::xtensor< double, 3 > &J, const xtl::span< const double > &detJ, const xt::xtensor< double, 3 > &K) const
Definition: finite-element.cpp:951
element::lagrange_variant lagrange_variant() const
Definition: finite-element.cpp:1138
FiniteElement(FiniteElement &&element)=default
Move constructor.
const xt::xtensor< double, 2 > & dual_matrix() const
Definition: finite-element.cpp:1097
std::map< cell::type, xt::xtensor< double, 3 > > entity_transformations() const
Definition: finite-element.cpp:1092
bool operator==(const FiniteElement &e) const
Definition: finite-element.cpp:738
bool dof_transformations_are_permutations() const
Definition: finite-element.cpp:822
maps::type map_type() const
Definition: finite-element.cpp:818
const xt::xtensor< double, 2 > & points() const
Definition: finite-element.cpp:921
void apply_inverse_dof_transformation_to_transpose(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1465
std::vector< std::tuple< std::vector< FiniteElement >, std::vector< int > > > get_tensor_product_representation() const
Definition: finite-element.cpp:1151
const std::vector< std::vector< std::vector< int > > > & entity_dofs() const
Definition: finite-element.cpp:843
const std::array< std::vector< xt::xtensor< double, 3 > >, 4 > & M() const
Definition: finite-element.cpp:1116
const xt::xtensor< double, 2 > & interpolation_matrix() const
Definition: finite-element.cpp:832
void apply_transpose_dof_transformation(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1185
const std::vector< std::vector< std::vector< int > > > & entity_closure_dofs() const
Definition: finite-element.cpp:855
const std::vector< std::vector< int > > & num_entity_dofs() const
Definition: finite-element.cpp:837
int degree() const
Definition: finite-element.cpp:807
void apply_inverse_transpose_dof_transformation_to_transpose(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1369
std::function< void(O &, const P &, const Q &, double, const R &)> map_fn() const
Definition: finite-element.h:424
const std::vector< int > & value_shape() const
Definition: finite-element.cpp:809
bool has_tensor_product_factorisation() const
Definition: finite-element.cpp:1133
~FiniteElement()=default
Destructor.
const std::vector< std::vector< int > > & num_entity_closure_dofs() const
Definition: finite-element.cpp:849
int dim() const
Definition: finite-element.cpp:814
void apply_inverse_transpose_dof_transformation(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1229
std::array< std::size_t, 4 > tabulate_shape(std::size_t nd, std::size_t num_points) const
Definition: finite-element.cpp:747
element::dpc_variant dpc_variant() const
Definition: finite-element.cpp:1143
element::family family() const
Definition: finite-element.cpp:816
const xt::xtensor< double, 2 > & wcoeffs() const
Definition: finite-element.cpp:1102
const xt::xtensor< double, 2 > & coefficient_matrix() const
Definition: finite-element.cpp:1123
void permute_dofs(const xtl::span< std::int32_t > &dofs, std::uint32_t cell_info) const
Definition: finite-element.cpp:980
bool dof_transformations_are_identity() const
Definition: finite-element.cpp:827
FiniteElement(element::family family, cell::type cell_type, int degree, const std::vector< std::size_t > &value_shape, const xt::xtensor< double, 2 > &wcoeffs, const std::array< std::vector< xt::xtensor< double, 2 >>, 4 > &x, const std::array< std::vector< xt::xtensor< double, 3 >>, 4 > &M, maps::type map_type, bool discontinuous, int highest_complete_degree, std::vector< std::tuple< std::vector< FiniteElement >, std::vector< int >>> tensor_factors={}, element::lagrange_variant lvariant=element::lagrange_variant::unset, element::dpc_variant dvariant=element::dpc_variant::unset)
Definition: finite-element.cpp:470
FiniteElement(const FiniteElement &element)=default
Copy constructor.
bool interpolation_is_identity() const
Definition: finite-element.cpp:1145
void apply_dof_transformation_to_transpose(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1322
void unpermute_dofs(const xtl::span< std::int32_t > &dofs, std::uint32_t cell_info) const
Definition: finite-element.cpp:1036
void apply_transpose_dof_transformation_to_transpose(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1417
void apply_dof_transformation(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1139
bool discontinuous() const
Definition: finite-element.cpp:820
void apply_inverse_dof_transformation(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1276
xt::xtensor< double, 3 > base_transformations() const
Definition: finite-element.cpp:860
type
Cell type.
Definition: cell.h:19
std::tuple< std::array< std::vector< xt::xtensor< double, 2 > >, 4 >, std::array< std::vector< xt::xtensor< double, 3 > >, 4 > > make_discontinuous(const std::array< std::vector< xt::xtensor< double, 2 >>, 4 > &x, const std::array< std::vector< xt::xtensor< double, 3 >>, 4 > &M, int tdim, int value_size)
Definition: finite-element.cpp:341
lagrange_variant
Variants of a Lagrange space that can be created.
Definition: element-families.h:12
dpc_variant
Variants of a DPC space that can be created.
Definition: element-families.h:30
family
Available element families.
Definition: element-families.h:43
void l2_piola(O &&r, const P &U, const Q &, double detJ, const R &)
L2 Piola map.
Definition: maps.h:63
void covariant_piola(O &&r, const P &U, const Q &, double, const R &K)
Covariant Piola map.
Definition: maps.h:71
void contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Contravariant Piola map.
Definition: maps.h:85
void double_contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Double contravariant Piola map.
Definition: maps.h:115
void double_covariant_piola(O &&r, const P &U, const Q &J, double, const R &K)
Double covariant Piola map.
Definition: maps.h:100
type
Map type.
Definition: maps.h:20
void apply_matrix(const std::tuple< std::vector< std::size_t >, std::vector< T >, xt::xtensor< T, 2 >> &matrix, const xtl::span< E > &data, std::size_t offset=0, std::size_t block_size=1)
Definition: precompute.h:319
void apply_matrix_to_transpose(const std::tuple< std::vector< std::size_t >, std::vector< T >, xt::xtensor< T, 2 >> &matrix, const xtl::span< E > &data, std::size_t offset=0, std::size_t block_size=1)
Definition: precompute.h:351
Basix: FEniCS runtime basis evaluation library.
Definition: cell.h:15
FiniteElement create_element(element::family family, cell::type cell, int degree, element::lagrange_variant lvariant, bool discontinuous)
Definition: finite-element.cpp:294
std::string version()
Definition: finite-element.cpp:1158
FiniteElement create_custom_element(cell::type cell_type, int degree, const std::vector< std::size_t > &value_shape, const xt::xtensor< double, 2 > &wcoeffs, const std::array< std::vector< xt::xtensor< double, 2 >>, 4 > &x, const std::array< std::vector< xt::xtensor< double, 3 >>, 4 > &M, maps::type map_type, bool discontinuous, int highest_complete_degree)
Definition: finite-element.cpp:397