Note: this is documentation for an old release. View the latest documentation at docs.fenicsproject.org/basix/v0.8.0/cpp/finite-element_8h_source.html

Basix 0.4.0

Home     Installation     Demos     C++ docs     Python docs

finite-element.h
1 // Copyright (c) 2020 Chris Richardson
2 // FEniCS Project
3 // SPDX-License-Identifier: MIT
4 
5 #pragma once
6 
7 #include "cell.h"
8 #include "element-families.h"
9 #include "maps.h"
10 #include "precompute.h"
11 #include <array>
12 #include <functional>
13 #include <numeric>
14 #include <string>
15 #include <vector>
16 #include <xtensor/xtensor.hpp>
17 #include <xtensor/xview.hpp>
18 #include <xtl/xspan.hpp>
19 
21 namespace basix
22 {
23 
24 namespace element
25 {
26 
39 std::tuple<std::array<std::vector<xt::xtensor<double, 2>>, 4>,
40  std::array<std::vector<xt::xtensor<double, 3>>, 4>>
41 make_discontinuous(const std::array<std::vector<xt::xtensor<double, 2>>, 4>& x,
42  const std::array<std::vector<xt::xtensor<double, 3>>, 4>& M,
43  int tdim, int value_size);
44 
45 } // namespace element
46 
48 
53 {
54 
55 public:
221  const std::vector<std::size_t>& value_shape,
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,
225  maps::type map_type, bool discontinuous, int highest_complete_degree,
226  std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
227  tensor_factors
228  = {},
229  element::lagrange_variant lvariant = element::lagrange_variant::unset,
230  element::dpc_variant dvariant = element::dpc_variant::unset);
231 
233  FiniteElement(const FiniteElement& element) = default;
234 
236  FiniteElement(FiniteElement&& element) = default;
237 
239  ~FiniteElement() = default;
240 
242  FiniteElement& operator=(const FiniteElement& element) = default;
243 
245  FiniteElement& operator=(FiniteElement&& element) = default;
246 
251  bool operator==(const FiniteElement& e) const;
252 
262  std::array<std::size_t, 4> tabulate_shape(std::size_t nd,
263  std::size_t num_points) const;
264 
286  xt::xtensor<double, 4> tabulate(int nd, const xt::xarray<double>& x) const;
287 
313  void tabulate(int nd, const xt::xarray<double>& x,
314  xt::xtensor<double, 4>& basis) const;
315 
318  cell::type cell_type() const;
319 
322  int degree() const;
323 
327  const std::vector<int>& value_shape() const;
328 
332  int dim() const;
333 
336  element::family family() const;
337 
341 
345 
348  maps::type map_type() const;
349 
353  bool discontinuous() const;
354 
358 
362 
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;
380 
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;
391 
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
425  {
426  switch (_map_type)
427  {
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) {
432  maps::l2_piola(u, U, J, detJ, K);
433  };
434  case maps::type::covariantPiola:
435  return [](O& u, const P& U, const Q& J, double detJ, const R& K) {
436  maps::covariant_piola(u, U, J, detJ, K);
437  };
438  case maps::type::contravariantPiola:
439  return [](O& u, const P& U, const Q& J, double detJ, const R& K) {
440  maps::contravariant_piola(u, U, J, detJ, K);
441  };
442  case maps::type::doubleCovariantPiola:
443  return [](O& u, const P& U, const Q& J, double detJ, const R& K) {
444  maps::double_covariant_piola(u, U, J, detJ, K);
445  };
446  case maps::type::doubleContravariantPiola:
447  return [](O& u, const P& U, const Q& J, double detJ, const R& K) {
448  maps::double_contravariant_piola(u, U, J, detJ, K);
449  };
450  default:
451  throw std::runtime_error("Map not implemented");
452  }
453  }
454 
467  const std::vector<std::vector<int>>& num_entity_dofs() const;
468 
474  const std::vector<std::vector<int>>& num_entity_closure_dofs() const;
475 
482  const std::vector<std::vector<std::vector<int>>>& entity_dofs() const;
483 
490  const std::vector<std::vector<std::vector<int>>>& entity_closure_dofs() const;
491 
571  xt::xtensor<double, 3> base_transformations() const;
572 
575  std::map<cell::type, xt::xtensor<double, 3>> entity_transformations() const;
576 
584  void permute_dofs(const xtl::span<std::int32_t>& dofs,
585  std::uint32_t cell_info) const;
586 
594  void unpermute_dofs(const xtl::span<std::int32_t>& dofs,
595  std::uint32_t cell_info) const;
596 
605  template <typename T>
606  void apply_dof_transformation(const xtl::span<T>& data, int block_size,
607  std::uint32_t cell_info) const;
608 
617  template <typename T>
618  void apply_transpose_dof_transformation(const xtl::span<T>& data,
619  int block_size,
620  std::uint32_t cell_info) const;
621 
630  template <typename T>
632  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const;
633 
642  template <typename T>
643  void apply_inverse_dof_transformation(const xtl::span<T>& data,
644  int block_size,
645  std::uint32_t cell_info) const;
646 
655  template <typename T>
656  void apply_dof_transformation_to_transpose(const xtl::span<T>& data,
657  int block_size,
658  std::uint32_t cell_info) const;
659 
668  template <typename T>
670  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const;
671 
680  template <typename T>
682  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const;
683 
692  template <typename T>
694  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const;
695 
700  const xt::xtensor<double, 2>& points() const;
701 
751  const xt::xtensor<double, 2>& interpolation_matrix() const;
752 
758  const xt::xtensor<double, 2>& dual_matrix() const;
759 
791  const xt::xtensor<double, 2>& wcoeffs() const;
792 
796  const std::array<std::vector<xt::xtensor<double, 2>>, 4>& x() const;
797 
830  const std::array<std::vector<xt::xtensor<double, 3>>, 4>& M() const;
831 
837  const xt::xtensor<double, 2>& coefficient_matrix() const;
838 
844  std::array<int, 2> degree_bounds() const;
845 
849 
861  std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
863 
866  bool interpolation_is_identity() const;
867 
868 private:
869  // Cell type
870  cell::type _cell_type;
871 
872  // Topological dimension of the cell
873  std::size_t _cell_tdim;
874 
875  // Topological dimension of the cell
876  std::vector<std::vector<cell::type>> _cell_subentity_types;
877 
878  // Finite element family
879  element::family _family;
880 
881  // Lagrange variant
882  element::lagrange_variant _lagrange_variant;
883 
884  // Lagrange variant
885  element::dpc_variant _dpc_variant;
886 
887  // Degree
888  int _degree;
889 
890  // Value shape
891  std::vector<int> _value_shape;
892 
894  maps::type _map_type;
895 
896  // Shape function coefficient of expansion sets on cell. If shape
897  // function is given by @f$\psi_i = \sum_{k} \phi_{k}
898  // \alpha^{i}_{k}@f$, then _coeffs(i, j) = @f$\alpha^i_k@f$. ie
899  // _coeffs.row(i) are the expansion coefficients for shape function i
900  // (@f$\psi_{i}@f$).
901  xt::xtensor<double, 2> _coeffs;
902 
903  // Number of dofs associated with each cell (sub-)entity
904  //
905  // The dofs of an element are associated with entities of different
906  // topological dimension (vertices, edges, faces, cells). The dofs are
907  // listed in this order, with vertex dofs first. Each entry is the dof
908  // count on the associated entity, as listed by cell::topology.
909  std::vector<std::vector<int>> _num_edofs;
910 
911  // Number of dofs associated with the closure of each cell
912  // (sub-)entity
913  std::vector<std::vector<int>> _num_e_closure_dofs;
914 
915  // Dofs associated with each cell (sub-)entity
916  std::vector<std::vector<std::vector<int>>> _edofs;
917 
918  // Dofs associated with each cell (sub-)entity
919  std::vector<std::vector<std::vector<int>>> _e_closure_dofs;
920 
921  // Entity transformations
922  std::map<cell::type, xt::xtensor<double, 3>> _entity_transformations;
923 
924  // Set of points used for point evaluation
925  // Experimental - currently used for an implementation of
926  // "tabulate_dof_coordinates" Most useful for Lagrange. This may change or go
927  // away. For non-Lagrange elements, these points will be used in combination
928  // with _interpolation_matrix to perform interpolation
929  xt::xtensor<double, 2> _points;
930 
931  // Interpolation points on the cell. The shape is (entity_dim, num
932  // entities of given dimension, num_points, tdim)
933  std::array<std::vector<xt::xtensor<double, 2>>, 4> _x;
934 
936  xt::xtensor<double, 2> _matM;
937 
938  // Indicates whether or not the DOF transformations are all
939  // permutations
940  bool _dof_transformations_are_permutations;
941 
942  // Indicates whether or not the DOF transformations are all identity
943  bool _dof_transformations_are_identity;
944 
945  // The entity permutations (factorised). This will only be set if
946  // _dof_transformations_are_permutations is True and
947  // _dof_transformations_are_identity is False
948  std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm;
949 
950  // The reverse entity permutations (factorised). This will only be set
951  // if _dof_transformations_are_permutations is True and
952  // _dof_transformations_are_identity is False
953  std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_rev;
954 
955  // The entity transformations in precomputed form
956  std::map<cell::type,
957  std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
958  xt::xtensor<double, 2>>>>
959  _etrans;
960 
961  // The transposed entity transformations in precomputed form
962  std::map<cell::type,
963  std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
964  xt::xtensor<double, 2>>>>
965  _etransT;
966 
967  // The inverse entity transformations in precomputed form
968  std::map<cell::type,
969  std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
970  xt::xtensor<double, 2>>>>
971  _etrans_inv;
972 
973  // The inverse transpose entity transformations in precomputed form
974  std::map<cell::type,
975  std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
976  xt::xtensor<double, 2>>>>
977  _etrans_invT;
978 
979  // Indicates whether or not this is the discontinuous version of the
980  // element
981  bool _discontinuous;
982 
983  // The dual matrix
984  xt::xtensor<double, 2> _dual_matrix;
985 
986  // Polynomial degree bounds
987  // [0]: highest degree n such that Lagrange order n is a subspace of
988  // this space
989  // [1]: highest polynomial degree
990  std::array<int, 2> _degree_bounds;
991 
992  // Tensor product representation
993  // Entries of tuple are (list of elements on an interval, permutation
994  // of DOF numbers)
995  // @todo: For vector-valued elements, a tensor product type and a
996  // scaling factor may additionally be needed.
997  std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
998  _tensor_factors;
999 
1000  // Is the interpolation matrix an identity?
1001  bool _interpolation_is_identity;
1002 
1003  // The coefficients that define the polynomial set in terms of the orthonormal
1004  // polynomials
1005  xt::xtensor<double, 2> _wcoeffs;
1006 
1007  // Interpolation matrices for each entity
1008  std::array<std::vector<xt::xtensor<double, 3>>, 4> _M;
1009 };
1010 
1030 FiniteElement create_custom_element(
1031  cell::type cell_type, int degree,
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);
1037 
1047 FiniteElement create_element(element::family family, cell::type cell,
1048  int degree, element::lagrange_variant lvariant,
1049  bool discontinuous);
1050 
1061 FiniteElement create_element(element::family family, cell::type cell,
1062  int degree, element::lagrange_variant lvariant,
1063  element::dpc_variant dvariant, bool discontinuous);
1064 
1074 FiniteElement create_element(element::family family, cell::type cell,
1075  int degree, element::dpc_variant dvariant,
1076  bool discontinuous);
1077 
1088 FiniteElement create_element(element::family family, cell::type cell,
1089  int degree, bool discontinuous);
1090 
1098 FiniteElement create_element(element::family family, cell::type cell,
1099  int degree, element::lagrange_variant lvariant);
1100 
1110 FiniteElement create_element(element::family family, cell::type cell,
1111  int degree, element::lagrange_variant lvariant,
1112  element::dpc_variant dvariant);
1113 
1121 FiniteElement create_element(element::family family, cell::type cell,
1122  int degree, element::dpc_variant dvariant);
1123 
1130 FiniteElement create_element(element::family family, cell::type cell,
1131  int degree);
1132 
1135 std::string version();
1136 
1137 //-----------------------------------------------------------------------------
1138 template <typename T>
1139 void FiniteElement::apply_dof_transformation(const xtl::span<T>& data,
1140  int block_size,
1141  std::uint32_t cell_info) const
1142 {
1143  if (_dof_transformations_are_identity)
1144  return;
1145 
1146  if (_cell_tdim >= 2)
1147  {
1148  // This assumes 3 bits are used per face. This will need updating if
1149  // 3D cells with faces with more than 4 sides are implemented
1150  int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1151  int dofstart
1152  = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1153 
1154  // Transform DOFs on edges
1155  for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1156  {
1157  // Reverse an edge
1158  if (cell_info >> (face_start + e) & 1)
1159  precompute::apply_matrix(_etrans.at(cell::type::interval)[0], data,
1160  dofstart, block_size);
1161  dofstart += _num_edofs[1][e];
1162  }
1163 
1164  if (_cell_tdim == 3)
1165  {
1166  // Permute DOFs on faces
1167  for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1168  {
1169  // Reflect a face
1170  if (cell_info >> (3 * f) & 1)
1171  precompute::apply_matrix(_etrans.at(_cell_subentity_types[2][f])[1],
1172  data, dofstart, block_size);
1173 
1174  // Rotate a face
1175  for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1176  precompute::apply_matrix(_etrans.at(_cell_subentity_types[2][f])[0],
1177  data, dofstart, block_size);
1178  dofstart += _num_edofs[2][f];
1179  }
1180  }
1181  }
1182 }
1183 //-----------------------------------------------------------------------------
1184 template <typename T>
1186  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const
1187 {
1188  if (_dof_transformations_are_identity)
1189  return;
1190 
1191  if (_cell_tdim >= 2)
1192  {
1193  // This assumes 3 bits are used per face. This will need updating if
1194  // 3D cells with faces with more than 4 sides are implemented
1195  int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1196  int dofstart
1197  = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1198 
1199  // Transform DOFs on edges
1200  for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1201  {
1202  // Reverse an edge
1203  if (cell_info >> (face_start + e) & 1)
1204  precompute::apply_matrix(_etransT.at(cell::type::interval)[0], data,
1205  dofstart, block_size);
1206  dofstart += _num_edofs[1][e];
1207  }
1208 
1209  if (_cell_tdim == 3)
1210  {
1211  // Permute DOFs on faces
1212  for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1213  {
1214  // Rotate a face
1215  for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1216  precompute::apply_matrix(_etransT.at(_cell_subentity_types[2][f])[0],
1217  data, dofstart, block_size);
1218  // Reflect a face
1219  if (cell_info >> (3 * f) & 1)
1220  precompute::apply_matrix(_etransT.at(_cell_subentity_types[2][f])[1],
1221  data, dofstart, block_size);
1222  dofstart += _num_edofs[2][f];
1223  }
1224  }
1225  }
1226 }
1227 //-----------------------------------------------------------------------------
1228 template <typename T>
1230  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const
1231 {
1232  if (_dof_transformations_are_identity)
1233  return;
1234 
1235  if (_cell_tdim >= 2)
1236  {
1237  // This assumes 3 bits are used per face. This will need updating if
1238  // 3D cells with faces with more than 4 sides are implemented
1239  int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1240  int dofstart
1241  = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1242 
1243  // Transform DOFs on edges
1244  for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1245  {
1246  // Reverse an edge
1247  if (cell_info >> (face_start + e) & 1)
1248  precompute::apply_matrix(_etrans_invT.at(cell::type::interval)[0], data,
1249  dofstart, block_size);
1250  dofstart += _num_edofs[1][e];
1251  }
1252 
1253  if (_cell_tdim == 3)
1254  {
1255  // Permute DOFs on faces
1256  for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1257  {
1258  // Reflect a face
1259  if (cell_info >> (3 * f) & 1)
1261  _etrans_invT.at(_cell_subentity_types[2][f])[1], data, dofstart,
1262  block_size);
1263 
1264  // Rotate a face
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,
1268  block_size);
1269  dofstart += _num_edofs[2][f];
1270  }
1271  }
1272  }
1273 }
1274 //-----------------------------------------------------------------------------
1275 template <typename T>
1277  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const
1278 {
1279  if (_dof_transformations_are_identity)
1280  return;
1281 
1282  if (_cell_tdim >= 2)
1283  {
1284  // This assumes 3 bits are used per face. This will need updating if
1285  // 3D cells with faces with more than 4 sides are implemented
1286  int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1287  int dofstart
1288  = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1289 
1290  // Transform DOFs on edges
1291  for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1292  {
1293  // Reverse an edge
1294  if (cell_info >> (face_start + e) & 1)
1295  precompute::apply_matrix(_etrans_inv.at(cell::type::interval)[0], data,
1296  dofstart, block_size);
1297  dofstart += _num_edofs[1][e];
1298  }
1299 
1300  if (_cell_tdim == 3)
1301  {
1302  // Permute DOFs on faces
1303  for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1304  {
1305  // Rotate a face
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,
1309  block_size);
1310  // Reflect a face
1311  if (cell_info >> (3 * f) & 1)
1313  _etrans_inv.at(_cell_subentity_types[2][f])[1], data, dofstart,
1314  block_size);
1315  dofstart += _num_edofs[2][f];
1316  }
1317  }
1318  }
1319 }
1320 //-----------------------------------------------------------------------------
1321 template <typename T>
1323  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const
1324 {
1325  if (_dof_transformations_are_identity)
1326  return;
1327 
1328  if (_cell_tdim >= 2)
1329  {
1330  // This assumes 3 bits are used per face. This will need updating if
1331  // 3D cells with faces with more than 4 sides are implemented
1332  int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1333  int dofstart
1334  = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1335 
1336  // Transform DOFs on edges
1337  for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1338  {
1339  // Reverse an edge
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];
1344  }
1345 
1346  if (_cell_tdim == 3)
1347  {
1348  // Permute DOFs on faces
1349  for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1350  {
1351  // Reflect a face
1352  if (cell_info >> (3 * f) & 1)
1354  _etrans.at(_cell_subentity_types[2][f])[1], data, dofstart,
1355  block_size);
1356 
1357  // Rotate a face
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,
1361  block_size);
1362  dofstart += _num_edofs[2][f];
1363  }
1364  }
1365  }
1366 }
1367 //-----------------------------------------------------------------------------
1368 template <typename T>
1370  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const
1371 {
1372  if (_dof_transformations_are_identity)
1373  return;
1374 
1375  if (_cell_tdim >= 2)
1376  {
1377  // This assumes 3 bits are used per face. This will need updating if
1378  // 3D cells with faces with more than 4 sides are implemented
1379  int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1380  int dofstart
1381  = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1382 
1383  // Transform DOFs on edges
1384  for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1385  {
1386  // Reverse an edge
1387  if (cell_info >> (face_start + e) & 1)
1389  _etrans_invT.at(cell::type::interval)[0], data, dofstart,
1390  block_size);
1391  dofstart += _num_edofs[1][e];
1392  }
1393 
1394  if (_cell_tdim == 3)
1395  {
1396  // Permute DOFs on faces
1397  for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1398  {
1399  // Reflect a face
1400  if (cell_info >> (3 * f) & 1)
1402  _etrans_invT.at(_cell_subentity_types[2][f])[1], data, dofstart,
1403  block_size);
1404 
1405  // Rotate a face
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,
1409  block_size);
1410  dofstart += _num_edofs[2][f];
1411  }
1412  }
1413  }
1414 }
1415 //-----------------------------------------------------------------------------
1416 template <typename T>
1418  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const
1419 {
1420  if (_dof_transformations_are_identity)
1421  return;
1422 
1423  if (_cell_tdim >= 2)
1424  {
1425  // This assumes 3 bits are used per face. This will need updating if
1426  // 3D cells with faces with more than 4 sides are implemented
1427  int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1428  int dofstart
1429  = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1430 
1431  // Transform DOFs on edges
1432  for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1433  {
1434  // Reverse an edge
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];
1439  }
1440 
1441  if (_cell_tdim == 3)
1442  {
1443  // Permute DOFs on faces
1444  for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1445  {
1446  // Rotate a face
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,
1450  block_size);
1451 
1452  // Reflect a face
1453  if (cell_info >> (3 * f) & 1)
1455  _etransT.at(_cell_subentity_types[2][f])[1], data, dofstart,
1456  block_size);
1457 
1458  dofstart += _num_edofs[2][f];
1459  }
1460  }
1461  }
1462 }
1463 //-----------------------------------------------------------------------------
1464 template <typename T>
1466  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const
1467 {
1468  if (_dof_transformations_are_identity)
1469  return;
1470 
1471  if (_cell_tdim >= 2)
1472  {
1473  // This assumes 3 bits are used per face. This will need updating if
1474  // 3D cells with faces with more than 4 sides are implemented
1475  int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1476  int dofstart
1477  = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1478 
1479  // Transform DOFs on edges
1480  for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1481  {
1482  // Reverse an edge
1483  if (cell_info >> (face_start + e) & 1)
1485  _etrans_inv.at(cell::type::interval)[0], data, dofstart,
1486  block_size);
1487  dofstart += _num_edofs[1][e];
1488  }
1489 
1490  if (_cell_tdim == 3)
1491  {
1492  // Permute DOFs on faces
1493  for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1494  {
1495  // Rotate a face
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,
1499  block_size);
1500 
1501  // Reflect a face
1502  if (cell_info >> (3 * f) & 1)
1504  _etrans_inv.at(_cell_subentity_types[2][f])[1], data, dofstart,
1505  block_size);
1506 
1507  dofstart += _num_edofs[2][f];
1508  }
1509  }
1510  }
1511 }
1512 //-----------------------------------------------------------------------------
1513 
1514 } // namespace basix
basix::FiniteElement::map_type
maps::type map_type() const
Definition: finite-element.cpp:818
basix::FiniteElement::permute_dofs
void permute_dofs(const xtl::span< std::int32_t > &dofs, std::uint32_t cell_info) const
Definition: finite-element.cpp:980
basix::FiniteElement::dof_transformations_are_identity
bool dof_transformations_are_identity() const
Definition: finite-element.cpp:827
basix::FiniteElement::M
const std::array< std::vector< xt::xtensor< double, 3 > >, 4 > & M() const
Definition: finite-element.cpp:1116
basix::FiniteElement::tabulate
xt::xtensor< double, 4 > tabulate(int nd, const xt::xarray< double > &x) const
Definition: finite-element.cpp:761
basix::FiniteElement::lagrange_variant
element::lagrange_variant lagrange_variant() const
Definition: finite-element.cpp:1138
basix::element::lagrange_variant
lagrange_variant
Variants of a Lagrange space that can be created.
Definition: element-families.h:11
basix::FiniteElement::apply_transpose_dof_transformation
void apply_transpose_dof_transformation(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1185
basix::FiniteElement::apply_dof_transformation
void apply_dof_transformation(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1139
basix::FiniteElement::FiniteElement
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
basix::FiniteElement::degree
int degree() const
Definition: finite-element.cpp:807
basix::FiniteElement::~FiniteElement
~FiniteElement()=default
Destructor.
basix::FiniteElement::dpc_variant
element::dpc_variant dpc_variant() const
Definition: finite-element.cpp:1143
basix::maps::l2_piola
void l2_piola(O &&r, const P &U, const Q &, double detJ, const R &)
L2 Piola map.
Definition: maps.h:63
basix::create_custom_element
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
basix::FiniteElement::num_entity_dofs
const std::vector< std::vector< int > > & num_entity_dofs() const
Definition: finite-element.cpp:837
basix::cell::type
type
Cell type.
Definition: cell.h:18
basix::FiniteElement::points
const xt::xtensor< double, 2 > & points() const
Definition: finite-element.cpp:921
basix::FiniteElement::push_forward
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
basix::FiniteElement::has_tensor_product_factorisation
bool has_tensor_product_factorisation() const
Definition: finite-element.cpp:1133
basix::FiniteElement::entity_transformations
std::map< cell::type, xt::xtensor< double, 3 > > entity_transformations() const
Definition: finite-element.cpp:1092
basix::FiniteElement::cell_type
cell::type cell_type() const
Definition: finite-element.cpp:805
basix::FiniteElement::apply_inverse_transpose_dof_transformation_to_transpose
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
basix::FiniteElement::entity_dofs
const std::vector< std::vector< std::vector< int > > > & entity_dofs() const
Definition: finite-element.cpp:843
basix::maps::double_contravariant_piola
void double_contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Double contravariant Piola map.
Definition: maps.h:115
basix::FiniteElement::discontinuous
bool discontinuous() const
Definition: finite-element.cpp:820
basix::FiniteElement::apply_inverse_dof_transformation
void apply_inverse_dof_transformation(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1276
basix::FiniteElement::apply_transpose_dof_transformation_to_transpose
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
basix::create_element
FiniteElement create_element(element::family family, cell::type cell, int degree, element::lagrange_variant lvariant, bool discontinuous)
Definition: finite-element.cpp:294
basix::FiniteElement::dof_transformations_are_permutations
bool dof_transformations_are_permutations() const
Definition: finite-element.cpp:822
basix::FiniteElement::family
element::family family() const
Definition: finite-element.cpp:816
basix::FiniteElement::apply_inverse_dof_transformation_to_transpose
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
basix::FiniteElement::get_tensor_product_representation
std::vector< std::tuple< std::vector< FiniteElement >, std::vector< int > > > get_tensor_product_representation() const
Definition: finite-element.cpp:1151
basix::element::make_discontinuous
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
basix::FiniteElement::dim
int dim() const
Definition: finite-element.cpp:814
basix::FiniteElement::operator==
bool operator==(const FiniteElement &e) const
Definition: finite-element.cpp:738
basix::FiniteElement::apply_dof_transformation_to_transpose
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
basix::FiniteElement::tabulate_shape
std::array< std::size_t, 4 > tabulate_shape(std::size_t nd, std::size_t num_points) const
Definition: finite-element.cpp:747
basix::maps::type
type
Map type.
Definition: maps.h:19
basix::FiniteElement::value_shape
const std::vector< int > & value_shape() const
Definition: finite-element.cpp:809
basix::FiniteElement::base_transformations
xt::xtensor< double, 3 > base_transformations() const
Definition: finite-element.cpp:860
basix::FiniteElement::interpolation_is_identity
bool interpolation_is_identity() const
Definition: finite-element.cpp:1145
basix::FiniteElement::num_entity_closure_dofs
const std::vector< std::vector< int > > & num_entity_closure_dofs() const
Definition: finite-element.cpp:849
basix::maps::double_covariant_piola
void double_covariant_piola(O &&r, const P &U, const Q &J, double, const R &K)
Double covariant Piola map.
Definition: maps.h:100
basix::precompute::apply_matrix_to_transpose
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::FiniteElement
A finite element.
Definition: finite-element.h:52
basix::FiniteElement::interpolation_matrix
const xt::xtensor< double, 2 > & interpolation_matrix() const
Definition: finite-element.cpp:832
basix::FiniteElement::wcoeffs
const xt::xtensor< double, 2 > & wcoeffs() const
Definition: finite-element.cpp:1102
basix::maps::contravariant_piola
void contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Contravariant Piola map.
Definition: maps.h:85
basix::FiniteElement::coefficient_matrix
const xt::xtensor< double, 2 > & coefficient_matrix() const
Definition: finite-element.cpp:1123
basix::FiniteElement::dual_matrix
const xt::xtensor< double, 2 > & dual_matrix() const
Definition: finite-element.cpp:1097
basix::version
std::string version()
Definition: finite-element.cpp:1158
basix::FiniteElement::entity_closure_dofs
const std::vector< std::vector< std::vector< int > > > & entity_closure_dofs() const
Definition: finite-element.cpp:855
basix::FiniteElement::map_fn
std::function< void(O &, const P &, const Q &, double, const R &)> map_fn() const
Definition: finite-element.h:424
basix
Basix: FEniCS runtime basis evaluation library.
Definition: cell.h:14
basix::FiniteElement::unpermute_dofs
void unpermute_dofs(const xtl::span< std::int32_t > &dofs, std::uint32_t cell_info) const
Definition: finite-element.cpp:1036
basix::FiniteElement::degree_bounds
std::array< int, 2 > degree_bounds() const
Definition: finite-element.cpp:1128
basix::precompute::apply_matrix
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
basix::FiniteElement::operator=
FiniteElement & operator=(const FiniteElement &element)=default
Assignment operator.
basix::FiniteElement::x
const std::array< std::vector< xt::xtensor< double, 2 > >, 4 > & x() const
Definition: finite-element.cpp:1110
basix::FiniteElement::pull_back
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
basix::maps::covariant_piola
void covariant_piola(O &&r, const P &U, const Q &, double, const R &K)
Covariant Piola map.
Definition: maps.h:71
basix::element::dpc_variant
dpc_variant
Variants of a DPC space that can be created.
Definition: element-families.h:29
basix::FiniteElement::apply_inverse_transpose_dof_transformation
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
basix::element::family
family
Available element families.
Definition: element-families.h:42