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:
319 const std::vector<std::size_t>& value_shape,
321 const std::array<std::vector<mdspan_t<const F, 2>>, 4>& x,
322 const std::array<std::vector<mdspan_t<const F, 4>>, 4>& M,
328 std::vector<int> dof_ordering = {});
329
331 FiniteElement(const FiniteElement& element) = default;
332
334 FiniteElement(FiniteElement&& element) = default;
335
337 ~FiniteElement() = default;
338
340 FiniteElement& operator=(const FiniteElement& element) = default;
341
343 FiniteElement& operator=(FiniteElement&& element) = default;
344
349 bool operator==(const FiniteElement& e) const;
350
352 std::size_t hash() const;
353
363 std::array<std::size_t, 4> tabulate_shape(std::size_t nd,
364 std::size_t num_points) const
365 {
366 std::size_t ndsize = 1;
367 for (std::size_t i = 1; i <= nd; ++i)
368 ndsize *= (_cell_tdim + i);
369 for (std::size_t i = 1; i <= nd; ++i)
370 ndsize /= i;
371 std::size_t vs = std::accumulate(_value_shape.begin(), _value_shape.end(),
372 1, std::multiplies{});
373 std::size_t ndofs = _coeffs.second[0];
374 return {ndsize, num_points, ndofs, vs};
375 }
376
398 std::pair<std::vector<F>, std::array<std::size_t, 4>>
399 tabulate(int nd, impl::mdspan_t<const F, 2> x) const;
400
424 std::pair<std::vector<F>, std::array<std::size_t, 4>>
425 tabulate(int nd, std::span<const F> x,
426 std::array<std::size_t, 2> shape) const;
427
453 void tabulate(int nd, impl::mdspan_t<const F, 2> x,
454 mdspan_t<F, 4> basis) const;
455
481 void tabulate(int nd, std::span<const F> x, std::array<std::size_t, 2> xshape,
482 std::span<F> basis) const;
483
486 cell::type cell_type() const { return _cell_type; }
487
490 polyset::type polyset_type() const { return _poly_type; }
491
494 int degree() const { return _degree; }
495
500 int embedded_superdegree() const { return _embedded_superdegree; }
501
505 int embedded_subdegree() const { return _embedded_subdegree; }
506
510 const std::vector<std::size_t>& value_shape() const { return _value_shape; }
511
515 int dim() const { return _coeffs.second[0]; }
516
519 element::family family() const { return _family; }
520
524 {
525 return _lagrange_variant;
526 }
527
530 element::dpc_variant dpc_variant() const { return _dpc_variant; }
531
534 maps::type map_type() const { return _map_type; }
535
538 sobolev::space sobolev_space() const { return _sobolev_space; }
539
543 bool discontinuous() const { return _discontinuous; }
544
548 {
549 return _dof_transformations_are_permutations;
550 }
551
555 {
556 return _dof_transformations_are_identity;
557 }
558
572 std::pair<std::vector<F>, std::array<std::size_t, 3>>
573 push_forward(impl::mdspan_t<const F, 3> U, impl::mdspan_t<const F, 3> J,
574 std::span<const F> detJ, impl::mdspan_t<const F, 3> K) const;
575
583 std::pair<std::vector<F>, std::array<std::size_t, 3>>
584 pull_back(impl::mdspan_t<const F, 3> u, impl::mdspan_t<const F, 3> J,
585 std::span<const F> detJ, impl::mdspan_t<const F, 3> K) const;
586
618 template <typename O, typename P, typename Q, typename R>
619 std::function<void(O&, const P&, const Q&, F, const R&)> map_fn() const
620 {
621 switch (_map_type)
622 {
623 case maps::type::identity:
624 return [](O& u, const P& U, const Q&, F, const R&)
625 {
626 assert(U.extent(0) == u.extent(0));
627 assert(U.extent(1) == u.extent(1));
628 for (std::size_t i = 0; i < U.extent(0); ++i)
629 for (std::size_t j = 0; j < U.extent(1); ++j)
630 u(i, j) = U(i, j);
631 };
632 case maps::type::covariantPiola:
633 return [](O& u, const P& U, const Q& J, F detJ, const R& K)
634 { maps::covariant_piola(u, U, J, detJ, K); };
635 case maps::type::contravariantPiola:
636 return [](O& u, const P& U, const Q& J, F detJ, const R& K)
638 case maps::type::doubleCovariantPiola:
639 return [](O& u, const P& U, const Q& J, F detJ, const R& K)
641 case maps::type::doubleContravariantPiola:
642 return [](O& u, const P& U, const Q& J, F detJ, const R& K)
644 default:
645 throw std::runtime_error("Map not implemented");
646 }
647 }
648
655 const std::vector<std::vector<std::vector<int>>>& entity_dofs() const
656 {
657 return _edofs;
658 }
659
667 const std::vector<std::vector<std::vector<int>>>& entity_closure_dofs() const
668 {
669 return _e_closure_dofs;
670 }
671
753 std::pair<std::vector<F>, std::array<std::size_t, 3>>
754 base_transformations() const;
755
760 std::map<cell::type, std::pair<std::vector<F>, std::array<std::size_t, 3>>>
762 {
763 return _entity_transformations;
764 }
765
773 void permute_dofs(std::span<std::int32_t> dofs, std::uint32_t cell_info) const
774 {
775 if (!_dof_transformations_are_permutations)
776 {
777 throw std::runtime_error(
778 "The DOF transformations for this element are not permutations");
779 }
780
781 if (_dof_transformations_are_identity)
782 return;
783
785 }
786
794 void unpermute_dofs(std::span<std::int32_t> dofs,
795 std::uint32_t cell_info) const
796 {
797 if (!_dof_transformations_are_permutations)
798 {
799 throw std::runtime_error(
800 "The DOF transformations for this element are not permutations");
801 }
802 if (_dof_transformations_are_identity)
803 return;
804
806 }
807
816 template <typename T>
817 void pre_apply_dof_transformation(std::span<T> data, int block_size,
818 std::uint32_t cell_info) const;
819
828 template <typename T>
830 std::uint32_t cell_info) const;
831
840 template <typename T>
842 std::span<T> data, int block_size, std::uint32_t cell_info) const;
843
852 template <typename T>
854 std::uint32_t cell_info) const;
855
864 template <typename T>
866 int block_size,
867 std::uint32_t cell_info) const;
868
877 template <typename T>
878 void post_apply_dof_transformation(std::span<T> data, int block_size,
879 std::uint32_t cell_info) const;
880
889 template <typename T>
891 std::uint32_t cell_info) const;
892
902 template <typename T>
904 std::span<T> data, int block_size, std::uint32_t cell_info) const;
905
910 const std::pair<std::vector<F>, std::array<std::size_t, 2>>& points() const
911 {
912 return _points;
913 }
914
967 const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
969 {
970 return _matM;
971 }
972
978 const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
980 {
981 return _dual_matrix;
982 }
983
1018 const std::pair<std::vector<F>, std::array<std::size_t, 2>>& wcoeffs() const
1019 {
1020 return _wcoeffs;
1021 }
1022
1027 const std::array<
1028 std::vector<std::pair<std::vector<F>, std::array<std::size_t, 2>>>, 4>&
1029 x() const
1030 {
1031 return _x;
1032 }
1033
1070 const std::array<
1071 std::vector<std::pair<std::vector<F>, std::array<std::size_t, 4>>>, 4>&
1072 M() const
1073 {
1074 return _M;
1075 }
1076
1080 const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
1082 {
1083 return _coeffs;
1084 }
1085
1096 {
1097 return _tensor_factors.size() > 0;
1098 }
1099
1111 std::vector<std::vector<FiniteElement<F>>>
1113 {
1115 throw std::runtime_error("Element has no tensor product representation.");
1116 return _tensor_factors;
1117 }
1118
1121 bool interpolation_is_identity() const { return _interpolation_is_identity; }
1122
1124 int interpolation_nderivs() const { return _interpolation_nderivs; }
1125
1127 const std::vector<int>& dof_ordering() const { return _dof_ordering; }
1128
1129private:
1130 // Data permutation
1131 // @param data Data to be permuted
1132 // @param block_size
1133 // @param cell_info Cell bitmap selecting required permutations
1134 // @param eperm Permutation to use
1135 // @param post Whether reflect is pre- or post- rotation.
1136 template <typename T, bool post>
1137 void permute_data(
1138 std::span<T> data, int block_size, std::uint32_t cell_info,
1139 const std::map<cell::type, std::vector<std::vector<std::size_t>>>& eperm)
1140 const;
1141
1142 using array2_t = std::pair<std::vector<F>, std::array<std::size_t, 2>>;
1143 using array3_t = std::pair<std::vector<F>, std::array<std::size_t, 3>>;
1144 using trans_data_t
1145 = std::vector<std::pair<std::vector<std::size_t>, array2_t>>;
1146
1153 template <typename T, bool post, typename OP>
1154 void
1155 transform_data(std::span<T> data, int block_size, std::uint32_t cell_info,
1156 const std::map<cell::type, trans_data_t>& etrans, OP op) const;
1157
1158 // Cell type
1159 cell::type _cell_type;
1160
1161 // Polyset type
1162 polyset::type _poly_type;
1163
1164 // Topological dimension of the cell
1165 std::size_t _cell_tdim;
1166
1167 // Topological dimension of the cell
1168 std::vector<std::vector<cell::type>> _cell_subentity_types;
1169
1170 // Finite element family
1171 element::family _family;
1172
1173 // Lagrange variant
1174 element::lagrange_variant _lagrange_variant;
1175
1176 // DPC variant
1177 element::dpc_variant _dpc_variant;
1178
1179 // Degree that was input when creating the element
1180 int _degree;
1181
1182 // Degree
1183 int _interpolation_nderivs;
1184
1185 // Highest degree polynomial in element's polyset
1186 int _embedded_superdegree;
1187
1188 // Highest degree space that is a subspace of element's polyset
1189 int _embedded_subdegree;
1190
1191 // Value shape
1192 std::vector<std::size_t> _value_shape;
1193
1195 maps::type _map_type;
1196
1198 sobolev::space _sobolev_space;
1199
1200 // Shape function coefficient of expansion sets on cell. If shape
1201 // function is given by @f$\psi_i = \sum_{k} \phi_{k}
1202 // \alpha^{i}_{k}@f$, then _coeffs(i, j) = @f$\alpha^i_k@f$. ie
1203 // _coeffs.row(i) are the expansion coefficients for shape function i
1204 // (@f$\psi_{i}@f$).
1205 std::pair<std::vector<F>, std::array<std::size_t, 2>> _coeffs;
1206
1207 // Dofs associated with each cell (sub-)entity
1208 std::vector<std::vector<std::vector<int>>> _edofs;
1209
1210 // Dofs associated with the closdure of each cell (sub-)entity
1211 std::vector<std::vector<std::vector<int>>> _e_closure_dofs;
1212
1213 // Entity transformations
1214 std::map<cell::type, array3_t> _entity_transformations;
1215
1216 // Set of points used for point evaluation
1217 // Experimental - currently used for an implementation of
1218 // "tabulate_dof_coordinates" Most useful for Lagrange. This may change or go
1219 // away. For non-Lagrange elements, these points will be used in combination
1220 // with _interpolation_matrix to perform interpolation
1221 std::pair<std::vector<F>, std::array<std::size_t, 2>> _points;
1222
1223 // Interpolation points on the cell. The shape is (entity_dim, num
1224 // entities of given dimension, num_points, tdim)
1225 std::array<std::vector<std::pair<std::vector<F>, std::array<std::size_t, 2>>>,
1226 4>
1227 _x;
1228
1230 std::pair<std::vector<F>, std::array<std::size_t, 2>> _matM;
1231
1232 // Indicates whether or not the DOF transformations are all
1233 // permutations
1234 bool _dof_transformations_are_permutations;
1235
1236 // Indicates whether or not the DOF transformations are all identity
1237 bool _dof_transformations_are_identity;
1238
1239 // The entity permutations (factorised). This will only be set if
1240 // _dof_transformations_are_permutations is True and
1241 // _dof_transformations_are_identity is False
1242 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm;
1243
1244 // The reverse entity permutations (factorised). This will only be set
1245 // if _dof_transformations_are_permutations is True and
1246 // _dof_transformations_are_identity is False
1247 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_rev;
1248
1249 // The entity transformations in precomputed form
1250 std::map<cell::type, trans_data_t> _etrans;
1251
1252 // The transposed entity transformations in precomputed form
1253 std::map<cell::type, trans_data_t> _etransT;
1254
1255 // The inverse entity transformations in precomputed form
1256 std::map<cell::type, trans_data_t> _etrans_inv;
1257
1258 // The inverse transpose entity transformations in precomputed form
1259 std::map<cell::type, trans_data_t> _etrans_invT;
1260
1261 // Indicates whether or not this is the discontinuous version of the
1262 // element
1263 bool _discontinuous;
1264
1265 // The dual matrix
1266 std::pair<std::vector<F>, std::array<std::size_t, 2>> _dual_matrix;
1267
1268 // Dof reordering for different element dof layout compatibility.
1269 // The reference basix layout is ordered by entity, i.e. dofs on
1270 // vertices, followed by edges, faces, then internal dofs.
1271 // _dof_ordering stores the map to the new order required, e.g.
1272 // for a P2 triangle, _dof_ordering=[0 3 5 1 2 4] will place
1273 // dofs 0, 3, 5 on the vertices and 1, 2, 4, on the edges.
1274 std::vector<int> _dof_ordering;
1275
1276 // Tensor product representation
1277 // Entries of tuple are (list of elements on an interval, permutation
1278 // of DOF numbers)
1279 // @todo: For vector-valued elements, a tensor product type and a
1280 // scaling factor may additionally be needed.
1281 std::vector<std::vector<FiniteElement>> _tensor_factors;
1282
1283 // Is the interpolation matrix an identity?
1284 bool _interpolation_is_identity;
1285
1286 // The coefficients that define the polynomial set in terms of the
1287 // orthonormal polynomials
1288 std::pair<std::vector<F>, std::array<std::size_t, 2>> _wcoeffs;
1289
1290 // Interpolation matrices for each entity
1291 using array4_t
1292 = std::vector<std::pair<std::vector<F>, std::array<std::size_t, 4>>>;
1293 std::array<array4_t, 4> _M;
1294 // std::array<
1295 // std::vector<std::pair<std::vector<F>, std::array<std::size_t,
1296 // 4>>>, 4> _M;
1297};
1298
1324template <std::floating_point T>
1325FiniteElement<T> create_custom_element(
1326 cell::type cell_type, const std::vector<std::size_t>& value_shape,
1327 impl::mdspan_t<const T, 2> wcoeffs,
1328 const std::array<std::vector<impl::mdspan_t<const T, 2>>, 4>& x,
1329 const std::array<std::vector<impl::mdspan_t<const T, 4>>, 4>& M,
1330 int interpolation_nderivs, maps::type map_type,
1331 sobolev::space sobolev_space, bool discontinuous, int embedded_subdegree,
1332 int embedded_superdegree, polyset::type poly_type);
1333
1345template <std::floating_point T>
1346FiniteElement<T> create_element(element::family family, cell::type cell,
1347 int degree, element::lagrange_variant lvariant,
1348 element::dpc_variant dvariant,
1349 bool discontinuous,
1350 std::vector<int> dof_ordering = {});
1351
1363std::vector<int> tp_dof_ordering(element::family family, cell::type cell,
1364 int degree, element::lagrange_variant lvariant,
1365 element::dpc_variant dvariant,
1366 bool discontinuous);
1367
1380template <std::floating_point T>
1381std::vector<std::vector<FiniteElement<T>>>
1382tp_factors(element::family family, cell::type cell, int degree,
1384 bool discontinuous, std::vector<int> dof_ordering);
1385
1397template <std::floating_point T>
1398FiniteElement<T>
1399create_tp_element(element::family family, cell::type cell, int degree,
1401 element::dpc_variant dvariant, bool discontinuous);
1402
1405std::string version();
1406
1407//-----------------------------------------------------------------------------
1408template <std::floating_point F>
1409template <typename T, bool post>
1410void FiniteElement<F>::permute_data(
1411 std::span<T> data, int block_size, std::uint32_t cell_info,
1412 const std::map<cell::type, std::vector<std::vector<std::size_t>>>& eperm)
1413 const
1414{
1415 if (_cell_tdim >= 2)
1416 {
1417 // This assumes 3 bits are used per face. This will need updating if 3D
1418 // cells with faces with more than 4 sides are implemented
1419 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1420
1421 // Permute DOFs on edges
1422 {
1423 auto& trans = eperm.at(cell::type::interval)[0];
1424 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1425 {
1426 // Reverse an edge
1427 if (cell_info >> (face_start + e) & 1)
1428 precompute::pre_apply_permutation_mapped(trans, data, _edofs[1][e],
1429 block_size);
1430 }
1431 }
1432
1433 if (_cell_tdim == 3)
1434 {
1435 // Permute DOFs on faces
1436 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1437 {
1438 auto& trans = eperm.at(_cell_subentity_types[2][f]);
1439
1440 // Reflect a face (pre rotate)
1441 if (!post and cell_info >> (3 * f) & 1)
1442 {
1443 precompute::pre_apply_permutation_mapped(trans[1], data, _edofs[2][f],
1444 block_size);
1445 }
1446
1447 // Rotate a face
1448 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1449 {
1450 precompute::pre_apply_permutation_mapped(trans[0], data, _edofs[2][f],
1451 block_size);
1452 }
1453
1454 // Reflect a face (post rotate)
1455 if (post and cell_info >> (3 * f) & 1)
1456 {
1457 precompute::pre_apply_permutation_mapped(trans[1], data, _edofs[2][f],
1458 block_size);
1459 }
1460 }
1461 }
1462 }
1463}
1464//-----------------------------------------------------------------------------
1465template <std::floating_point F>
1466template <typename T, bool post, typename OP>
1467void FiniteElement<F>::transform_data(
1468 std::span<T> data, int block_size, std::uint32_t cell_info,
1469 const std::map<cell::type, trans_data_t>& etrans, OP op) const
1470{
1471 if (_cell_tdim >= 2)
1472 {
1473 // This assumes 3 bits are used per face. This will need updating if
1474 // 3D cells with faces with more than 4 sides are implemented
1475 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1476 int dofstart = 0;
1477 for (auto& edofs0 : _edofs[0])
1478 dofstart += edofs0.size();
1479
1480 // Transform DOFs on edges
1481 {
1482 auto& [v_size_t, matrix] = etrans.at(cell::type::interval)[0];
1483 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1484 {
1485 // Reverse an edge
1486 if (cell_info >> (face_start + e) & 1)
1487 {
1488 op(std::span(v_size_t),
1489 mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
1490 dofstart, block_size);
1491 }
1492 dofstart += _edofs[1][e].size();
1493 }
1494 }
1495
1496 if (_cell_tdim == 3)
1497 {
1498 // Permute DOFs on faces
1499 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1500 {
1501 auto& trans = etrans.at(_cell_subentity_types[2][f]);
1502
1503 // Reflect a face (pre rotation)
1504 if (!post and cell_info >> (3 * f) & 1)
1505 {
1506 const auto& m = trans[1];
1507 const auto& v_size_t = std::get<0>(m);
1508 const auto& matrix = std::get<1>(m);
1509 op(std::span(v_size_t),
1510 mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
1511 dofstart, block_size);
1512 }
1513
1514 // Rotate a face
1515 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1516 {
1517 const auto& m = trans[0];
1518 const auto& v_size_t = std::get<0>(m);
1519 const auto& matrix = std::get<1>(m);
1520 op(std::span(v_size_t),
1521 mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
1522 dofstart, block_size);
1523 }
1524
1525 // Reflect a face (post rotation)
1526 if (post and cell_info >> (3 * f) & 1)
1527 {
1528 const auto& m = trans[1];
1529 const auto& v_size_t = std::get<0>(m);
1530 const auto& matrix = std::get<1>(m);
1531 op(std::span(v_size_t),
1532 mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
1533 dofstart, block_size);
1534 }
1535
1536 dofstart += _edofs[2][f].size();
1537 }
1538 }
1539 }
1540}
1541//-----------------------------------------------------------------------------
1542template <std::floating_point F>
1543template <typename T>
1545 std::span<T> data, int block_size, std::uint32_t cell_info) const
1546{
1547 if (_dof_transformations_are_identity)
1548 return;
1549
1550 if (_dof_transformations_are_permutations)
1551 {
1553 }
1554 else
1555 {
1557 precompute::pre_apply_matrix<F, T>);
1558 }
1559}
1560//-----------------------------------------------------------------------------
1561template <std::floating_point F>
1562template <typename T>
1564 std::span<T> data, int block_size, std::uint32_t cell_info) const
1565{
1566 if (_dof_transformations_are_identity)
1567 return;
1568
1569 if (_dof_transformations_are_permutations)
1570 {
1572 }
1573 else
1574 {
1576 precompute::pre_apply_matrix<F, T>);
1577 }
1578}
1579//-----------------------------------------------------------------------------
1580template <std::floating_point F>
1581template <typename T>
1583 std::span<T> data, int block_size, std::uint32_t cell_info) const
1584{
1585 if (_dof_transformations_are_identity)
1586 return;
1587
1588 if (_dof_transformations_are_permutations)
1589 {
1591 }
1592 else
1593 {
1595 precompute::pre_apply_matrix<F, T>);
1596 }
1597}
1598//-----------------------------------------------------------------------------
1599template <std::floating_point F>
1600template <typename T>
1602 std::span<T> data, int block_size, std::uint32_t cell_info) const
1603{
1604 if (_dof_transformations_are_identity)
1605 return;
1606
1607 if (_dof_transformations_are_permutations)
1608 {
1610 }
1611 else
1612 {
1614 precompute::pre_apply_matrix<F, T>);
1615 }
1616}
1617//-----------------------------------------------------------------------------
1618template <std::floating_point F>
1619template <typename T>
1621 std::span<T> data, int block_size, std::uint32_t cell_info) const
1622{
1623 if (_dof_transformations_are_identity)
1624 return;
1625
1626 if (_dof_transformations_are_permutations)
1627 {
1628 assert(data.size() % block_size == 0);
1629 const int step = data.size() / block_size;
1630 for (int i = 0; i < block_size; ++i)
1631 {
1632 std::span<T> dblock(data.data() + i * step, step);
1634 }
1635 }
1636 else
1637 {
1639 precompute::post_apply_tranpose_matrix<F, T>);
1640 }
1641}
1642//-----------------------------------------------------------------------------
1643template <std::floating_point F>
1644template <typename T>
1646 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 (_dof_transformations_are_permutations)
1652 {
1653 assert(data.size() % block_size == 0);
1654 const int step = data.size() / block_size;
1655 for (int i = 0; i < block_size; ++i)
1656 {
1657 std::span<T> dblock(data.data() + i * step, step);
1659 }
1660 }
1661 else
1662 {
1664 precompute::post_apply_tranpose_matrix<F, T>);
1665 }
1666}
1667//-----------------------------------------------------------------------------
1668template <std::floating_point F>
1669template <typename T>
1671 std::span<T> data, int block_size, std::uint32_t cell_info) const
1672{
1673 if (_dof_transformations_are_identity)
1674 return;
1675
1676 if (_dof_transformations_are_permutations)
1677 {
1678 assert(data.size() % block_size == 0);
1679 const int step = data.size() / block_size;
1680 for (int i = 0; i < block_size; ++i)
1681 {
1682 std::span<T> dblock(data.data() + i * step, step);
1683 permute_data<T, true>(dblock, 1, cell_info, _eperm_rev);
1684 }
1685 }
1686 else
1687 {
1689 precompute::post_apply_tranpose_matrix<F, T>);
1690 }
1691}
1692//-----------------------------------------------------------------------------
1693template <std::floating_point F>
1694template <typename T>
1696 std::span<T> data, int block_size, std::uint32_t cell_info) const
1697{
1698 if (_dof_transformations_are_identity)
1699 return;
1700
1701 if (_dof_transformations_are_permutations)
1702 {
1703 assert(data.size() % block_size == 0);
1704 const int step = data.size() / block_size;
1705 for (int i = 0; i < block_size; ++i)
1706 {
1707 std::span<T> dblock(data.data() + i * step, step);
1708 permute_data<T, true>(dblock, 1, cell_info, _eperm_rev);
1709 }
1710 }
1711 else
1712 {
1714 precompute::post_apply_tranpose_matrix<F, T>);
1715 }
1716}
1717//-----------------------------------------------------------------------------
1718
1719} // namespace basix
A finite element.
Definition finite-element.h:139
void pre_apply_transpose_dof_transformation(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1563
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 embedded_subdegree, int embedded_superdegree, element::lagrange_variant lvariant, element::dpc_variant dvariant, std::vector< int > dof_ordering={})
Construct a finite element.
void post_apply_inverse_transpose_dof_transformation(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1695
const std::array< std::vector< std::pair< std::vector< F >, std::array< std::size_t, 4 > > >, 4 > & M() const
Definition finite-element.h:1072
void pre_apply_inverse_dof_transformation(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1601
std::pair< std::vector< F >, std::array< std::size_t, 3 > > base_transformations() const
Get the base transformations.
Definition finite-element.cpp:1339
void pre_apply_inverse_transpose_dof_transformation(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1582
bool dof_transformations_are_identity() const
Definition finite-element.h:554
bool operator==(const FiniteElement &e) const
Definition finite-element.cpp:1157
std::map< cell::type, std::pair< std::vector< F >, std::array< std::size_t, 3 > > > entity_transformations() const
Definition finite-element.h:761
void pre_apply_dof_transformation(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1544
const std::vector< int > & dof_ordering() const
Get dof layout.
Definition finite-element.h:1127
int embedded_subdegree() const
Definition finite-element.h:505
const std::vector< std::vector< std::vector< int > > > & entity_dofs() const
Definition finite-element.h:655
FiniteElement(FiniteElement &&element)=default
Move constructor.
bool has_tensor_product_factorisation() const
Definition finite-element.h:1095
int interpolation_nderivs() const
The number of derivatives needed when interpolating.
Definition finite-element.h:1124
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:1247
int embedded_superdegree() const
Definition finite-element.h:500
int dim() const
Definition finite-element.h:515
std::size_t hash() const
Get a unique hash of this element.
Definition finite-element.cpp:1197
const std::vector< std::size_t > & value_shape() const
Definition finite-element.h:510
element::lagrange_variant lagrange_variant() const
Definition finite-element.h:523
const std::array< std::vector< std::pair< std::vector< F >, std::array< std::size_t, 2 > > >, 4 > & x() const
Definition finite-element.h:1029
int degree() const
Definition finite-element.h:494
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & points() const
Definition finite-element.h:910
sobolev::space sobolev_space() const
Definition finite-element.h:538
const std::vector< std::vector< std::vector< int > > > & entity_closure_dofs() const
Definition finite-element.h:667
FiniteElement(const FiniteElement &element)=default
Copy constructor.
void post_apply_inverse_dof_transformation(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1645
void post_apply_transpose_dof_transformation(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1620
std::array< std::size_t, 4 > tabulate_shape(std::size_t nd, std::size_t num_points) const
Definition finite-element.h:363
FiniteElement & operator=(FiniteElement &&element)=default
Move assignment operator.
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:968
polyset::type polyset_type() const
Definition finite-element.h:490
void permute_dofs(std::span< std::int32_t > dofs, std::uint32_t cell_info) const
Definition finite-element.h:773
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & wcoeffs() const
Definition finite-element.h:1018
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:1448
bool dof_transformations_are_permutations() const
Definition finite-element.h:547
void post_apply_dof_transformation(std::span< T > data, int block_size, std::uint32_t cell_info) const
Definition finite-element.h:1670
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:1408
cell::type cell_type() const
Definition finite-element.h:486
bool interpolation_is_identity() const
Definition finite-element.h:1121
bool discontinuous() const
Definition finite-element.h:543
std::vector< std::vector< FiniteElement< F > > > get_tensor_product_representation() const
Definition finite-element.h:1112
maps::type map_type() const
Definition finite-element.h:534
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & dual_matrix() const
Definition finite-element.h:979
element::dpc_variant dpc_variant() const
Definition finite-element.h:530
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & coefficient_matrix() const
Definition finite-element.h:1081
FiniteElement & operator=(const FiniteElement &element)=default
Assignment operator.
void unpermute_dofs(std::span< std::int32_t > dofs, std::uint32_t cell_info) const
Definition finite-element.h:794
element::family family() const
Definition finite-element.h:519
~FiniteElement()=default
Destructor.
std::function< void(O &, const P &, const Q &, F, const R &)> map_fn() const
Definition finite-element.h:619
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:520
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:32
family
Available element families.
Definition element-families.h:45
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 pre_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:193
std::vector< int > tp_dof_ordering(element::family family, cell::type cell, int degree, element::lagrange_variant lvariant, element::dpc_variant dvariant, bool discontinuous)
Definition finite-element.cpp:399
std::vector< std::vector< FiniteElement< T > > > tp_factors(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:349
std::string version()
Definition finite-element.cpp:1486
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 embedded_subdegree, int embedded_superdegree, polyset::type poly_type)
Definition finite-element.cpp:609
FiniteElement< T > create_tp_element(element::family family, cell::type cell, int degree, element::lagrange_variant lvariant, element::dpc_variant dvariant, bool discontinuous)
Definition finite-element.cpp:330