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