# Basix 0.4.0

## HomeInstallationDemosC++ docsPython docs

finite-element.h
1 // Copyright (c) 2020 Chris Richardson
2 // FEniCS Project
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