Basix
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 "mdspan.hpp"
11#include "precompute.h"
12#include "sobolev-spaces.h"
13#include <array>
14#include <cstdint>
15#include <functional>
16#include <map>
17#include <numeric>
18#include <span>
19#include <string>
20#include <tuple>
21#include <vector>
22
24namespace basix
25{
26
27namespace impl
28{
29using mdspan2_t
30 = std::experimental::mdspan<double,
31 std::experimental::dextents<std::size_t, 2>>;
32using mdspan4_t
33 = std::experimental::mdspan<double,
34 std::experimental::dextents<std::size_t, 4>>;
35using cmdspan2_t
36 = std::experimental::mdspan<const double,
37 std::experimental::dextents<std::size_t, 2>>;
38using cmdspan3_t
39 = std::experimental::mdspan<const double,
40 std::experimental::dextents<std::size_t, 3>>;
41using cmdspan4_t
42 = std::experimental::mdspan<const double,
43 std::experimental::dextents<std::size_t, 4>>;
44
45using mdarray2_t
46 = std::experimental::mdarray<double,
47 std::experimental::dextents<std::size_t, 2>>;
48using mdarray4_t
49 = std::experimental::mdarray<double,
50 std::experimental::dextents<std::size_t, 4>>;
51
54inline std::array<std::vector<cmdspan2_t>, 4>
55to_mdspan(std::array<std::vector<mdarray2_t>, 4>& x)
56{
57 std::array<std::vector<cmdspan2_t>, 4> x1;
58 for (std::size_t i = 0; i < x.size(); ++i)
59 for (std::size_t j = 0; j < x[i].size(); ++j)
60 x1[i].emplace_back(x[i][j].data(), x[i][j].extents());
61
62 return x1;
63}
64
67inline std::array<std::vector<cmdspan4_t>, 4>
68to_mdspan(std::array<std::vector<mdarray4_t>, 4>& M)
69{
70 std::array<std::vector<cmdspan4_t>, 4> M1;
71 for (std::size_t i = 0; i < M.size(); ++i)
72 for (std::size_t j = 0; j < M[i].size(); ++j)
73 M1[i].emplace_back(M[i][j].data(), M[i][j].extents());
74
75 return M1;
76}
77
80inline std::array<std::vector<cmdspan2_t>, 4>
81to_mdspan(const std::array<std::vector<std::vector<double>>, 4>& x,
82 const std::array<std::vector<std::array<std::size_t, 2>>, 4>& shape)
83{
84 std::array<std::vector<cmdspan2_t>, 4> x1;
85 for (std::size_t i = 0; i < x.size(); ++i)
86 for (std::size_t j = 0; j < x[i].size(); ++j)
87 x1[i].push_back(cmdspan2_t(x[i][j].data(), shape[i][j]));
88
89 return x1;
90}
91
94inline std::array<std::vector<cmdspan4_t>, 4>
95to_mdspan(const std::array<std::vector<std::vector<double>>, 4>& M,
96 const std::array<std::vector<std::array<std::size_t, 4>>, 4>& shape)
97{
98 std::array<std::vector<cmdspan4_t>, 4> M1;
99 for (std::size_t i = 0; i < M.size(); ++i)
100 for (std::size_t j = 0; j < M[i].size(); ++j)
101 M1[i].push_back(cmdspan4_t(M[i][j].data(), shape[i][j]));
102
103 return M1;
104}
105
106} // namespace impl
107
108namespace element
109{
111using cmdspan2_t = impl::cmdspan2_t;
113using cmdspan4_t = impl::cmdspan4_t;
114
130std::tuple<std::array<std::vector<std::vector<double>>, 4>,
131 std::array<std::vector<std::array<std::size_t, 2>>, 4>,
132 std::array<std::vector<std::vector<double>>, 4>,
133 std::array<std::vector<std::array<std::size_t, 4>>, 4>>
134make_discontinuous(const std::array<std::vector<cmdspan2_t>, 4>& x,
135 const std::array<std::vector<cmdspan4_t>, 4>& M,
136 std::size_t tdim, std::size_t value_size);
137
138} // namespace element
139
146{
147 using cmdspan2_t
148 = std::experimental::mdspan<const double,
149 std::experimental::dextents<std::size_t, 2>>;
150 using cmdspan4_t
151 = std::experimental::mdspan<const double,
152 std::experimental::dextents<std::size_t, 4>>;
153
154public:
327 const std::vector<std::size_t>& value_shape, const cmdspan2_t& wcoeffs,
328 const std::array<std::vector<cmdspan2_t>, 4>& x,
329 const std::array<std::vector<cmdspan4_t>, 4>& M,
334 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
335 tensor_factors
336 = {});
337
341 const std::vector<std::size_t>& value_shape, const cmdspan2_t& wcoeffs,
342 const std::array<std::vector<cmdspan2_t>, 4>& x,
343 const std::array<std::vector<cmdspan4_t>, 4>& M,
348 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
349 tensor_factors
350 = {});
351
355 const std::vector<std::size_t>& value_shape, const cmdspan2_t& wcoeffs,
356 const std::array<std::vector<cmdspan2_t>, 4>& x,
357 const std::array<std::vector<cmdspan4_t>, 4>& M,
361 element::dpc_variant dvariant,
362 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
363 tensor_factors
364 = {});
365
369 const std::vector<std::size_t>& value_shape, const cmdspan2_t& wcoeffs,
370 const std::array<std::vector<cmdspan2_t>, 4>& x,
371 const std::array<std::vector<cmdspan4_t>, 4>& M,
375 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
376 tensor_factors
377 = {});
378
380 FiniteElement(const FiniteElement& element) = default;
381
383 FiniteElement(FiniteElement&& element) = default;
384
386 ~FiniteElement() = default;
387
389 FiniteElement& operator=(const FiniteElement& element) = default;
390
392 FiniteElement& operator=(FiniteElement&& element) = default;
393
398 bool operator==(const FiniteElement& e) const;
399
409 std::array<std::size_t, 4> tabulate_shape(std::size_t nd,
410 std::size_t num_points) const;
411
433 std::pair<std::vector<double>, std::array<std::size_t, 4>>
434 tabulate(int nd, impl::cmdspan2_t x) const;
435
459 std::pair<std::vector<double>, std::array<std::size_t, 4>>
460 tabulate(int nd, const std::span<const double>& x,
461 std::array<std::size_t, 2> shape) const;
462
488 void tabulate(int nd, impl::cmdspan2_t x, impl::mdspan4_t basis) const;
489
515 void tabulate(int nd, const std::span<const double>& x,
516 std::array<std::size_t, 2> xshape,
517 const std::span<double>& basis) const;
518
521 cell::type cell_type() const;
522
525 int degree() const;
526
531 int highest_degree() const;
532
536 int highest_complete_degree() const;
537
541 const std::vector<std::size_t>& value_shape() const;
542
546 int dim() const;
547
550 element::family family() const;
551
555
559
562 maps::type map_type() const;
563
567
571 bool discontinuous() const;
572
576
580
594 std::pair<std::vector<double>, std::array<std::size_t, 3>>
595 push_forward(impl::cmdspan3_t U, impl::cmdspan3_t J,
596 std::span<const double> detJ, impl::cmdspan3_t K) const;
597
605 std::pair<std::vector<double>, std::array<std::size_t, 3>>
606 pull_back(impl::cmdspan3_t u, impl::cmdspan3_t J,
607 std::span<const double> detJ, impl::cmdspan3_t K) const;
608
640 template <typename O, typename P, typename Q, typename R>
641 std::function<void(O&, const P&, const Q&, double, const R&)> map_fn() const
642 {
643 switch (_map_type)
644 {
645 case maps::type::identity:
646 return [](O& u, const P& U, const Q&, double, const R&)
647 {
648 assert(U.extent(0) == u.extent(0));
649 assert(U.extent(1) == u.extent(1));
650 for (std::size_t i = 0; i < U.extent(0); ++i)
651 for (std::size_t j = 0; j < U.extent(1); ++j)
652 u(i, j) = U(i, j);
653 };
654 case maps::type::covariantPiola:
655 return [](O& u, const P& U, const Q& J, double detJ, const R& K)
656 { maps::covariant_piola(u, U, J, detJ, K); };
657 case maps::type::contravariantPiola:
658 return [](O& u, const P& U, const Q& J, double detJ, const R& K)
659 { maps::contravariant_piola(u, U, J, detJ, K); };
660 case maps::type::doubleCovariantPiola:
661 return [](O& u, const P& U, const Q& J, double detJ, const R& K)
662 { maps::double_covariant_piola(u, U, J, detJ, K); };
663 case maps::type::doubleContravariantPiola:
664 return [](O& u, const P& U, const Q& J, double detJ, const R& K)
665 { maps::double_contravariant_piola(u, U, J, detJ, K); };
666 default:
667 throw std::runtime_error("Map not implemented");
668 }
669 }
670
677 const std::vector<std::vector<std::vector<int>>>& entity_dofs() const;
678
686 const std::vector<std::vector<std::vector<int>>>& entity_closure_dofs() const;
687
769 std::pair<std::vector<double>, std::array<std::size_t, 3>>
770 base_transformations() const;
771
776 std::map<cell::type,
777 std::pair<std::vector<double>, std::array<std::size_t, 3>>>
779
787 void permute_dofs(const std::span<std::int32_t>& dofs,
788 std::uint32_t cell_info) const;
789
797 void unpermute_dofs(const std::span<std::int32_t>& dofs,
798 std::uint32_t cell_info) const;
799
808 template <typename T>
809 void apply_dof_transformation(const std::span<T>& data, int block_size,
810 std::uint32_t cell_info) const;
811
820 template <typename T>
821 void apply_transpose_dof_transformation(const std::span<T>& data,
822 int block_size,
823 std::uint32_t cell_info) const;
824
833 template <typename T>
835 const std::span<T>& data, int block_size, std::uint32_t cell_info) const;
836
845 template <typename T>
846 void apply_inverse_dof_transformation(const std::span<T>& data,
847 int block_size,
848 std::uint32_t cell_info) const;
849
858 template <typename T>
859 void apply_dof_transformation_to_transpose(const std::span<T>& data,
860 int block_size,
861 std::uint32_t cell_info) const;
862
871 template <typename T>
873 const std::span<T>& data, int block_size, std::uint32_t cell_info) const;
874
884 template <typename T>
886 const std::span<T>& data, int block_size, std::uint32_t cell_info) const;
887
896 template <typename T>
898 const std::span<T>& data, int block_size, std::uint32_t cell_info) const;
899
904 const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
905 points() const;
906
959 const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
960 interpolation_matrix() const;
961
967 const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
968 dual_matrix() const;
969
1004 const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
1005 wcoeffs() const;
1006
1011 const std::array<
1012 std::vector<std::pair<std::vector<double>, std::array<std::size_t, 2>>>,
1013 4>&
1014 x() const;
1015
1052 const std::array<
1053 std::vector<std::pair<std::vector<double>, std::array<std::size_t, 4>>>,
1054 4>&
1055 M() const;
1056
1062 const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
1063 coefficient_matrix() const;
1064
1073
1085 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
1087
1090 bool interpolation_is_identity() const;
1091
1093 int interpolation_nderivs() const;
1094
1095private:
1096 // Cell type
1097 cell::type _cell_type;
1098
1099 // Topological dimension of the cell
1100 std::size_t _cell_tdim;
1101
1102 // Topological dimension of the cell
1103 std::vector<std::vector<cell::type>> _cell_subentity_types;
1104
1105 // Finite element family
1106 element::family _family;
1107
1108 // Lagrange variant
1109 element::lagrange_variant _lagrange_variant;
1110
1111 // Lagrange variant
1112 element::dpc_variant _dpc_variant;
1113
1114 // Degree that was input when creating the element
1115 int _degree;
1116
1117 // Degree
1118 int _interpolation_nderivs;
1119
1120 // Highest degree polynomial in element's polyset
1121 int _highest_degree;
1122
1123 // Highest degree space that is a subspace of element's polyset
1124 int _highest_complete_degree;
1125
1126 // Value shape
1127 std::vector<std::size_t> _value_shape;
1128
1130 maps::type _map_type;
1131
1133 sobolev::space _sobolev_space;
1134
1135 // Shape function coefficient of expansion sets on cell. If shape
1136 // function is given by @f$\psi_i = \sum_{k} \phi_{k}
1137 // \alpha^{i}_{k}@f$, then _coeffs(i, j) = @f$\alpha^i_k@f$. ie
1138 // _coeffs.row(i) are the expansion coefficients for shape function i
1139 // (@f$\psi_{i}@f$).
1140 std::pair<std::vector<double>, std::array<std::size_t, 2>> _coeffs;
1141
1142 // Dofs associated with each cell (sub-)entity
1143 std::vector<std::vector<std::vector<int>>> _edofs;
1144
1145 // Dofs associated with the closdure of each cell (sub-)entity
1146 std::vector<std::vector<std::vector<int>>> _e_closure_dofs;
1147
1148 using array2_t = std::pair<std::vector<double>, std::array<std::size_t, 2>>;
1149 using array3_t = std::pair<std::vector<double>, std::array<std::size_t, 3>>;
1150
1151 // Entity transformations
1152 std::map<cell::type, array3_t> _entity_transformations;
1153
1154 // Set of points used for point evaluation
1155 // Experimental - currently used for an implementation of
1156 // "tabulate_dof_coordinates" Most useful for Lagrange. This may change or go
1157 // away. For non-Lagrange elements, these points will be used in combination
1158 // with _interpolation_matrix to perform interpolation
1159 std::pair<std::vector<double>, std::array<std::size_t, 2>> _points;
1160
1161 // Interpolation points on the cell. The shape is (entity_dim, num
1162 // entities of given dimension, num_points, tdim)
1163 std::array<
1164 std::vector<std::pair<std::vector<double>, std::array<std::size_t, 2>>>,
1165 4>
1166 _x;
1167
1169 std::pair<std::vector<double>, std::array<std::size_t, 2>> _matM;
1170
1171 // Indicates whether or not the DOF transformations are all
1172 // permutations
1173 bool _dof_transformations_are_permutations;
1174
1175 // Indicates whether or not the DOF transformations are all identity
1176 bool _dof_transformations_are_identity;
1177
1178 // The entity permutations (factorised). This will only be set if
1179 // _dof_transformations_are_permutations is True and
1180 // _dof_transformations_are_identity is False
1181 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm;
1182
1183 // The reverse entity permutations (factorised). This will only be set
1184 // if _dof_transformations_are_permutations is True and
1185 // _dof_transformations_are_identity is False
1186 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_rev;
1187
1188 using trans_data_t
1189 = std::vector<std::pair<std::vector<std::size_t>, array2_t>>;
1190
1191 // The entity transformations in precomputed form
1192 std::map<cell::type, trans_data_t> _etrans;
1193
1194 // The transposed entity transformations in precomputed form
1195 std::map<cell::type, trans_data_t> _etransT;
1196
1197 // The inverse entity transformations in precomputed form
1198 std::map<cell::type, trans_data_t> _etrans_inv;
1199
1200 // The inverse transpose entity transformations in precomputed form
1201 std::map<cell::type, trans_data_t> _etrans_invT;
1202
1203 // Indicates whether or not this is the discontinuous version of the
1204 // element
1205 bool _discontinuous;
1206
1207 // The dual matrix
1208 std::pair<std::vector<double>, std::array<std::size_t, 2>> _dual_matrix;
1209
1210 // Tensor product representation
1211 // Entries of tuple are (list of elements on an interval, permutation
1212 // of DOF numbers)
1213 // @todo: For vector-valued elements, a tensor product type and a
1214 // scaling factor may additionally be needed.
1215 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
1216 _tensor_factors;
1217
1218 // Is the interpolation matrix an identity?
1219 bool _interpolation_is_identity;
1220
1221 // The coefficients that define the polynomial set in terms of the
1222 // orthonormal polynomials
1223 std::pair<std::vector<double>, std::array<std::size_t, 2>> _wcoeffs;
1224
1225 // Interpolation matrices for each entity
1226 using array4_t
1227 = std::vector<std::pair<std::vector<double>, std::array<std::size_t, 4>>>;
1228 std::array<array4_t, 4> _M;
1229 // std::array<
1230 // std::vector<std::pair<std::vector<double>, std::array<std::size_t,
1231 // 4>>>, 4> _M;
1232};
1233
1258FiniteElement
1260 const std::vector<std::size_t>& value_shape,
1261 const impl::cmdspan2_t& wcoeffs,
1262 const std::array<std::vector<impl::cmdspan2_t>, 4>& x,
1263 const std::array<std::vector<impl::cmdspan4_t>, 4>& M,
1264 int interpolation_nderivs, maps::type map_type,
1265 sobolev::space sobolev_space, bool discontinuous,
1266 int highest_complete_degree, int highest_degree);
1267
1277FiniteElement create_element(element::family family, cell::type cell,
1278 int degree, element::lagrange_variant lvariant,
1279 bool discontinuous);
1280
1291FiniteElement create_element(element::family family, cell::type cell,
1292 int degree, element::lagrange_variant lvariant,
1293 element::dpc_variant dvariant, bool discontinuous);
1294
1304FiniteElement create_element(element::family family, cell::type cell,
1305 int degree, element::dpc_variant dvariant,
1306 bool discontinuous);
1307
1318FiniteElement create_element(element::family family, cell::type cell,
1319 int degree, bool discontinuous);
1320
1328FiniteElement create_element(element::family family, cell::type cell,
1329 int degree, element::lagrange_variant lvariant);
1330
1340FiniteElement create_element(element::family family, cell::type cell,
1341 int degree, element::lagrange_variant lvariant,
1342 element::dpc_variant dvariant);
1343
1351FiniteElement create_element(element::family family, cell::type cell,
1352 int degree, element::dpc_variant dvariant);
1353
1360FiniteElement create_element(element::family family, cell::type cell,
1361 int degree);
1362
1365std::string version();
1366
1367//-----------------------------------------------------------------------------
1368template <typename T>
1369void FiniteElement::apply_dof_transformation(const std::span<T>& data,
1370 int block_size,
1371 std::uint32_t cell_info) const
1372{
1373 if (_dof_transformations_are_identity)
1374 return;
1375
1376 if (_cell_tdim >= 2)
1377 {
1378 // This assumes 3 bits are used per face. This will need updating if
1379 // 3D cells with faces with more than 4 sides are implemented
1380 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1381 int dofstart = 0;
1382 for (auto& edofs0 : _edofs[0])
1383 dofstart += edofs0.size();
1384
1385 // Transform DOFs on edges
1386 {
1387 auto& [v_size_t, matrix] = _etrans.at(cell::type::interval)[0];
1388 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1389 {
1390 // Reverse an edge
1391 if (cell_info >> (face_start + e) & 1)
1392 {
1394 std::span(v_size_t),
1395 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1396 block_size);
1397 }
1398 dofstart += _edofs[1][e].size();
1399 }
1400 }
1401
1402 if (_cell_tdim == 3)
1403 {
1404 // Permute DOFs on faces
1405 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1406 {
1407 auto& trans = _etrans.at(_cell_subentity_types[2][f]);
1408
1409 // Reflect a face
1410 if (cell_info >> (3 * f) & 1)
1411 {
1412 const auto& m = trans[1];
1413 const auto& v_size_t = std::get<0>(m);
1414 const auto& matrix = std::get<1>(m);
1416 std::span(v_size_t),
1417 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1418 block_size);
1419 }
1420
1421 // Rotate a face
1422 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1423 {
1424 const auto& m = trans[0];
1425 const auto& v_size_t = std::get<0>(m);
1426 const auto& matrix = std::get<1>(m);
1428 std::span(v_size_t),
1429 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1430 block_size);
1431 }
1432
1433 dofstart += _edofs[2][f].size();
1434 }
1435 }
1436 }
1437}
1438//-----------------------------------------------------------------------------
1439template <typename T>
1441 const std::span<T>& data, int block_size, std::uint32_t cell_info) const
1442{
1443 if (_dof_transformations_are_identity)
1444 return;
1445
1446 if (_cell_tdim >= 2)
1447 {
1448 // This assumes 3 bits are used per face. This will need updating if
1449 // 3D cells with faces with more than 4 sides are implemented
1450 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1451 int dofstart = 0;
1452 for (auto& edofs0 : _edofs[0])
1453 dofstart += edofs0.size();
1454
1455 // Transform DOFs on edges
1456 {
1457 auto& [v_size_t, matrix] = _etransT.at(cell::type::interval)[0];
1458 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1459 {
1460 // Reverse an edge
1461 if (cell_info >> (face_start + e) & 1)
1462 {
1464 std::span(v_size_t),
1465 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1466 block_size);
1467 }
1468 dofstart += _edofs[1][e].size();
1469 }
1470 }
1471
1472 if (_cell_tdim == 3)
1473 {
1474 // Permute DOFs on faces
1475 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1476 {
1477 auto& trans = _etransT.at(_cell_subentity_types[2][f]);
1478
1479 // Rotate a face
1480 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1481 {
1482 auto& m = trans[0];
1483 auto& v_size_t = std::get<0>(m);
1484 auto& matrix = std::get<1>(m);
1486 std::span(v_size_t),
1487 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1488 block_size);
1489 }
1490
1491 // Reflect a face
1492 if (cell_info >> (3 * f) & 1)
1493 {
1494 auto& m = trans[1];
1495 auto& v_size_t = std::get<0>(m);
1496 auto& matrix = std::get<1>(m);
1498 std::span(v_size_t),
1499 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1500 block_size);
1501 }
1502 dofstart += _edofs[2][f].size();
1503 }
1504 }
1505 }
1506}
1507//-----------------------------------------------------------------------------
1508template <typename T>
1510 const std::span<T>& data, int block_size, std::uint32_t cell_info) const
1511{
1512 if (_dof_transformations_are_identity)
1513 return;
1514
1515 if (_cell_tdim >= 2)
1516 {
1517 // This assumes 3 bits are used per face. This will need updating if
1518 // 3D cells with faces with more than 4 sides are implemented
1519 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1520 int dofstart = 0;
1521 for (auto& edofs0 : _edofs[0])
1522 dofstart += edofs0.size();
1523
1524 // Transform DOFs on edges
1525 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1526 {
1527 // Reverse an edge
1528 if (cell_info >> (face_start + e) & 1)
1529 {
1530 auto& [v_size_t, matrix] = _etrans_invT.at(cell::type::interval)[0];
1531 precompute::apply_matrix(std::span(v_size_t),
1532 cmdspan2_t(matrix.first.data(), matrix.second),
1533 data, dofstart, block_size);
1534 }
1535 dofstart += _edofs[1][e].size();
1536 }
1537
1538 if (_cell_tdim == 3)
1539 {
1540 // Permute DOFs on faces
1541 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1542 {
1543 auto& trans = _etrans_invT.at(_cell_subentity_types[2][f]);
1544
1545 // Reflect a face
1546 if (cell_info >> (3 * f) & 1)
1547 {
1548 auto& m = trans[1];
1549 auto& v_size_t = std::get<0>(m);
1550 auto& matrix = std::get<1>(m);
1552 std::span(v_size_t),
1553 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1554 block_size);
1555 }
1556
1557 // Rotate a face
1558 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1559 {
1560 const auto& m = trans[0];
1561 const auto& v_size_t = std::get<0>(m);
1562 const auto& matrix = std::get<1>(m);
1564 std::span(v_size_t),
1565 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1566 block_size);
1567 }
1568
1569 dofstart += _edofs[2][f].size();
1570 }
1571 }
1572 }
1573}
1574//-----------------------------------------------------------------------------
1575template <typename T>
1577 const std::span<T>& data, int block_size, std::uint32_t cell_info) const
1578{
1579 if (_dof_transformations_are_identity)
1580 return;
1581
1582 if (_cell_tdim >= 2)
1583 {
1584 // This assumes 3 bits are used per face. This will need updating if
1585 // 3D cells with faces with more than 4 sides are implemented
1586 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1587 int dofstart = 0;
1588 for (auto& edofs0 : _edofs[0])
1589 dofstart += edofs0.size();
1590
1591 // Transform DOFs on edges
1592 {
1593 auto& [v_size_t, matrix] = _etrans_inv.at(cell::type::interval)[0];
1594 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1595 {
1596 // Reverse an edge
1597 if (cell_info >> (face_start + e) & 1)
1598 {
1600 std::span(v_size_t),
1601 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1602 block_size);
1603 }
1604 dofstart += _edofs[1][e].size();
1605 }
1606 }
1607
1608 if (_cell_tdim == 3)
1609 {
1610 // Permute DOFs on faces
1611 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1612 {
1613 auto& trans = _etrans_inv.at(_cell_subentity_types[2][f]);
1614
1615 // Rotate a face
1616 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1617 {
1618 auto& m = trans[0];
1619 auto& v_size_t = std::get<0>(m);
1620 auto& matrix = std::get<1>(m);
1622 std::span(v_size_t),
1623 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1624 block_size);
1625 }
1626
1627 // Reflect a face
1628 if (cell_info >> (3 * f) & 1)
1629 {
1630 auto& m = trans[1];
1631 auto& v_size_t = std::get<0>(m);
1632 auto& matrix = std::get<1>(m);
1634 std::span(v_size_t),
1635 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1636 block_size);
1637 }
1638
1639 dofstart += _edofs[2][f].size();
1640 }
1641 }
1642 }
1643}
1644//-----------------------------------------------------------------------------
1645template <typename T>
1647 const std::span<T>& data, int block_size, std::uint32_t cell_info) const
1648{
1649 if (_dof_transformations_are_identity)
1650 return;
1651
1652 if (_cell_tdim >= 2)
1653 {
1654 // This assumes 3 bits are used per face. This will need updating if
1655 // 3D cells with faces with more than 4 sides are implemented
1656 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1657 int dofstart = 0;
1658 for (auto& edofs0 : _edofs[0])
1659 dofstart += edofs0.size();
1660
1661 // Transform DOFs on edges
1662 {
1663 auto& [v_size_t, matrix] = _etrans.at(cell::type::interval)[0];
1664 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1665 {
1666 // Reverse an edge
1667 if (cell_info >> (face_start + e) & 1)
1668 {
1670 std::span(v_size_t),
1671 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1672 block_size);
1673 }
1674
1675 dofstart += _edofs[1][e].size();
1676 }
1677 }
1678
1679 if (_cell_tdim == 3)
1680 {
1681 // Permute DOFs on faces
1682 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1683 {
1684 auto& trans = _etrans.at(_cell_subentity_types[2][f]);
1685
1686 // Reflect a face
1687 if (cell_info >> (3 * f) & 1)
1688 {
1689 auto& m = trans[1];
1690 auto& v_size_t = std::get<0>(m);
1691 auto& matrix = std::get<1>(m);
1693 std::span(v_size_t),
1694 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1695 block_size);
1696 }
1697
1698 // Rotate a face
1699 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1700 {
1701 auto& m = trans[0];
1702 auto& v_size_t = std::get<0>(m);
1703 auto& matrix = std::get<1>(m);
1705 std::span(v_size_t),
1706 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1707 block_size);
1708 }
1709 dofstart += _edofs[2][f].size();
1710 }
1711 }
1712 }
1713}
1714//-----------------------------------------------------------------------------
1715template <typename T>
1717 const std::span<T>& data, int block_size, std::uint32_t cell_info) const
1718{
1719 if (_dof_transformations_are_identity)
1720 return;
1721
1722 if (_cell_tdim >= 2)
1723 {
1724 // This assumes 3 bits are used per face. This will need updating if
1725 // 3D cells with faces with more than 4 sides are implemented
1726 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1727 int dofstart = 0;
1728 for (auto& edofs0 : _edofs[0])
1729 dofstart += edofs0.size();
1730
1731 // Transform DOFs on edges
1732 {
1733 auto& [v_size_t, matrix] = _etrans_invT.at(cell::type::interval)[0];
1734 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1735 {
1736 // Reverse an edge
1737 if (cell_info >> (face_start + e) & 1)
1738 {
1740 std::span(v_size_t),
1741 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1742 block_size);
1743 }
1744 dofstart += _edofs[1][e].size();
1745 }
1746 }
1747
1748 if (_cell_tdim == 3)
1749 {
1750 // Permute DOFs on faces
1751 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1752 {
1753 auto& trans = _etrans_invT.at(_cell_subentity_types[2][f]);
1754
1755 // Reflect a face
1756 if (cell_info >> (3 * f) & 1)
1757 {
1758 auto& m = trans[1];
1759 auto& v_size_t = std::get<0>(m);
1760 auto& matrix = std::get<1>(m);
1762 std::span(v_size_t),
1763 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1764 block_size);
1765 }
1766
1767 // Rotate a face
1768 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1769 {
1770 auto& m = trans[0];
1771 auto& v_size_t = std::get<0>(m);
1772 auto& matrix = std::get<1>(m);
1774 std::span(v_size_t),
1775 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1776 block_size);
1777 }
1778 dofstart += _edofs[2][f].size();
1779 }
1780 }
1781 }
1782}
1783//-----------------------------------------------------------------------------
1784template <typename T>
1786 const std::span<T>& data, int block_size, std::uint32_t cell_info) const
1787{
1788 if (_dof_transformations_are_identity)
1789 return;
1790
1791 if (_cell_tdim >= 2)
1792 {
1793 // This assumes 3 bits are used per face. This will need updating if
1794 // 3D cells with faces with more than 4 sides are implemented
1795 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1796 int dofstart = 0;
1797 for (auto& edofs0 : _edofs[0])
1798 dofstart += edofs0.size();
1799
1800 // Transform DOFs on edges
1801 {
1802 auto& [v_size_t, matrix] = _etransT.at(cell::type::interval)[0];
1803 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1804 {
1805 // Reverse an edge
1806 if (cell_info >> (face_start + e) & 1)
1807 {
1809 std::span(v_size_t),
1810 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1811 block_size);
1812 }
1813 dofstart += _edofs[1][e].size();
1814 }
1815 }
1816
1817 if (_cell_tdim == 3)
1818 {
1819 // Permute DOFs on faces
1820 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1821 {
1822 auto& trans = _etransT.at(_cell_subentity_types[2][f]);
1823
1824 // Rotate a face
1825 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1826 {
1827 auto& m = trans[0];
1828 auto& v_size_t = std::get<0>(m);
1829 auto& matrix = std::get<1>(m);
1831 std::span(v_size_t),
1832 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1833 block_size);
1834 }
1835
1836 // Reflect a face
1837 if (cell_info >> (3 * f) & 1)
1838 {
1839 auto& m = trans[1];
1840 auto& v_size_t = std::get<0>(m);
1841 auto& matrix = std::get<1>(m);
1843 std::span(v_size_t),
1844 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1845 block_size);
1846 }
1847 dofstart += _edofs[2][f].size();
1848 }
1849 }
1850 }
1851}
1852//-----------------------------------------------------------------------------
1853template <typename T>
1855 const std::span<T>& data, int block_size, std::uint32_t cell_info) const
1856{
1857 if (_dof_transformations_are_identity)
1858 return;
1859
1860 if (_cell_tdim >= 2)
1861 {
1862 // This assumes 3 bits are used per face. This will need updating if
1863 // 3D cells with faces with more than 4 sides are implemented
1864 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1865 int dofstart = 0;
1866 for (auto& edofs0 : _edofs[0])
1867 dofstart += edofs0.size();
1868
1869 // Transform DOFs on edges
1870 {
1871 auto& [v_size_t, matrix] = _etrans_inv.at(cell::type::interval)[0];
1872 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1873 {
1874 // Reverse an edge
1875 if (cell_info >> (face_start + e) & 1)
1876 {
1878 std::span(v_size_t),
1879 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1880 block_size);
1881 }
1882 dofstart += _edofs[1][e].size();
1883 }
1884 }
1885
1886 if (_cell_tdim == 3)
1887 {
1888 // Permute DOFs on faces
1889 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1890 {
1891 auto& trans = _etrans_inv.at(_cell_subentity_types[2][f]);
1892
1893 // Rotate a face
1894 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1895 {
1896 auto& m = trans[0];
1897 auto& v_size_t = std::get<0>(m);
1898 auto& matrix = std::get<1>(m);
1900 std::span(v_size_t),
1901 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1902 block_size);
1903 }
1904
1905 // Reflect a face
1906 if (cell_info >> (3 * f) & 1)
1907 {
1908 auto& m = trans[1];
1909 auto& v_size_t = std::get<0>(m);
1910 auto& matrix = std::get<1>(m);
1912 std::span(v_size_t),
1913 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1914 block_size);
1915 }
1916
1917 dofstart += _edofs[2][f].size();
1918 }
1919 }
1920 }
1921}
1922//-----------------------------------------------------------------------------
1923
1924} // namespace basix
A finite element.
Definition: finite-element.h:146
std::map< cell::type, std::pair< std::vector< double >, std::array< std::size_t, 3 > > > entity_transformations() const
Definition: finite-element.cpp:1427
std::pair< std::vector< double >, std::array< std::size_t, 4 > > tabulate(int nd, impl::cmdspan2_t x) const
Compute basis values and derivatives at set of points.
Definition: finite-element.cpp:1049
void apply_inverse_dof_transformation_to_transpose(const std::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1854
cell::type cell_type() const
Definition: finite-element.cpp:1123
sobolev::space sobolev_space() const
Definition: finite-element.cpp:1145
int highest_complete_degree() const
Definition: finite-element.cpp:1129
FiniteElement & operator=(const FiniteElement &element)=default
Assignment operator.
element::lagrange_variant lagrange_variant() const
Definition: finite-element.cpp:1469
FiniteElement(FiniteElement &&element)=default
Move constructor.
bool operator==(const FiniteElement &e) const
Definition: finite-element.cpp:998
const std::array< std::vector< std::pair< std::vector< double >, std::array< std::size_t, 2 > > >, 4 > & x() const
Definition: finite-element.cpp:1446
void apply_dof_transformation(const std::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1369
bool dof_transformations_are_permutations() const
Definition: finite-element.cpp:1149
maps::type map_type() const
Definition: finite-element.cpp:1143
std::vector< std::tuple< std::vector< FiniteElement >, std::vector< int > > > get_tensor_product_representation() const
Definition: finite-element.cpp:1487
const std::vector< std::vector< std::vector< int > > > & entity_dofs() const
Definition: finite-element.cpp:1166
const std::pair< std::vector< double >, std::array< std::size_t, 2 > > & dual_matrix() const
Definition: finite-element.cpp:1433
void apply_inverse_transpose_dof_transformation_to_transpose(const std::span< T > &data, int block_size, std::uint32_t cell_info) const
Apply inverse transpose DOF transformations to some transposed data.
Definition: finite-element.h:1716
int interpolation_nderivs() const
The number of derivatives needed when interpolating.
Definition: finite-element.cpp:1481
FiniteElement & operator=(FiniteElement &&element)=default
Move assignment operator.
const std::vector< std::vector< std::vector< int > > > & entity_closure_dofs() const
Definition: finite-element.cpp:1172
const std::pair< std::vector< double >, std::array< std::size_t, 2 > > & points() const
Definition: finite-element.cpp:1246
int degree() const
Definition: finite-element.cpp:1125
void apply_transpose_dof_transformation(const std::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1440
const std::vector< std::size_t > & value_shape() const
Definition: finite-element.cpp:1134
const std::pair< std::vector< double >, std::array< std::size_t, 2 > > & coefficient_matrix() const
Definition: finite-element.cpp:1459
bool has_tensor_product_factorisation() const
Definition: finite-element.cpp:1464
~FiniteElement()=default
Destructor.
const std::pair< std::vector< double >, std::array< std::size_t, 2 > > & wcoeffs() const
Definition: finite-element.cpp:1439
int dim() const
Definition: finite-element.cpp:1139
int highest_degree() const
Definition: finite-element.cpp:1127
std::pair< std::vector< double >, std::array< std::size_t, 3 > > push_forward(impl::cmdspan3_t U, impl::cmdspan3_t J, std::span< const double > detJ, impl::cmdspan3_t K) const
Definition: finite-element.cpp:1252
std::array< std::size_t, 4 > tabulate_shape(std::size_t nd, std::size_t num_points) const
Definition: finite-element.cpp:1035
element::dpc_variant dpc_variant() const
Definition: finite-element.cpp:1474
void apply_transpose_dof_transformation_to_transpose(const std::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1785
void apply_dof_transformation_to_transpose(const std::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1646
element::family family() const
Definition: finite-element.cpp:1141
void permute_dofs(const std::span< std::int32_t > &dofs, std::uint32_t cell_info) const
Definition: finite-element.cpp:1319
bool dof_transformations_are_identity() const
Definition: finite-element.cpp:1154
void apply_inverse_dof_transformation(const std::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1576
std::pair< std::vector< double >, std::array< std::size_t, 3 > > pull_back(impl::cmdspan3_t u, impl::cmdspan3_t J, std::span< const double > detJ, impl::cmdspan3_t K) const
Definition: finite-element.cpp:1287
FiniteElement(const FiniteElement &element)=default
Copy constructor.
const std::pair< std::vector< double >, std::array< std::size_t, 2 > > & interpolation_matrix() const
Return a matrix of weights interpolation,.
Definition: finite-element.cpp:1160
FiniteElement(element::family family, cell::type cell_type, int degree, const std::vector< std::size_t > &value_shape, const cmdspan2_t &wcoeffs, const std::array< std::vector< cmdspan2_t >, 4 > &x, const std::array< std::vector< cmdspan4_t >, 4 > &M, int interpolation_nderivs, maps::type map_type, sobolev::space sobolev_space, bool discontinuous, int highest_complete_degree, int highest_degree, element::lagrange_variant lvariant, element::dpc_variant dvariant, std::vector< std::tuple< std::vector< FiniteElement >, std::vector< int > > > tensor_factors={})
Construct a finite element.
Definition: finite-element.cpp:606
bool interpolation_is_identity() const
Definition: finite-element.cpp:1476
std::pair< std::vector< double >, std::array< std::size_t, 3 > > base_transformations() const
Get the base transformations.
Definition: finite-element.cpp:1178
const std::array< std::vector< std::pair< std::vector< double >, std::array< std::size_t, 4 > > >, 4 > & M() const
Definition: finite-element.cpp:1453
std::function< void(O &, const P &, const Q &, double, const R &)> map_fn() const
Definition: finite-element.h:641
bool discontinuous() const
Definition: finite-element.cpp:1147
void apply_inverse_transpose_dof_transformation(const std::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1509
void unpermute_dofs(const std::span< std::int32_t > &dofs, std::uint32_t cell_info) const
Definition: finite-element.cpp:1373
type
Cell type.
Definition: cell.h:20
lagrange_variant
Variants of a Lagrange space that can be created.
Definition: element-families.h:12
dpc_variant
Definition: element-families.h:33
impl::cmdspan4_t cmdspan4_t
Typedef for mdspan.
Definition: finite-element.h:113
impl::cmdspan2_t cmdspan2_t
Typedef for mdspan.
Definition: finite-element.h:111
std::tuple< std::array< std::vector< std::vector< double > >, 4 >, std::array< std::vector< std::array< std::size_t, 2 > >, 4 >, std::array< std::vector< std::vector< double > >, 4 >, std::array< std::vector< std::array< std::size_t, 4 > >, 4 > > make_discontinuous(const std::array< std::vector< cmdspan2_t >, 4 > &x, const std::array< std::vector< cmdspan4_t >, 4 > &M, std::size_t tdim, std::size_t value_size)
Definition: finite-element.cpp:394
family
Available element families.
Definition: element-families.h:46
void covariant_piola(O &&r, const P &U, const Q &, double, const R &K)
Covariant Piola map.
Definition: maps.h:39
void contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Contravariant Piola map.
Definition: maps.h:58
void double_contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Double contravariant Piola map.
Definition: maps.h:107
void double_covariant_piola(O &&r, const P &U, const Q &J, double, const R &K)
Double covariant Piola map.
Definition: maps.h:79
type
Map type.
Definition: maps.h:17
void apply_matrix_to_transpose(const std::span< const std::size_t > &v_size_t, const std::experimental::mdspan< const T, std::experimental::dextents< std::size_t, 2 > > &M, const std::span< E > &data, std::size_t offset=0, std::size_t block_size=1)
Apply a (precomputed) matrix to some transposed data.
Definition: precompute.h:244
void apply_matrix(const std::span< const std::size_t > &v_size_t, const std::experimental::mdspan< const T, std::experimental::dextents< std::size_t, 2 > > &M, const std::span< E > &data, std::size_t offset=0, std::size_t block_size=1)
Apply a (precomputed) matrix.
Definition: precompute.h:207
space
Sobolev space type.
Definition: sobolev-spaces.h:13
Basix: FEniCS runtime basis evaluation library.
Definition: cell.h:16
FiniteElement create_element(element::family family, cell::type cell, int degree, element::lagrange_variant lvariant, bool discontinuous)
Definition: finite-element.cpp:345
std::string version()
Definition: finite-element.cpp:1494
FiniteElement create_custom_element(cell::type cell_type, const std::vector< std::size_t > &value_shape, const impl::cmdspan2_t &wcoeffs, const std::array< std::vector< impl::cmdspan2_t >, 4 > &x, const std::array< std::vector< impl::cmdspan4_t >, 4 > &M, int interpolation_nderivs, maps::type map_type, sobolev::space sobolev_space, bool discontinuous, int highest_complete_degree, int highest_degree)
Definition: finite-element.cpp:464