11 #include "element-families.h"
14 #include "precompute.h"
20 #include <xtensor/xadapt.hpp>
21 #include <xtensor/xcomplex.hpp>
22 #include <xtensor/xtensor.hpp>
23 #include <xtensor/xview.hpp>
24 #include <xtl/xspan.hpp>
186 cell::type cell_type,
const xt::xtensor<double, 2>& B,
187 const std::vector<std::vector<xt::xtensor<double, 3>>>& M,
188 const std::vector<std::vector<xt::xtensor<double, 2>>>& x,
int degree,
189 double kappa_tol = 0.0);
214 const xt::xtensor<double, 3>& coeffs,
215 const std::map<
cell::type, xt::xtensor<double, 3>>&
217 const std::array<std::vector<xt::xtensor<double, 2>>, 4>& x,
218 const std::array<std::vector<xt::xtensor<double, 3>>, 4>& M,
253 xt::xtensor<double, 4>
tabulate(
int nd,
const xt::xarray<double>& x)
const;
259 void tabulate(
int nd,
const xt::xarray<double>& x,
260 xt::xtensor<double, 4>& basis_data)
const;
286 element::family
family()
const;
314 xt::xtensor<double, 3>
316 const xt::xtensor<double, 3>& J,
317 const xtl::span<const double>& detJ,
318 const xt::xtensor<double, 3>& K)
const;
335 template <
typename O,
typename P,
typename Q,
typename S,
typename T>
343 for (std::size_t i = 0; i < U.shape(0); ++i)
345 auto _J = xt::view(J, i, xt::all(), xt::all());
346 auto _K = xt::view(K, i, xt::all(), xt::all());
347 auto _U = xt::view(U, i, xt::all(), xt::all());
348 auto _u = xt::view(u, i, xt::all(), xt::all());
359 xt::xtensor<double, 3>
map_pull_back(
const xt::xtensor<double, 3>& u,
360 const xt::xtensor<double, 3>& J,
361 const xtl::span<const double>& detJ,
362 const xt::xtensor<double, 3>& K)
const;
379 template <
typename O,
typename P,
typename Q,
typename S,
typename T>
384 for (std::size_t i = 0; i < u.shape(0); ++i)
386 auto _J = xt::view(J, i, xt::all(), xt::all());
387 auto _K = xt::view(K, i, xt::all(), xt::all());
388 auto _u = xt::view(u, i, xt::all(), xt::all());
389 auto _U = xt::view(U, i, xt::all(), xt::all());
421 const std::vector<std::vector<std::set<int>>>&
entity_dofs()
const;
518 std::uint32_t cell_info)
const;
524 std::uint32_t cell_info)
const;
530 template <
typename T>
532 std::uint32_t cell_info)
const;
538 template <
typename T>
541 std::uint32_t cell_info)
const;
547 template <
typename T>
549 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const;
555 template <
typename T>
558 std::uint32_t cell_info)
const;
564 template <
typename T>
567 std::uint32_t cell_info)
const;
573 template <
typename T>
575 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const;
581 template <
typename T>
583 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const;
589 template <
typename T>
591 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const;
597 const xt::xtensor<double, 2>&
points()
const;
616 template <
typename T>
618 const xtl::span<const T>& data,
const int block_size)
const;
628 std::size_t _cell_tdim;
631 std::vector<std::vector<cell::type>> _cell_subentity_types;
634 element::family _family;
640 std::vector<int> _value_shape;
650 xt::xtensor<double, 2> _coeffs;
658 std::vector<std::vector<int>> _num_edofs;
661 std::vector<std::vector<int>> _num_e_closure_dofs;
664 std::vector<std::vector<std::set<int>>> _edofs;
667 std::vector<std::vector<std::set<int>>> _e_closure_dofs;
670 std::map<cell::type, xt::xtensor<double, 3>> _entity_transformations;
677 xt::xtensor<double, 2> _points;
681 std::array<std::vector<xt::xtensor<double, 2>>, 4> _x;
684 xt::xtensor<double, 2> _matM;
687 std::array<std::vector<xt::xtensor<double, 3>>, 4> _matM_new;
690 bool _dof_transformations_are_permutations;
693 bool _dof_transformations_are_identity;
698 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm;
703 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_rev;
707 std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
708 xt::xtensor<double, 2>>>>
713 std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
714 xt::xtensor<double, 2>>>>
719 std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
720 xt::xtensor<double, 2>>>>
725 std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
726 xt::xtensor<double, 2>>>>
738 lattice::type lattice_type);
752 template <
typename T>
755 std::uint32_t cell_info)
const
757 if (_dof_transformations_are_identity)
764 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
766 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
769 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
772 if (cell_info >> (face_start + e) & 1)
774 dofstart, block_size);
775 dofstart += _num_edofs[1][e];
781 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
784 if (cell_info >> (3 * f) & 1)
786 data, dofstart, block_size);
789 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
791 data, dofstart, block_size);
792 dofstart += _num_edofs[2][f];
798 template <
typename T>
800 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
802 if (_dof_transformations_are_identity)
809 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
811 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
814 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
817 if (cell_info >> (face_start + e) & 1)
819 dofstart, block_size);
820 dofstart += _num_edofs[1][e];
826 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
829 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
831 data, dofstart, block_size);
833 if (cell_info >> (3 * f) & 1)
835 data, dofstart, block_size);
836 dofstart += _num_edofs[2][f];
842 template <
typename T>
844 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
846 if (_dof_transformations_are_identity)
853 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
855 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
858 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
861 if (cell_info >> (face_start + e) & 1)
863 dofstart, block_size);
864 dofstart += _num_edofs[1][e];
870 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
873 if (cell_info >> (3 * f) & 1)
875 _etrans_invT.at(_cell_subentity_types[2][f])[1], data, dofstart,
879 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
881 _etrans_invT.at(_cell_subentity_types[2][f])[0], data, dofstart,
883 dofstart += _num_edofs[2][f];
889 template <
typename T>
891 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
893 if (_dof_transformations_are_identity)
900 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
902 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
905 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
908 if (cell_info >> (face_start + e) & 1)
910 dofstart, block_size);
911 dofstart += _num_edofs[1][e];
917 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
920 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
922 _etrans_inv.at(_cell_subentity_types[2][f])[0], data, dofstart,
925 if (cell_info >> (3 * f) & 1)
927 _etrans_inv.at(_cell_subentity_types[2][f])[1], data, dofstart,
929 dofstart += _num_edofs[2][f];
935 template <
typename T>
937 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
939 if (_dof_transformations_are_identity)
946 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
948 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
951 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
954 if (cell_info >> (face_start + e) & 1)
956 _etrans.at(cell::type::interval)[0], data, dofstart, block_size);
957 dofstart += _num_edofs[1][e];
963 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
966 if (cell_info >> (3 * f) & 1)
968 _etrans.at(_cell_subentity_types[2][f])[1], data, dofstart,
972 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
974 _etrans.at(_cell_subentity_types[2][f])[0], data, dofstart,
976 dofstart += _num_edofs[2][f];
982 template <
typename T>
984 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
986 if (_dof_transformations_are_identity)
993 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
995 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
998 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1001 if (cell_info >> (face_start + e) & 1)
1003 _etrans_invT.at(cell::type::interval)[0], data, dofstart,
1005 dofstart += _num_edofs[1][e];
1008 if (_cell_tdim == 3)
1011 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1014 if (cell_info >> (3 * f) & 1)
1016 _etrans_invT.at(_cell_subentity_types[2][f])[1], data, dofstart,
1020 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1022 _etrans_invT.at(_cell_subentity_types[2][f])[0], data, dofstart,
1024 dofstart += _num_edofs[2][f];
1030 template <
typename T>
1032 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1034 if (_dof_transformations_are_identity)
1037 if (_cell_tdim >= 2)
1041 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1043 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1046 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1049 if (cell_info >> (face_start + e) & 1)
1051 _etransT.at(cell::type::interval)[0], data, dofstart, block_size);
1052 dofstart += _num_edofs[1][e];
1055 if (_cell_tdim == 3)
1058 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1061 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1063 _etransT.at(_cell_subentity_types[2][f])[0], data, dofstart,
1067 if (cell_info >> (3 * f) & 1)
1069 _etransT.at(_cell_subentity_types[2][f])[1], data, dofstart,
1072 dofstart += _num_edofs[2][f];
1078 template <
typename T>
1080 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1082 if (_dof_transformations_are_identity)
1085 if (_cell_tdim >= 2)
1089 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1091 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1094 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1097 if (cell_info >> (face_start + e) & 1)
1099 _etrans_inv.at(cell::type::interval)[0], data, dofstart,
1101 dofstart += _num_edofs[1][e];
1104 if (_cell_tdim == 3)
1107 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1110 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1112 _etrans_inv.at(_cell_subentity_types[2][f])[0], data, dofstart,
1116 if (cell_info >> (3 * f) & 1)
1118 _etrans_inv.at(_cell_subentity_types[2][f])[1], data, dofstart,
1121 dofstart += _num_edofs[2][f];
1127 template <
typename T>
1129 const xtl::span<const T>& data,
1130 const int block_size)
const
1132 if (block_size != 1)
1134 throw std::runtime_error(
1135 "Interpolation of blocked data not implemented yet.");
1138 const std::size_t rows =
dim();
1142 assert(Pi.size() % rows == 0);
1143 const std::size_t cols = Pi.size() / rows;
1144 for (std::size_t i = 0; i < rows; ++i)
1148 coefficients[i] = std::inner_product(std::next(Pi.data(), i * cols),
1149 std::next(Pi.data(), i * cols + cols),
1150 data.data(), T(0.0));
Definition: finite-element.h:195
FiniteElement & operator=(const FiniteElement &element)=default
Assignment operator.
cell::type cell_type() const
Definition: finite-element.cpp:488
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:548
FiniteElement(element::family family, cell::type cell_type, int degree, const std::vector< std::size_t > &value_shape, const xt::xtensor< double, 3 > &coeffs, const std::map< cell::type, xt::xtensor< double, 3 >> &entity_transformations, 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=maps::type::identity)
Definition: finite-element.cpp:259
int value_size() const
Definition: finite-element.cpp:492
FiniteElement(FiniteElement &&element)=default
Move constructor.
std::map< cell::type, xt::xtensor< double, 3 > > entity_transformations() const
Return the entity dof transformation matricess.
Definition: finite-element.cpp:780
const std::vector< std::vector< std::set< int > > > & entity_dofs() const
Definition: finite-element.cpp:530
bool dof_transformations_are_permutations() const
Definition: finite-element.cpp:509
xt::xtensor< double, 3 > map_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:668
const xt::xtensor< double, 2 > & points() const
Definition: finite-element.cpp:655
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:1079
const xt::xtensor< double, 2 > & interpolation_matrix() const
Definition: finite-element.cpp:519
void apply_transpose_dof_transformation(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:799
const std::vector< std::vector< int > > & num_entity_dofs() const
Definition: finite-element.cpp:524
int degree() const
Definition: finite-element.cpp:490
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:983
const std::vector< int > & value_shape() const
Definition: finite-element.cpp:498
~FiniteElement()=default
Destructor.
const std::vector< std::vector< int > > & num_entity_closure_dofs() const
Definition: finite-element.cpp:536
int num_points() const
Return the number of interpolation points.
Definition: finite-element.cpp:653
int dim() const
Definition: finite-element.cpp:503
void apply_inverse_transpose_dof_transformation(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:843
element::family family() const
Definition: finite-element.cpp:505
void map_push_forward_m(const O &U, const P &J, const Q &detJ, const S &K, T &&u) const
Definition: finite-element.h:336
maps::type map_type
Element map type.
Definition: finite-element.h:621
void permute_dofs(const xtl::span< std::int32_t > &dofs, std::uint32_t cell_info) const
Definition: finite-element.cpp:678
bool dof_transformations_are_identity() const
Definition: finite-element.cpp:514
const std::vector< std::vector< std::set< int > > > & entity_closure_dofs() const
Definition: finite-element.cpp:542
void interpolate(const xtl::span< T > &coefficients, const xtl::span< const T > &data, const int block_size) const
Definition: finite-element.h:1128
FiniteElement(const FiniteElement &element)=default
Copy constructor.
maps::type mapping_type() const
Definition: finite-element.cpp:507
void apply_dof_transformation_to_transpose(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:936
void unpermute_dofs(const xtl::span< std::int32_t > &dofs, std::uint32_t cell_info) const
Definition: finite-element.cpp:729
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:1031
void apply_dof_transformation(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:753
void map_pull_back_m(const O &u, const P &J, const Q &detJ, const S &K, T &&U) const
Definition: finite-element.h:380
xt::xtensor< double, 3 > map_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:657
void apply_inverse_dof_transformation(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:890
xt::xtensor< double, 3 > base_transformations() const
Definition: finite-element.cpp:595
type
Cell type.
Definition: cell.h:17
void apply_map(O &&u, const P &U, const Mat0 &J, double detJ, const Mat1 &K, maps::type map_type)
Definition: maps.h:130
type
Cell 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:311
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:340
Placeholder.
Definition: brezzi-douglas-marini.h:10
FiniteElement create_element(element::family family, cell::type cell, int degree, lattice::type lattice_type)
Definition: finite-element.cpp:146
std::string version()
Definition: finite-element.cpp:785
xt::xtensor< double, 3 > compute_expansion_coefficients(cell::type cell_type, const xt::xtensor< double, 2 > &B, const std::vector< std::vector< xt::xtensor< double, 3 >>> &M, const std::vector< std::vector< xt::xtensor< double, 2 >>> &x, int degree, double kappa_tol=0.0)
Definition: finite-element.cpp:159