Basix
Loading...
Searching...
No Matches
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 "polyset.h"
12#include "precompute.h"
13#include "sobolev-spaces.h"
14#include <array>
15#include <concepts>
16#include <cstdint>
17#include <functional>
18#include <map>
19#include <numeric>
20#include <span>
21#include <string>
22#include <tuple>
23#include <utility>
24#include <vector>
25
27namespace basix
28{
29
30namespace impl
31{
32template <typename T, std::size_t d>
33using mdspan_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
34 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, d>>;
35template <typename T, std::size_t d>
36using mdarray_t
37 = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE::mdarray<
38 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, d>>;
39
42template <typename T>
43std::array<std::vector<mdspan_t<const T, 2>>, 4>
44to_mdspan(std::array<std::vector<mdarray_t<T, 2>>, 4>& x)
45{
46 std::array<std::vector<mdspan_t<const T, 2>>, 4> x1;
47 for (std::size_t i = 0; i < x.size(); ++i)
48 for (std::size_t j = 0; j < x[i].size(); ++j)
49 x1[i].emplace_back(x[i][j].data(), x[i][j].extents());
50
51 return x1;
52}
53
56template <typename T>
57std::array<std::vector<mdspan_t<const T, 4>>, 4>
58to_mdspan(std::array<std::vector<mdarray_t<T, 4>>, 4>& M)
59{
60 std::array<std::vector<mdspan_t<const T, 4>>, 4> M1;
61 for (std::size_t i = 0; i < M.size(); ++i)
62 for (std::size_t j = 0; j < M[i].size(); ++j)
63 M1[i].emplace_back(M[i][j].data(), M[i][j].extents());
64
65 return M1;
66}
67
70template <typename T>
71std::array<std::vector<mdspan_t<const T, 2>>, 4>
72to_mdspan(const std::array<std::vector<std::vector<T>>, 4>& x,
73 const std::array<std::vector<std::array<std::size_t, 2>>, 4>& shape)
74{
75 std::array<std::vector<mdspan_t<const T, 2>>, 4> x1;
76 for (std::size_t i = 0; i < x.size(); ++i)
77 for (std::size_t j = 0; j < x[i].size(); ++j)
78 x1[i].push_back(mdspan_t<const T, 2>(x[i][j].data(), shape[i][j]));
79
80 return x1;
81}
82
85template <typename T>
86std::array<std::vector<mdspan_t<const T, 4>>, 4>
87to_mdspan(const std::array<std::vector<std::vector<T>>, 4>& M,
88 const std::array<std::vector<std::array<std::size_t, 4>>, 4>& shape)
89{
90 std::array<std::vector<mdspan_t<const T, 4>>, 4> M1;
91 for (std::size_t i = 0; i < M.size(); ++i)
92 for (std::size_t j = 0; j < M[i].size(); ++j)
93 M1[i].push_back(mdspan_t<const T, 4>(M[i][j].data(), shape[i][j]));
94
95 return M1;
96}
97
98} // namespace impl
99
100namespace element
101{
103template <typename T, std::size_t d>
104using mdspan_t = impl::mdspan_t<T, d>;
105
121template <std::floating_point T>
122std::tuple<std::array<std::vector<std::vector<T>>, 4>,
123 std::array<std::vector<std::array<std::size_t, 2>>, 4>,
124 std::array<std::vector<std::vector<T>>, 4>,
125 std::array<std::vector<std::array<std::size_t, 4>>, 4>>
126make_discontinuous(const std::array<std::vector<mdspan_t<const T, 2>>, 4>& x,
127 const std::array<std::vector<mdspan_t<const T, 4>>, 4>& M,
128 std::size_t tdim, std::size_t value_size);
129
130} // namespace element
131
137template <std::floating_point F>
139{
140 template <typename T, std::size_t d>
141 using mdspan_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
142 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, d>>;
143
144public:
321 int degree, const std::vector<std::size_t>& value_shape,
323 const std::array<std::vector<mdspan_t<const F, 2>>, 4>& x,
324 const std::array<std::vector<mdspan_t<const F, 4>>, 4>& M,
329 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
331 = {},
332 std::vector<int> dof_ordering = {});
333
335 FiniteElement(const FiniteElement& element) = default;
336
338 FiniteElement(FiniteElement&& element) = default;
339
341 ~FiniteElement() = default;
342
344 FiniteElement& operator=(const FiniteElement& element) = default;
345
347 FiniteElement& operator=(FiniteElement&& element) = default;
348
353 bool operator==(const FiniteElement& e) const;
354
364 std::array<std::size_t, 4> tabulate_shape(std::size_t nd,
365 std::size_t num_points) const
366 {
367 std::size_t ndsize = 1;
368 for (std::size_t i = 1; i <= nd; ++i)
369 ndsize *= (_cell_tdim + i);
370 for (std::size_t i = 1; i <= nd; ++i)
371 ndsize /= i;
372 std::size_t vs = std::accumulate(_value_shape.begin(), _value_shape.end(),
373 1, std::multiplies{});
374 std::size_t ndofs = _coeffs.second[0];
375 return {ndsize, num_points, ndofs, vs};
376 }
377
399 std::pair<std::vector<F>, std::array<std::size_t, 4>>
400 tabulate(int nd, impl::mdspan_t<const F, 2> x) const;
401
425 std::pair<std::vector<F>, std::array<std::size_t, 4>>
426 tabulate(int nd, std::span<const F> x,
427 std::array<std::size_t, 2> shape) const;
428
454 void tabulate(int nd, impl::mdspan_t<const F, 2> x,
455 mdspan_t<F, 4> basis) const;
456
482 void tabulate(int nd, std::span<const F> x, std::array<std::size_t, 2> xshape,
483 std::span<F> basis) const;
484
487 cell::type cell_type() const { return _cell_type; }
488
491 polyset::type polyset_type() const { return _poly_type; }
492
495 int degree() const { return _degree; }
496
501 int highest_degree() const { return _highest_degree; }
502
506 int highest_complete_degree() const { return _highest_complete_degree; }
507
511 const std::vector<std::size_t>& value_shape() const { return _value_shape; }
512
516 int dim() const { return _coeffs.second[0]; }
517
520 element::family family() const { return _family; }
521
525 {
526 return _lagrange_variant;
527 }
528
531 element::dpc_variant dpc_variant() const { return _dpc_variant; }
532
535 maps::type map_type() const { return _map_type; }
536
539 sobolev::space sobolev_space() const { return _sobolev_space; }
540
544 bool discontinuous() const { return _discontinuous; }
545
549 {
550 return _dof_transformations_are_permutations;
551 }
552
556 {
557 return _dof_transformations_are_identity;
558 }
559
573 std::pair<std::vector<F>, std::array<std::size_t, 3>>
574 push_forward(impl::mdspan_t<const F, 3> U, impl::mdspan_t<const F, 3> J,
575 std::span<const F> detJ, impl::mdspan_t<const F, 3> K) const;
576
584 std::pair<std::vector<F>, std::array<std::size_t, 3>>
585 pull_back(impl::mdspan_t<const F, 3> u, impl::mdspan_t<const F, 3> J,
586 std::span<const F> detJ, impl::mdspan_t<const F, 3> K) const;
587
619 template <typename O, typename P, typename Q, typename R>
620 std::function<void(O&, const P&, const Q&, F, const R&)> map_fn() const
621 {
622 switch (_map_type)
623 {
624 case maps::type::identity:
625 return [](O& u, const P& U, const Q&, F, const R&)
626 {
627 assert(U.extent(0) == u.extent(0));
628 assert(U.extent(1) == u.extent(1));
629 for (std::size_t i = 0; i < U.extent(0); ++i)
630 for (std::size_t j = 0; j < U.extent(1); ++j)
631 u(i, j) = U(i, j);
632 };
633 case maps::type::covariantPiola:
634 return [](O& u, const P& U, const Q& J, F detJ, const R& K)
635 { maps::covariant_piola(u, U, J, detJ, K); };
636 case maps::type::contravariantPiola:
637 return [](O& u, const P& U, const Q& J, F detJ, const R& K)
639 case maps::type::doubleCovariantPiola:
640 return [](O& u, const P& U, const Q& J, F detJ, const R& K)
642 case maps::type::doubleContravariantPiola:
643 return [](O& u, const P& U, const Q& J, F detJ, const R& K)
645 default:
646 throw std::runtime_error("Map not implemented");
647 }
648 }
649
656 const std::vector<std::vector<std::vector<int>>>& entity_dofs() const
657 {
658 return _edofs;
659 }
660
668 const std::vector<std::vector<std::vector<int>>>& entity_closure_dofs() const
669 {
670 return _e_closure_dofs;
671 }
672
754 std::pair<std::vector<F>, std::array<std::size_t, 3>>
755 base_transformations() const;
756
761 std::map<cell::type, std::pair<std::vector<F>, std::array<std::size_t, 3>>>
763 {
764 return _entity_transformations;
765 }
766
774 void permute_dofs(std::span<std::int32_t> dofs, std::uint32_t cell_info) const
775 {
776 if (!_dof_transformations_are_permutations)
777 {
778 throw std::runtime_error(
779 "The DOF transformations for this element are not permutations");
780 }
781
782 if (_dof_transformations_are_identity)
783 return;
784
786 }
787
795 void unpermute_dofs(std::span<std::int32_t> dofs,
796 std::uint32_t cell_info) const
797 {
798 if (!_dof_transformations_are_permutations)
799 {
800 throw std::runtime_error(
801 "The DOF transformations for this element are not permutations");
802 }
803 if (_dof_transformations_are_identity)
804 return;
805
807 }
808
817 template <typename T>
818 void apply_dof_transformation(std::span<T> data, int block_size,
819 std::uint32_t cell_info) const;
820
829 template <typename T>
831 std::uint32_t cell_info) const;
832
841 template <typename T>
842 void
844 std::uint32_t cell_info) const;
845
854 template <typename T>
855 void apply_inverse_dof_transformation(std::span<T> data, int block_size,
856 std::uint32_t cell_info) const;
857
866 template <typename T>
868 std::uint32_t cell_info) const;
869
878 template <typename T>
880 std::span<T> data, int block_size, std::uint32_t cell_info) const;
881
891 template <typename T>
893 std::span<T> data, int block_size, std::uint32_t cell_info) const;
894
903 template <typename T>
905 std::span<T> data, int block_size, std::uint32_t cell_info) const;
906
911 const std::pair<std::vector<F>, std::array<std::size_t, 2>>& points() const
912 {
913 return _points;
914 }
915
968 const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
970 {
971 return _matM;
972 }
973
979 const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
981 {
982 return _dual_matrix;
983 }
984
1019 const std::pair<std::vector<F>, std::array<std::size_t, 2>>& wcoeffs() const
1020 {
1021 return _wcoeffs;
1022 }
1023
1028 const std::array<
1029 std::vector<std::pair<std::vector<F>, std::array<std::size_t, 2>>>, 4>&
1030 x() const
1031 {
1032 return _x;
1033 }
1034
1071 const std::array<
1072 std::vector<std::pair<std::vector<F>, std::array<std::size_t, 4>>>, 4>&
1073 M() const
1074 {
1075 return _M;
1076 }
1077
1083 const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
1085 {
1086 return _coeffs;
1087 }
1088
1097 {
1098 return _tensor_factors.size() > 0;
1099 }
1100
1112 std::vector<std::tuple<std::vector<FiniteElement<F>>, std::vector<int>>>
1114 {
1116 throw std::runtime_error("Element has no tensor product representation.");
1117 return _tensor_factors;
1118 }
1119
1122 bool interpolation_is_identity() const { return _interpolation_is_identity; }
1123
1125 int interpolation_nderivs() const { return _interpolation_nderivs; }
1126
1128 const std::vector<int>& dof_ordering() const { return _dof_ordering; }
1129
1130private:
1131 // Data permutation
1132 // @param data Data to be permuted
1133 // @param block_size
1134 // @param cell_info Cell bitmap selecting required permutations
1135 // @param eperm Permutation to use
1136 // @param post Whether reflect is pre- or post- rotation.
1137 template <typename T, bool post>
1138 void permute_data(
1139 std::span<T> data, int block_size, std::uint32_t cell_info,
1140 const std::map<cell::type, std::vector<std::vector<std::size_t>>>& eperm)
1141 const;
1142
1143 using array2_t = std::pair<std::vector<F>, std::array<std::size_t, 2>>;
1144 using array3_t = std::pair<std::vector<F>, std::array<std::size_t, 3>>;
1145 using trans_data_t
1146 = std::vector<std::pair<std::vector<std::size_t>, array2_t>>;
1147
1154 template <typename T, bool post, typename OP>
1155 void
1156 transform_data(std::span<T> data, int block_size, std::uint32_t cell_info,
1157 const std::map<cell::type, trans_data_t>& etrans, OP op) const;
1158
1159 // Cell type
1160 cell::type _cell_type;
1161
1162 // Polyset type
1163 polyset::type _poly_type;
1164
1165 // Topological dimension of the cell
1166 std::size_t _cell_tdim;
1167
1168 // Topological dimension of the cell
1169 std::vector<std::vector<cell::type>> _cell_subentity_types;
1170
1171 // Finite element family
1172 element::family _family;
1173
1174 // Lagrange variant
1175 element::lagrange_variant _lagrange_variant;
1176
1177 // DPC variant
1178 element::dpc_variant _dpc_variant;
1179
1180 // Degree that was input when creating the element
1181 int _degree;
1182
1183 // Degree
1184 int _interpolation_nderivs;
1185
1186 // Highest degree polynomial in element's polyset
1187 int _highest_degree;
1188
1189 // Highest degree space that is a subspace of element's polyset
1190 int _highest_complete_degree;
1191
1192 // Value shape
1193 std::vector<std::size_t> _value_shape;
1194
1196 maps::type _map_type;
1197
1199 sobolev::space _sobolev_space;
1200
1201 // Shape function coefficient of expansion sets on cell. If shape
1202 // function is given by @f$\psi_i = \sum_{k} \phi_{k}
1203 // \alpha^{i}_{k}@f$, then _coeffs(i, j) = @f$\alpha^i_k@f$. ie
1204 // _coeffs.row(i) are the expansion coefficients for shape function i
1205 // (@f$\psi_{i}@f$).
1206 std::pair<std::vector<F>, std::array<std::size_t, 2>> _coeffs;
1207
1208 // Dofs associated with each cell (sub-)entity
1209 std::vector<std::vector<std::vector<int>>> _edofs;
1210
1211 // Dofs associated with the closdure of each cell (sub-)entity
1212 std::vector<std::vector<std::vector<int>>> _e_closure_dofs;
1213
1214 // Entity transformations
1215 std::map<cell::type, array3_t> _entity_transformations;
1216
1217 // Set of points used for point evaluation
1218 // Experimental - currently used for an implementation of
1219 // "tabulate_dof_coordinates" Most useful for Lagrange. This may change or go
1220 // away. For non-Lagrange elements, these points will be used in combination
1221 // with _interpolation_matrix to perform interpolation
1222 std::pair<std::vector<F>, std::array<std::size_t, 2>> _points;
1223
1224 // Interpolation points on the cell. The shape is (entity_dim, num
1225 // entities of given dimension, num_points, tdim)
1226 std::array<std::vector<std::pair<std::vector<F>, std::array<std::size_t, 2>>>,
1227 4>
1228 _x;
1229
1231 std::pair<std::vector<F>, std::array<std::size_t, 2>> _matM;
1232
1233 // Indicates whether or not the DOF transformations are all
1234 // permutations
1235 bool _dof_transformations_are_permutations;
1236
1237 // Indicates whether or not the DOF transformations are all identity
1238 bool _dof_transformations_are_identity;
1239
1240 // The entity permutations (factorised). This will only be set if
1241 // _dof_transformations_are_permutations is True and
1242 // _dof_transformations_are_identity is False
1243 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm;
1244
1245 // The reverse entity permutations (factorised). This will only be set
1246 // if _dof_transformations_are_permutations is True and
1247 // _dof_transformations_are_identity is False
1248 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_rev;
1249
1250 // The entity transformations in precomputed form
1251 std::map<cell::type, trans_data_t> _etrans;
1252
1253 // The transposed entity transformations in precomputed form
1254 std::map<cell::type, trans_data_t> _etransT;
1255
1256 // The inverse entity transformations in precomputed form
1257 std::map<cell::type, trans_data_t> _etrans_inv;
1258
1259 // The inverse transpose entity transformations in precomputed form
1260 std::map<cell::type, trans_data_t> _etrans_invT;
1261
1262 // Indicates whether or not this is the discontinuous version of the
1263 // element
1264 bool _discontinuous;
1265
1266 // The dual matrix
1267 std::pair<std::vector<F>, std::array<std::size_t, 2>> _dual_matrix;
1268
1269 // Tensor product representation
1270 // Entries of tuple are (list of elements on an interval, permutation
1271 // of DOF numbers)
1272 // @todo: For vector-valued elements, a tensor product type and a
1273 // scaling factor may additionally be needed.
1274 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
1275 _tensor_factors;
1276
1277 // Dof reordering for different element dof layout compatibility.
1278 // The reference basix layout is ordered by entity, i.e. dofs on
1279 // vertices, followed by edges, faces, then internal dofs.
1280 // _dof_ordering stores the map to the new order required, e.g.
1281 // for a P2 triangle, _dof_ordering=[0 3 5 1 2 4] will place
1282 // dofs 0, 3, 5 on the vertices and 1, 2, 4, on the edges.
1283 std::vector<int> _dof_ordering;
1284
1285 // Is the interpolation matrix an identity?
1286 bool _interpolation_is_identity;
1287
1288 // The coefficients that define the polynomial set in terms of the
1289 // orthonormal polynomials
1290 std::pair<std::vector<F>, std::array<std::size_t, 2>> _wcoeffs;
1291
1292 // Interpolation matrices for each entity
1293 using array4_t
1294 = std::vector<std::pair<std::vector<F>, std::array<std::size_t, 4>>>;
1295 std::array<array4_t, 4> _M;
1296 // std::array<
1297 // std::vector<std::pair<std::vector<F>, std::array<std::size_t,
1298 // 4>>>, 4> _M;
1299};
1300
1326template <std::floating_point T>
1327FiniteElement<T> create_custom_element(
1328 cell::type cell_type, const std::vector<std::size_t>& value_shape,
1329 impl::mdspan_t<const T, 2> wcoeffs,
1330 const std::array<std::vector<impl::mdspan_t<const T, 2>>, 4>& x,
1331 const std::array<std::vector<impl::mdspan_t<const T, 4>>, 4>& M,
1332 int interpolation_nderivs, maps::type map_type,
1333 sobolev::space sobolev_space, bool discontinuous,
1334 int highest_complete_degree, int highest_degree, polyset::type poly_type);
1335
1347template <std::floating_point T>
1348FiniteElement<T> create_element(element::family family, cell::type cell,
1349 int degree, element::lagrange_variant lvariant,
1350 element::dpc_variant dvariant,
1351 bool discontinuous,
1352 std::vector<int> dof_ordering = {});
1353
1356std::string version();
1357
1358//-----------------------------------------------------------------------------
1359template <std::floating_point F>
1360template <typename T, bool post>
1361void FiniteElement<F>::permute_data(
1362 std::span<T> data, int block_size, std::uint32_t cell_info,
1363 const std::map<cell::type, std::vector<std::vector<std::size_t>>>& eperm)
1364 const
1365{
1366 if (_cell_tdim >= 2)
1367 {
1368 // This assumes 3 bits are used per face. This will need updating if 3D
1369 // cells with faces with more than 4 sides are implemented
1370 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1371
1372 // Permute DOFs on edges
1373 {
1374 auto& trans = eperm.at(cell::type::interval)[0];
1375 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1376 {
1377 // Reverse an edge
1378 if (cell_info >> (face_start + e) & 1)
1379 precompute::apply_permutation_mapped(trans, data, _edofs[1][e],
1380 block_size);
1381 }
1382 }
1383
1384 if (_cell_tdim == 3)
1385 {
1386 // Permute DOFs on faces
1387 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1388 {
1389 auto& trans = eperm.at(_cell_subentity_types[2][f]);
1390
1391 // Reflect a face (pre rotate)
1392 if (!post and cell_info >> (3 * f) & 1)
1393 precompute::apply_permutation_mapped(trans[1], data, _edofs[2][f],
1394 block_size);
1395
1396 // Rotate a face
1397 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1398 precompute::apply_permutation_mapped(trans[0], data, _edofs[2][f],
1399 block_size);
1400
1401 // Reflect a face (post rotate)
1402 if (post and cell_info >> (3 * f) & 1)
1403 precompute::apply_permutation_mapped(trans[1], data, _edofs[2][f],
1404 block_size);
1405 }
1406 }
1407 }
1408}
1409//-----------------------------------------------------------------------------
1410template <std::floating_point F>
1411template <typename T, bool post, typename OP>
1412void FiniteElement<F>::transform_data(
1413 std::span<T> data, int block_size, std::uint32_t cell_info,
1414 const std::map<cell::type, trans_data_t>& etrans, OP op) const
1415{
1416 if (_cell_tdim >= 2)
1417 {
1418 // This assumes 3 bits are used per face. This will need updating if
1419 // 3D cells with faces with more than 4 sides are implemented
1420 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1421 int dofstart = 0;
1422 for (auto& edofs0 : _edofs[0])
1423 dofstart += edofs0.size();
1424
1425 // Transform DOFs on edges
1426 {
1427 auto& [v_size_t, matrix] = etrans.at(cell::type::interval)[0];
1428 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1429 {
1430 // Reverse an edge
1431 if (cell_info >> (face_start + e) & 1)
1432 {
1433 op(std::span(v_size_t),
1434 mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
1435 dofstart, block_size);
1436 }
1437 dofstart += _edofs[1][e].size();
1438 }
1439 }
1440
1441 if (_cell_tdim == 3)
1442 {
1443 // Permute DOFs on faces
1444 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1445 {
1446 auto& trans = etrans.at(_cell_subentity_types[2][f]);
1447
1448 // Reflect a face (pre rotation)
1449 if (!post and cell_info >> (3 * f) & 1)
1450 {
1451 const auto& m = trans[1];
1452 const auto& v_size_t = std::get<0>(m);
1453 const auto& matrix = std::get<1>(m);
1454 op(std::span(v_size_t),
1455 mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
1456 dofstart, block_size);
1457 }
1458
1459 // Rotate a face
1460 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1461 {
1462 const auto& m = trans[0];
1463 const auto& v_size_t = std::get<0>(m);
1464 const auto& matrix = std::get<1>(m);
1465 op(std::span(v_size_t),
1466 mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
1467 dofstart, block_size);
1468 }
1469
1470 // Reflect a face (post rotation)
1471 if (post and cell_info >> (3 * f) & 1)
1472 {
1473 const auto& m = trans[1];
1474 const auto& v_size_t = std::get<0>(m);
1475 const auto& matrix = std::get<1>(m);
1476 op(std::span(v_size_t),
1477 mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
1478 dofstart, block_size);
1479 }
1480
1481 dofstart += _edofs[2][f].size();
1482 }
1483 }
1484 }
1485}
1486//-----------------------------------------------------------------------------
1487template <std::floating_point F>
1488template <typename T>
1490 int block_size,
1491 std::uint32_t cell_info) const
1492{
1493 if (_dof_transformations_are_identity)
1494 return;
1495
1496 if (_dof_transformations_are_permutations)
1498 else
1500 precompute::apply_matrix<F, T>);
1501}
1502//-----------------------------------------------------------------------------
1503template <std::floating_point F>
1504template <typename T>
1506 std::span<T> data, int block_size, std::uint32_t cell_info) const
1507{
1508 if (_dof_transformations_are_identity)
1509 return;
1510
1511 if (_dof_transformations_are_permutations)
1513 else
1515 precompute::apply_matrix<F, T>);
1516}
1517//-----------------------------------------------------------------------------
1518template <std::floating_point F>
1519template <typename T>
1521 std::span<T> data, int block_size, std::uint32_t cell_info) const
1522{
1523 if (_dof_transformations_are_identity)
1524 return;
1525
1526 if (_dof_transformations_are_permutations)
1528 else
1530 precompute::apply_matrix<F, T>);
1531}
1532//-----------------------------------------------------------------------------
1533template <std::floating_point F>
1534template <typename T>
1536 std::span<T> data, int block_size, std::uint32_t cell_info) const
1537{
1538 if (_dof_transformations_are_identity)
1539 return;
1540
1541 if (_dof_transformations_are_permutations)
1543 else
1545 precompute::apply_matrix<F, T>);
1546}
1547//-----------------------------------------------------------------------------
1548template <std::floating_point F>
1549template <typename T>
1551 std::span<T> data, int block_size, std::uint32_t cell_info) const
1552{
1553 if (_dof_transformations_are_identity)
1554 return;
1555
1556 if (_dof_transformations_are_permutations)
1557 {
1558 assert(data.size() % block_size == 0);
1559 const int step = data.size() / block_size;
1560 for (int i = 0; i < block_size; ++i)
1561 {
1562 std::span<T> dblock(data.data() + i * step, step);
1564 }
1565 }
1566 else
1568 precompute::apply_matrix_to_transpose<F, T>);
1569}
1570//-----------------------------------------------------------------------------
1571template <std::floating_point F>
1572template <typename T>
1574 std::span<T> data, int block_size, std::uint32_t cell_info) const
1575{
1576 if (_dof_transformations_are_identity)
1577 return;
1578
1579 if (_dof_transformations_are_permutations)
1580 {
1581 assert(data.size() % block_size == 0);
1582 const int step = data.size() / block_size;
1583 for (int i = 0; i < block_size; ++i)
1584 {
1585 std::span<T> dblock(data.data() + i * step, step);
1587 }
1588 }
1589 else
1591 precompute::apply_matrix_to_transpose<F, T>);
1592}
1593//-----------------------------------------------------------------------------
1594template <std::floating_point F>
1595template <typename T>
1597 std::span<T> data, int block_size, std::uint32_t cell_info) const
1598{
1599 if (_dof_transformations_are_identity)
1600 return;
1601
1602 if (_dof_transformations_are_permutations)
1603 {
1604 assert(data.size() % block_size == 0);
1605 const int step = data.size() / block_size;
1606 for (int i = 0; i < block_size; ++i)
1607 {
1608 std::span<T> dblock(data.data() + i * step, step);
1609 permute_data<T, true>(dblock, 1, cell_info, _eperm_rev);
1610 }
1611 }
1612 else
1614 precompute::apply_matrix_to_transpose<F, T>);
1615}
1616//-----------------------------------------------------------------------------
1617template <std::floating_point F>
1618template <typename T>
1620 std::span<T> data, int block_size, std::uint32_t cell_info) const
1621{
1622 if (_dof_transformations_are_identity)
1623 return;
1624
1625 if (_dof_transformations_are_permutations)
1626 {
1627 assert(data.size() % block_size == 0);
1628 const int step = data.size() / block_size;
1629 for (int i = 0; i < block_size; ++i)
1630 {
1631 std::span<T> dblock(data.data() + i * step, step);
1632 permute_data<T, true>(dblock, 1, cell_info, _eperm_rev);
1633 }
1634 }
1635 else
1636 {
1638 precompute::apply_matrix_to_transpose<F, T>);
1639 }
1640}
1641//-----------------------------------------------------------------------------
1642
1643} // namespace basix
A finite element.
Definition finite-element.h:139
const std::array< std::vector< std::pair< std::vector< F >, std::array< std::size_t, 4 > > >, 4 > & M() const
Definition finite-element.h:1073
void apply_inverse_transpose_dof_transformation_to_transpose(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:1573
void apply_inverse_transpose_dof_transformation(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1520
std::pair< std::vector< F >, std::array< std::size_t, 3 > > base_transformations() const
Get the base transformations.
Definition finite-element.cpp:1088
int highest_complete_degree() const
Definition finite-element.h:506
bool dof_transformations_are_identity() const
Definition finite-element.h:555
bool operator==(const FiniteElement &e) const
Definition finite-element.cpp:959
std::map< cell::type, std::pair< std::vector< F >, std::array< std::size_t, 3 > > > entity_transformations() const
Definition finite-element.h:762
const std::vector< int > & dof_ordering() const
Get dof layout.
Definition finite-element.h:1128
const std::vector< std::vector< std::vector< int > > > & entity_dofs() const
Definition finite-element.h:656
FiniteElement(FiniteElement &&element)=default
Move constructor.
bool has_tensor_product_factorisation() const
Definition finite-element.h:1096
int interpolation_nderivs() const
The number of derivatives needed when interpolating.
Definition finite-element.h:1125
std::pair< std::vector< F >, std::array< std::size_t, 4 > > tabulate(int nd, impl::mdspan_t< const F, 2 > x) const
Compute basis values and derivatives at set of points.
Definition finite-element.cpp:997
int dim() const
Definition finite-element.h:516
void apply_inverse_dof_transformation(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1535
const std::vector< std::size_t > & value_shape() const
Definition finite-element.h:511
element::lagrange_variant lagrange_variant() const
Definition finite-element.h:524
const std::array< std::vector< std::pair< std::vector< F >, std::array< std::size_t, 2 > > >, 4 > & x() const
Definition finite-element.h:1030
int degree() const
Definition finite-element.h:495
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & points() const
Definition finite-element.h:911
sobolev::space sobolev_space() const
Definition finite-element.h:539
const std::vector< std::vector< std::vector< int > > > & entity_closure_dofs() const
Definition finite-element.h:668
FiniteElement(const FiniteElement &element)=default
Copy constructor.
std::vector< std::tuple< std::vector< FiniteElement< F > >, std::vector< int > > > get_tensor_product_representation() const
Definition finite-element.h:1113
std::array< std::size_t, 4 > tabulate_shape(std::size_t nd, std::size_t num_points) const
Definition finite-element.h:364
void apply_dof_transformation_to_transpose(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1550
FiniteElement & operator=(FiniteElement &&element)=default
Move assignment operator.
void apply_inverse_dof_transformation_to_transpose(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1619
void apply_dof_transformation(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1489
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & interpolation_matrix() const
Return a matrix of weights interpolation,.
Definition finite-element.h:969
polyset::type polyset_type() const
Definition finite-element.h:491
void permute_dofs(std::span< std::int32_t > dofs, std::uint32_t cell_info) const
Definition finite-element.h:774
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & wcoeffs() const
Definition finite-element.h:1019
std::pair< std::vector< F >, std::array< std::size_t, 3 > > pull_back(impl::mdspan_t< const F, 3 > u, impl::mdspan_t< const F, 3 > J, std::span< const F > detJ, impl::mdspan_t< const F, 3 > K) const
Definition finite-element.cpp:1197
bool dof_transformations_are_permutations() const
Definition finite-element.h:548
std::pair< std::vector< F >, std::array< std::size_t, 3 > > push_forward(impl::mdspan_t< const F, 3 > U, impl::mdspan_t< const F, 3 > J, std::span< const F > detJ, impl::mdspan_t< const F, 3 > K) const
Definition finite-element.cpp:1157
cell::type cell_type() const
Definition finite-element.h:487
bool interpolation_is_identity() const
Definition finite-element.h:1122
bool discontinuous() const
Definition finite-element.h:544
FiniteElement(element::family family, cell::type cell_type, polyset::type poly_type, int degree, const std::vector< std::size_t > &value_shape, mdspan_t< const F, 2 > wcoeffs, const std::array< std::vector< mdspan_t< const F, 2 > >, 4 > &x, const std::array< std::vector< mdspan_t< const F, 4 > >, 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={}, std::vector< int > dof_ordering={})
Construct a finite element.
maps::type map_type() const
Definition finite-element.h:535
void apply_transpose_dof_transformation_to_transpose(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1596
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & dual_matrix() const
Definition finite-element.h:980
element::dpc_variant dpc_variant() const
Definition finite-element.h:531
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & coefficient_matrix() const
Definition finite-element.h:1084
FiniteElement & operator=(const FiniteElement &element)=default
Assignment operator.
int highest_degree() const
Definition finite-element.h:501
void unpermute_dofs(std::span< std::int32_t > dofs, std::uint32_t cell_info) const
Definition finite-element.h:795
element::family family() const
Definition finite-element.h:520
void apply_transpose_dof_transformation(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1505
~FiniteElement()=default
Destructor.
std::function< void(O &, const P &, const Q &, F, const R &)> map_fn() const
Definition finite-element.h:620
type
Cell type.
Definition cell.h:21
std::tuple< std::array< std::vector< std::vector< T > >, 4 >, std::array< std::vector< std::array< std::size_t, 2 > >, 4 >, std::array< std::vector< std::vector< T > >, 4 >, std::array< std::vector< std::array< std::size_t, 4 > >, 4 > > make_discontinuous(const std::array< std::vector< mdspan_t< const T, 2 > >, 4 > &x, const std::array< std::vector< mdspan_t< const T, 4 > >, 4 > &M, std::size_t tdim, std::size_t value_size)
Definition finite-element.cpp:329
lagrange_variant
Variants of a Lagrange space that can be created.
Definition element-families.h:12
impl::mdspan_t< T, d > mdspan_t
Typedef for mdspan.
Definition finite-element.h:104
dpc_variant
Definition element-families.h:33
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:60
void contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Contravariant Piola map.
Definition maps.h:80
void double_contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Double contravariant Piola map.
Definition maps.h:134
void double_covariant_piola(O &&r, const P &U, const Q &J, double, const R &K)
Double covariant Piola map.
Definition maps.h:102
type
Map type.
Definition maps.h:38
type
Cell type.
Definition polyset.h:136
void apply_permutation_mapped(std::span< const std::size_t > perm, std::span< E > data, std::span< const int > emap, std::size_t block_size=1)
Permutation of mapped data.
Definition precompute.h:156
space
Sobolev space type.
Definition sobolev-spaces.h:13
Basix: FEniCS runtime basis evaluation library.
Definition cell.h:17
FiniteElement< T > create_element(element::family family, cell::type cell, int degree, element::lagrange_variant lvariant, element::dpc_variant dvariant, bool discontinuous, std::vector< int > dof_ordering={})
Definition finite-element.cpp:189
std::string version()
Definition finite-element.cpp:1235
FiniteElement< T > create_custom_element(cell::type cell_type, const std::vector< std::size_t > &value_shape, impl::mdspan_t< const T, 2 > wcoeffs, const std::array< std::vector< impl::mdspan_t< const T, 2 > >, 4 > &x, const std::array< std::vector< impl::mdspan_t< const T, 4 > >, 4 > &M, int interpolation_nderivs, maps::type map_type, sobolev::space sobolev_space, bool discontinuous, int highest_complete_degree, int highest_degree, polyset::type poly_type)
Definition finite-element.cpp:418