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 <array>
13#include <functional>
14#include <map>
15#include <numeric>
16#include <span>
17#include <string>
18#include <tuple>
19#include <vector>
20
22namespace basix
23{
24
25namespace impl
26{
27using mdspan2_t
28 = std::experimental::mdspan<double,
29 std::experimental::dextents<std::size_t, 2>>;
30using mdspan4_t
31 = std::experimental::mdspan<double,
32 std::experimental::dextents<std::size_t, 4>>;
33using cmdspan2_t
34 = std::experimental::mdspan<const double,
35 std::experimental::dextents<std::size_t, 2>>;
36using cmdspan3_t
37 = std::experimental::mdspan<const double,
38 std::experimental::dextents<std::size_t, 3>>;
39using cmdspan4_t
40 = std::experimental::mdspan<const double,
41 std::experimental::dextents<std::size_t, 4>>;
42
43using mdarray2_t
44 = std::experimental::mdarray<double,
45 std::experimental::dextents<std::size_t, 2>>;
46using mdarray4_t
47 = std::experimental::mdarray<double,
48 std::experimental::dextents<std::size_t, 4>>;
49
52inline std::array<std::vector<cmdspan2_t>, 4>
53to_mdspan(std::array<std::vector<mdarray2_t>, 4>& x)
54{
55 std::array<std::vector<cmdspan2_t>, 4> x1;
56 for (std::size_t i = 0; i < x.size(); ++i)
57 for (std::size_t j = 0; j < x[i].size(); ++j)
58 x1[i].emplace_back(x[i][j].data(), x[i][j].extents());
59
60 return x1;
61}
62
65inline std::array<std::vector<cmdspan4_t>, 4>
66to_mdspan(std::array<std::vector<mdarray4_t>, 4>& M)
67{
68 std::array<std::vector<cmdspan4_t>, 4> M1;
69 for (std::size_t i = 0; i < M.size(); ++i)
70 for (std::size_t j = 0; j < M[i].size(); ++j)
71 M1[i].emplace_back(M[i][j].data(), M[i][j].extents());
72
73 return M1;
74}
75
78inline std::array<std::vector<cmdspan2_t>, 4>
79to_mdspan(const std::array<std::vector<std::vector<double>>, 4>& x,
80 const std::array<std::vector<std::array<std::size_t, 2>>, 4>& shape)
81{
82 std::array<std::vector<cmdspan2_t>, 4> x1;
83 for (std::size_t i = 0; i < x.size(); ++i)
84 for (std::size_t j = 0; j < x[i].size(); ++j)
85 x1[i].push_back(cmdspan2_t(x[i][j].data(), shape[i][j]));
86
87 return x1;
88}
89
92inline std::array<std::vector<cmdspan4_t>, 4>
93to_mdspan(const std::array<std::vector<std::vector<double>>, 4>& M,
94 const std::array<std::vector<std::array<std::size_t, 4>>, 4>& shape)
95{
96 std::array<std::vector<cmdspan4_t>, 4> M1;
97 for (std::size_t i = 0; i < M.size(); ++i)
98 for (std::size_t j = 0; j < M[i].size(); ++j)
99 M1[i].push_back(cmdspan4_t(M[i][j].data(), shape[i][j]));
100
101 return M1;
102}
103
104} // namespace impl
105
106namespace element
107{
109using cmdspan2_t = impl::cmdspan2_t;
111using cmdspan4_t = impl::cmdspan4_t;
112
128std::tuple<std::array<std::vector<std::vector<double>>, 4>,
129 std::array<std::vector<std::array<std::size_t, 2>>, 4>,
130 std::array<std::vector<std::vector<double>>, 4>,
131 std::array<std::vector<std::array<std::size_t, 4>>, 4>>
132make_discontinuous(const std::array<std::vector<cmdspan2_t>, 4>& x,
133 const std::array<std::vector<cmdspan4_t>, 4>& M,
134 std::size_t tdim, std::size_t value_size);
135
136} // namespace element
137
144{
145 using cmdspan2_t
146 = std::experimental::mdspan<const double,
147 std::experimental::dextents<std::size_t, 2>>;
148 using cmdspan4_t
149 = std::experimental::mdspan<const double,
150 std::experimental::dextents<std::size_t, 4>>;
151
152public:
324 const std::vector<std::size_t>& value_shape, const cmdspan2_t& wcoeffs,
325 const std::array<std::vector<cmdspan2_t>, 4>& x,
326 const std::array<std::vector<cmdspan4_t>, 4>& M,
330 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
331 tensor_factors
332 = {});
333
337 const std::vector<std::size_t>& value_shape, const cmdspan2_t& wcoeffs,
338 const std::array<std::vector<cmdspan2_t>, 4>& x,
339 const std::array<std::vector<cmdspan4_t>, 4>& M,
343 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
344 tensor_factors
345 = {});
346
350 const std::vector<std::size_t>& value_shape, const cmdspan2_t& wcoeffs,
351 const std::array<std::vector<cmdspan2_t>, 4>& x,
352 const std::array<std::vector<cmdspan4_t>, 4>& M,
355 element::dpc_variant dvariant,
356 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
357 tensor_factors
358 = {});
359
363 const std::vector<std::size_t>& value_shape, const cmdspan2_t& wcoeffs,
364 const std::array<std::vector<cmdspan2_t>, 4>& x,
365 const std::array<std::vector<cmdspan4_t>, 4>& M,
366 int interpolation_nderivs, maps::type map_type, bool disccontinuous,
368 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
369 tensor_factors
370 = {});
371
373 FiniteElement(const FiniteElement& element) = default;
374
376 FiniteElement(FiniteElement&& element) = default;
377
379 ~FiniteElement() = default;
380
382 FiniteElement& operator=(const FiniteElement& element) = default;
383
385 FiniteElement& operator=(FiniteElement&& element) = default;
386
391 bool operator==(const FiniteElement& e) const;
392
402 std::array<std::size_t, 4> tabulate_shape(std::size_t nd,
403 std::size_t num_points) const;
404
426 std::pair<std::vector<double>, std::array<std::size_t, 4>>
427 tabulate(int nd, impl::cmdspan2_t x) const;
428
452 std::pair<std::vector<double>, std::array<std::size_t, 4>>
453 tabulate(int nd, const std::span<const double>& x,
454 std::array<std::size_t, 2> shape) const;
455
481 void tabulate(int nd, impl::cmdspan2_t x, impl::mdspan4_t basis) const;
482
508 void tabulate(int nd, const std::span<const double>& x,
509 std::array<std::size_t, 2> xshape,
510 const std::span<double>& basis) const;
511
514 cell::type cell_type() const;
515
518 int degree() const;
519
524 int highest_degree() const;
525
529 int highest_complete_degree() const;
530
534 const std::vector<std::size_t>& value_shape() const;
535
539 int dim() const;
540
543 element::family family() const;
544
548
552
555 maps::type map_type() const;
556
560 bool discontinuous() const;
561
565
569
583 std::pair<std::vector<double>, std::array<std::size_t, 3>>
584 push_forward(impl::cmdspan3_t U, impl::cmdspan3_t J,
585 std::span<const double> detJ, impl::cmdspan3_t K) const;
586
594 std::pair<std::vector<double>, std::array<std::size_t, 3>>
595 pull_back(impl::cmdspan3_t u, impl::cmdspan3_t J,
596 std::span<const double> detJ, impl::cmdspan3_t K) const;
597
629 template <typename O, typename P, typename Q, typename R>
630 std::function<void(O&, const P&, const Q&, double, const R&)> map_fn() const
631 {
632 switch (_map_type)
633 {
634 case maps::type::identity:
635 return [](O& u, const P& U, const Q&, double, const R&)
636 {
637 assert(U.extent(0) == u.extent(0));
638 assert(U.extent(1) == u.extent(1));
639 for (std::size_t i = 0; i < U.extent(0); ++i)
640 for (std::size_t j = 0; j < U.extent(1); ++j)
641 u(i, j) = U(i, j);
642 };
643 case maps::type::covariantPiola:
644 return [](O& u, const P& U, const Q& J, double detJ, const R& K)
645 { maps::covariant_piola(u, U, J, detJ, K); };
646 case maps::type::contravariantPiola:
647 return [](O& u, const P& U, const Q& J, double detJ, const R& K)
648 { maps::contravariant_piola(u, U, J, detJ, K); };
649 case maps::type::doubleCovariantPiola:
650 return [](O& u, const P& U, const Q& J, double detJ, const R& K)
651 { maps::double_covariant_piola(u, U, J, detJ, K); };
652 case maps::type::doubleContravariantPiola:
653 return [](O& u, const P& U, const Q& J, double detJ, const R& K)
654 { maps::double_contravariant_piola(u, U, J, detJ, K); };
655 default:
656 throw std::runtime_error("Map not implemented");
657 }
658 }
659
666 const std::vector<std::vector<std::vector<int>>>& entity_dofs() const;
667
675 const std::vector<std::vector<std::vector<int>>>& entity_closure_dofs() const;
676
758 std::pair<std::vector<double>, std::array<std::size_t, 3>>
759 base_transformations() const;
760
765 std::map<cell::type,
766 std::pair<std::vector<double>, std::array<std::size_t, 3>>>
768
776 void permute_dofs(const std::span<std::int32_t>& dofs,
777 std::uint32_t cell_info) const;
778
786 void unpermute_dofs(const std::span<std::int32_t>& dofs,
787 std::uint32_t cell_info) const;
788
797 template <typename T>
798 void apply_dof_transformation(const std::span<T>& data, int block_size,
799 std::uint32_t cell_info) const;
800
809 template <typename T>
810 void apply_transpose_dof_transformation(const std::span<T>& data,
811 int block_size,
812 std::uint32_t cell_info) const;
813
822 template <typename T>
824 const std::span<T>& data, int block_size, std::uint32_t cell_info) const;
825
834 template <typename T>
835 void apply_inverse_dof_transformation(const std::span<T>& data,
836 int block_size,
837 std::uint32_t cell_info) const;
838
847 template <typename T>
848 void apply_dof_transformation_to_transpose(const std::span<T>& data,
849 int block_size,
850 std::uint32_t cell_info) const;
851
860 template <typename T>
862 const std::span<T>& data, int block_size, std::uint32_t cell_info) const;
863
873 template <typename T>
875 const std::span<T>& data, int block_size, std::uint32_t cell_info) const;
876
885 template <typename T>
887 const std::span<T>& data, int block_size, std::uint32_t cell_info) const;
888
893 const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
894 points() const;
895
948 const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
949 interpolation_matrix() const;
950
956 const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
957 dual_matrix() const;
958
993 const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
994 wcoeffs() const;
995
1000 const std::array<
1001 std::vector<std::pair<std::vector<double>, std::array<std::size_t, 2>>>,
1002 4>&
1003 x() const;
1004
1041 const std::array<
1042 std::vector<std::pair<std::vector<double>, std::array<std::size_t, 4>>>,
1043 4>&
1044 M() const;
1045
1051 const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
1052 coefficient_matrix() const;
1053
1062
1074 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
1076
1079 bool interpolation_is_identity() const;
1080
1082 int interpolation_nderivs() const;
1083
1084private:
1085 // Cell type
1086 cell::type _cell_type;
1087
1088 // Topological dimension of the cell
1089 std::size_t _cell_tdim;
1090
1091 // Topological dimension of the cell
1092 std::vector<std::vector<cell::type>> _cell_subentity_types;
1093
1094 // Finite element family
1095 element::family _family;
1096
1097 // Lagrange variant
1098 element::lagrange_variant _lagrange_variant;
1099
1100 // Lagrange variant
1101 element::dpc_variant _dpc_variant;
1102
1103 // Degree that was input when creating the element
1104 int _degree;
1105
1106 // Degree
1107 int _interpolation_nderivs;
1108
1109 // Highest degree polynomial in element's polyset
1110 int _highest_degree;
1111
1112 // Highest degree space that is a subspace of element's polyset
1113 int _highest_complete_degree;
1114
1115 // Value shape
1116 std::vector<std::size_t> _value_shape;
1117
1119 maps::type _map_type;
1120
1121 // Shape function coefficient of expansion sets on cell. If shape
1122 // function is given by @f$\psi_i = \sum_{k} \phi_{k}
1123 // \alpha^{i}_{k}@f$, then _coeffs(i, j) = @f$\alpha^i_k@f$. ie
1124 // _coeffs.row(i) are the expansion coefficients for shape function i
1125 // (@f$\psi_{i}@f$).
1126 std::pair<std::vector<double>, std::array<std::size_t, 2>> _coeffs;
1127
1128 // Dofs associated with each cell (sub-)entity
1129 std::vector<std::vector<std::vector<int>>> _edofs;
1130
1131 // Dofs associated with the closdure of each cell (sub-)entity
1132 std::vector<std::vector<std::vector<int>>> _e_closure_dofs;
1133
1134 using array2_t = std::pair<std::vector<double>, std::array<std::size_t, 2>>;
1135 using array3_t = std::pair<std::vector<double>, std::array<std::size_t, 3>>;
1136
1137 // Entity transformations
1138 std::map<cell::type, array3_t> _entity_transformations;
1139
1140 // Set of points used for point evaluation
1141 // Experimental - currently used for an implementation of
1142 // "tabulate_dof_coordinates" Most useful for Lagrange. This may change or go
1143 // away. For non-Lagrange elements, these points will be used in combination
1144 // with _interpolation_matrix to perform interpolation
1145 std::pair<std::vector<double>, std::array<std::size_t, 2>> _points;
1146
1147 // Interpolation points on the cell. The shape is (entity_dim, num
1148 // entities of given dimension, num_points, tdim)
1149 std::array<
1150 std::vector<std::pair<std::vector<double>, std::array<std::size_t, 2>>>,
1151 4>
1152 _x;
1153
1155 std::pair<std::vector<double>, std::array<std::size_t, 2>> _matM;
1156
1157 // Indicates whether or not the DOF transformations are all
1158 // permutations
1159 bool _dof_transformations_are_permutations;
1160
1161 // Indicates whether or not the DOF transformations are all identity
1162 bool _dof_transformations_are_identity;
1163
1164 // The entity permutations (factorised). This will only be set if
1165 // _dof_transformations_are_permutations is True and
1166 // _dof_transformations_are_identity is False
1167 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm;
1168
1169 // The reverse entity permutations (factorised). This will only be set
1170 // if _dof_transformations_are_permutations is True and
1171 // _dof_transformations_are_identity is False
1172 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_rev;
1173
1174 using trans_data_t
1175 = std::vector<std::pair<std::vector<std::size_t>, array2_t>>;
1176
1177 // The entity transformations in precomputed form
1178 std::map<cell::type, trans_data_t> _etrans;
1179
1180 // The transposed entity transformations in precomputed form
1181 std::map<cell::type, trans_data_t> _etransT;
1182
1183 // The inverse entity transformations in precomputed form
1184 std::map<cell::type, trans_data_t> _etrans_inv;
1185
1186 // The inverse transpose entity transformations in precomputed form
1187 std::map<cell::type, trans_data_t> _etrans_invT;
1188
1189 // Indicates whether or not this is the discontinuous version of the
1190 // element
1191 bool _discontinuous;
1192
1193 // The dual matrix
1194 std::pair<std::vector<double>, std::array<std::size_t, 2>> _dual_matrix;
1195
1196 // Tensor product representation
1197 // Entries of tuple are (list of elements on an interval, permutation
1198 // of DOF numbers)
1199 // @todo: For vector-valued elements, a tensor product type and a
1200 // scaling factor may additionally be needed.
1201 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
1202 _tensor_factors;
1203
1204 // Is the interpolation matrix an identity?
1205 bool _interpolation_is_identity;
1206
1207 // The coefficients that define the polynomial set in terms of the
1208 // orthonormal polynomials
1209 std::pair<std::vector<double>, std::array<std::size_t, 2>> _wcoeffs;
1210
1211 // Interpolation matrices for each entity
1212 using array4_t
1213 = std::vector<std::pair<std::vector<double>, std::array<std::size_t, 4>>>;
1214 std::array<array4_t, 4> _M;
1215 // std::array<
1216 // std::vector<std::pair<std::vector<double>, std::array<std::size_t,
1217 // 4>>>, 4> _M;
1218};
1219
1243FiniteElement create_custom_element(
1244 cell::type cell_type, const std::vector<std::size_t>& value_shape,
1245 const impl::cmdspan2_t& wcoeffs,
1246 const std::array<std::vector<impl::cmdspan2_t>, 4>& x,
1247 const std::array<std::vector<impl::cmdspan4_t>, 4>& M,
1248 int interpolation_nderivs, maps::type map_type, bool discontinuous,
1249 int highest_complete_degree, int highest_degree);
1250
1260FiniteElement create_element(element::family family, cell::type cell,
1261 int degree, element::lagrange_variant lvariant,
1262 bool discontinuous);
1263
1274FiniteElement create_element(element::family family, cell::type cell,
1275 int degree, element::lagrange_variant lvariant,
1276 element::dpc_variant dvariant, bool discontinuous);
1277
1287FiniteElement create_element(element::family family, cell::type cell,
1288 int degree, element::dpc_variant dvariant,
1289 bool discontinuous);
1290
1301FiniteElement create_element(element::family family, cell::type cell,
1302 int degree, bool discontinuous);
1303
1311FiniteElement create_element(element::family family, cell::type cell,
1312 int degree, element::lagrange_variant lvariant);
1313
1323FiniteElement create_element(element::family family, cell::type cell,
1324 int degree, element::lagrange_variant lvariant,
1325 element::dpc_variant dvariant);
1326
1334FiniteElement create_element(element::family family, cell::type cell,
1335 int degree, element::dpc_variant dvariant);
1336
1343FiniteElement create_element(element::family family, cell::type cell,
1344 int degree);
1345
1348std::string version();
1349
1350//-----------------------------------------------------------------------------
1351template <typename T>
1352void FiniteElement::apply_dof_transformation(const std::span<T>& data,
1353 int block_size,
1354 std::uint32_t cell_info) const
1355{
1356 if (_dof_transformations_are_identity)
1357 return;
1358
1359 if (_cell_tdim >= 2)
1360 {
1361 // This assumes 3 bits are used per face. This will need updating if
1362 // 3D cells with faces with more than 4 sides are implemented
1363 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1364 int dofstart = 0;
1365 for (auto& edofs0 : _edofs[0])
1366 dofstart += edofs0.size();
1367
1368 // Transform DOFs on edges
1369 {
1370 auto& [v_size_t, matrix] = _etrans.at(cell::type::interval)[0];
1371 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1372 {
1373 // Reverse an edge
1374 if (cell_info >> (face_start + e) & 1)
1375 {
1377 std::span(v_size_t),
1378 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1379 block_size);
1380 }
1381 dofstart += _edofs[1][e].size();
1382 }
1383 }
1384
1385 if (_cell_tdim == 3)
1386 {
1387 // Permute DOFs on faces
1388 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1389 {
1390 auto& trans = _etrans.at(_cell_subentity_types[2][f]);
1391
1392 // Reflect a face
1393 if (cell_info >> (3 * f) & 1)
1394 {
1395 const auto& m = trans[1];
1396 const auto& v_size_t = std::get<0>(m);
1397 const auto& matrix = std::get<1>(m);
1399 std::span(v_size_t),
1400 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1401 block_size);
1402 }
1403
1404 // Rotate a face
1405 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1406 {
1407 const auto& m = trans[0];
1408 const auto& v_size_t = std::get<0>(m);
1409 const auto& matrix = std::get<1>(m);
1411 std::span(v_size_t),
1412 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1413 block_size);
1414 }
1415
1416 dofstart += _edofs[2][f].size();
1417 }
1418 }
1419 }
1420}
1421//-----------------------------------------------------------------------------
1422template <typename T>
1424 const std::span<T>& data, int block_size, std::uint32_t cell_info) const
1425{
1426 if (_dof_transformations_are_identity)
1427 return;
1428
1429 if (_cell_tdim >= 2)
1430 {
1431 // This assumes 3 bits are used per face. This will need updating if
1432 // 3D cells with faces with more than 4 sides are implemented
1433 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1434 int dofstart = 0;
1435 for (auto& edofs0 : _edofs[0])
1436 dofstart += edofs0.size();
1437
1438 // Transform DOFs on edges
1439 {
1440 auto& [v_size_t, matrix] = _etransT.at(cell::type::interval)[0];
1441 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1442 {
1443 // Reverse an edge
1444 if (cell_info >> (face_start + e) & 1)
1445 {
1447 std::span(v_size_t),
1448 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1449 block_size);
1450 }
1451 dofstart += _edofs[1][e].size();
1452 }
1453 }
1454
1455 if (_cell_tdim == 3)
1456 {
1457 // Permute DOFs on faces
1458 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1459 {
1460 auto& trans = _etransT.at(_cell_subentity_types[2][f]);
1461
1462 // Rotate a face
1463 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1464 {
1465 auto& m = trans[0];
1466 auto& v_size_t = std::get<0>(m);
1467 auto& matrix = std::get<1>(m);
1469 std::span(v_size_t),
1470 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1471 block_size);
1472 }
1473
1474 // Reflect a face
1475 if (cell_info >> (3 * f) & 1)
1476 {
1477 auto& m = trans[1];
1478 auto& v_size_t = std::get<0>(m);
1479 auto& matrix = std::get<1>(m);
1481 std::span(v_size_t),
1482 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1483 block_size);
1484 }
1485 dofstart += _edofs[2][f].size();
1486 }
1487 }
1488 }
1489}
1490//-----------------------------------------------------------------------------
1491template <typename T>
1493 const std::span<T>& data, int block_size, std::uint32_t cell_info) const
1494{
1495 if (_dof_transformations_are_identity)
1496 return;
1497
1498 if (_cell_tdim >= 2)
1499 {
1500 // This assumes 3 bits are used per face. This will need updating if
1501 // 3D cells with faces with more than 4 sides are implemented
1502 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1503 int dofstart = 0;
1504 for (auto& edofs0 : _edofs[0])
1505 dofstart += edofs0.size();
1506
1507 // Transform DOFs on edges
1508 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1509 {
1510 // Reverse an edge
1511 if (cell_info >> (face_start + e) & 1)
1512 {
1513 auto& [v_size_t, matrix] = _etrans_invT.at(cell::type::interval)[0];
1514 precompute::apply_matrix(std::span(v_size_t),
1515 cmdspan2_t(matrix.first.data(), matrix.second),
1516 data, dofstart, block_size);
1517 }
1518 dofstart += _edofs[1][e].size();
1519 }
1520
1521 if (_cell_tdim == 3)
1522 {
1523 // Permute DOFs on faces
1524 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1525 {
1526 auto& trans = _etrans_invT.at(_cell_subentity_types[2][f]);
1527
1528 // Reflect a face
1529 if (cell_info >> (3 * f) & 1)
1530 {
1531 auto& m = trans[1];
1532 auto& v_size_t = std::get<0>(m);
1533 auto& matrix = std::get<1>(m);
1535 std::span(v_size_t),
1536 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1537 block_size);
1538 }
1539
1540 // Rotate a face
1541 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1542 {
1543 const auto& m = trans[0];
1544 const auto& v_size_t = std::get<0>(m);
1545 const auto& matrix = std::get<1>(m);
1547 std::span(v_size_t),
1548 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1549 block_size);
1550 }
1551
1552 dofstart += _edofs[2][f].size();
1553 }
1554 }
1555 }
1556}
1557//-----------------------------------------------------------------------------
1558template <typename T>
1560 const std::span<T>& data, int block_size, std::uint32_t cell_info) const
1561{
1562 if (_dof_transformations_are_identity)
1563 return;
1564
1565 if (_cell_tdim >= 2)
1566 {
1567 // This assumes 3 bits are used per face. This will need updating if
1568 // 3D cells with faces with more than 4 sides are implemented
1569 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1570 int dofstart = 0;
1571 for (auto& edofs0 : _edofs[0])
1572 dofstart += edofs0.size();
1573
1574 // Transform DOFs on edges
1575 {
1576 auto& [v_size_t, matrix] = _etrans_inv.at(cell::type::interval)[0];
1577 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1578 {
1579 // Reverse an edge
1580 if (cell_info >> (face_start + e) & 1)
1581 {
1583 std::span(v_size_t),
1584 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1585 block_size);
1586 }
1587 dofstart += _edofs[1][e].size();
1588 }
1589 }
1590
1591 if (_cell_tdim == 3)
1592 {
1593 // Permute DOFs on faces
1594 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1595 {
1596 auto& trans = _etrans_inv.at(_cell_subentity_types[2][f]);
1597
1598 // Rotate a face
1599 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1600 {
1601 auto& m = trans[0];
1602 auto& v_size_t = std::get<0>(m);
1603 auto& matrix = std::get<1>(m);
1605 std::span(v_size_t),
1606 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1607 block_size);
1608 }
1609
1610 // Reflect a face
1611 if (cell_info >> (3 * f) & 1)
1612 {
1613 auto& m = trans[1];
1614 auto& v_size_t = std::get<0>(m);
1615 auto& matrix = std::get<1>(m);
1617 std::span(v_size_t),
1618 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1619 block_size);
1620 }
1621
1622 dofstart += _edofs[2][f].size();
1623 }
1624 }
1625 }
1626}
1627//-----------------------------------------------------------------------------
1628template <typename T>
1630 const std::span<T>& data, int block_size, std::uint32_t cell_info) const
1631{
1632 if (_dof_transformations_are_identity)
1633 return;
1634
1635 if (_cell_tdim >= 2)
1636 {
1637 // This assumes 3 bits are used per face. This will need updating if
1638 // 3D cells with faces with more than 4 sides are implemented
1639 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1640 int dofstart = 0;
1641 for (auto& edofs0 : _edofs[0])
1642 dofstart += edofs0.size();
1643
1644 // Transform DOFs on edges
1645 {
1646 auto& [v_size_t, matrix] = _etrans.at(cell::type::interval)[0];
1647 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1648 {
1649 // Reverse an edge
1650 if (cell_info >> (face_start + e) & 1)
1651 {
1653 std::span(v_size_t),
1654 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1655 block_size);
1656 }
1657
1658 dofstart += _edofs[1][e].size();
1659 }
1660 }
1661
1662 if (_cell_tdim == 3)
1663 {
1664 // Permute DOFs on faces
1665 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1666 {
1667 auto& trans = _etrans.at(_cell_subentity_types[2][f]);
1668
1669 // Reflect a face
1670 if (cell_info >> (3 * f) & 1)
1671 {
1672 auto& m = trans[1];
1673 auto& v_size_t = std::get<0>(m);
1674 auto& matrix = std::get<1>(m);
1676 std::span(v_size_t),
1677 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1678 block_size);
1679 }
1680
1681 // Rotate a face
1682 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1683 {
1684 auto& m = trans[0];
1685 auto& v_size_t = std::get<0>(m);
1686 auto& matrix = std::get<1>(m);
1688 std::span(v_size_t),
1689 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1690 block_size);
1691 }
1692 dofstart += _edofs[2][f].size();
1693 }
1694 }
1695 }
1696}
1697//-----------------------------------------------------------------------------
1698template <typename T>
1700 const std::span<T>& data, int block_size, std::uint32_t cell_info) const
1701{
1702 if (_dof_transformations_are_identity)
1703 return;
1704
1705 if (_cell_tdim >= 2)
1706 {
1707 // This assumes 3 bits are used per face. This will need updating if
1708 // 3D cells with faces with more than 4 sides are implemented
1709 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1710 int dofstart = 0;
1711 for (auto& edofs0 : _edofs[0])
1712 dofstart += edofs0.size();
1713
1714 // Transform DOFs on edges
1715 {
1716 auto& [v_size_t, matrix] = _etrans_invT.at(cell::type::interval)[0];
1717 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1718 {
1719 // Reverse an edge
1720 if (cell_info >> (face_start + e) & 1)
1721 {
1723 std::span(v_size_t),
1724 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1725 block_size);
1726 }
1727 dofstart += _edofs[1][e].size();
1728 }
1729 }
1730
1731 if (_cell_tdim == 3)
1732 {
1733 // Permute DOFs on faces
1734 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1735 {
1736 auto& trans = _etrans_invT.at(_cell_subentity_types[2][f]);
1737
1738 // Reflect a face
1739 if (cell_info >> (3 * f) & 1)
1740 {
1741 auto& m = trans[1];
1742 auto& v_size_t = std::get<0>(m);
1743 auto& matrix = std::get<1>(m);
1745 std::span(v_size_t),
1746 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1747 block_size);
1748 }
1749
1750 // Rotate a face
1751 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1752 {
1753 auto& m = trans[0];
1754 auto& v_size_t = std::get<0>(m);
1755 auto& matrix = std::get<1>(m);
1757 std::span(v_size_t),
1758 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1759 block_size);
1760 }
1761 dofstart += _edofs[2][f].size();
1762 }
1763 }
1764 }
1765}
1766//-----------------------------------------------------------------------------
1767template <typename T>
1769 const std::span<T>& data, int block_size, std::uint32_t cell_info) const
1770{
1771 if (_dof_transformations_are_identity)
1772 return;
1773
1774 if (_cell_tdim >= 2)
1775 {
1776 // This assumes 3 bits are used per face. This will need updating if
1777 // 3D cells with faces with more than 4 sides are implemented
1778 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1779 int dofstart = 0;
1780 for (auto& edofs0 : _edofs[0])
1781 dofstart += edofs0.size();
1782
1783 // Transform DOFs on edges
1784 {
1785 auto& [v_size_t, matrix] = _etransT.at(cell::type::interval)[0];
1786 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1787 {
1788 // Reverse an edge
1789 if (cell_info >> (face_start + e) & 1)
1790 {
1792 std::span(v_size_t),
1793 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1794 block_size);
1795 }
1796 dofstart += _edofs[1][e].size();
1797 }
1798 }
1799
1800 if (_cell_tdim == 3)
1801 {
1802 // Permute DOFs on faces
1803 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1804 {
1805 auto& trans = _etransT.at(_cell_subentity_types[2][f]);
1806
1807 // Rotate a face
1808 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1809 {
1810 auto& m = trans[0];
1811 auto& v_size_t = std::get<0>(m);
1812 auto& matrix = std::get<1>(m);
1814 std::span(v_size_t),
1815 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1816 block_size);
1817 }
1818
1819 // Reflect a face
1820 if (cell_info >> (3 * f) & 1)
1821 {
1822 auto& m = trans[1];
1823 auto& v_size_t = std::get<0>(m);
1824 auto& matrix = std::get<1>(m);
1826 std::span(v_size_t),
1827 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1828 block_size);
1829 }
1830 dofstart += _edofs[2][f].size();
1831 }
1832 }
1833 }
1834}
1835//-----------------------------------------------------------------------------
1836template <typename T>
1838 const std::span<T>& data, int block_size, std::uint32_t cell_info) const
1839{
1840 if (_dof_transformations_are_identity)
1841 return;
1842
1843 if (_cell_tdim >= 2)
1844 {
1845 // This assumes 3 bits are used per face. This will need updating if
1846 // 3D cells with faces with more than 4 sides are implemented
1847 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1848 int dofstart = 0;
1849 for (auto& edofs0 : _edofs[0])
1850 dofstart += edofs0.size();
1851
1852 // Transform DOFs on edges
1853 {
1854 auto& [v_size_t, matrix] = _etrans_inv.at(cell::type::interval)[0];
1855 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1856 {
1857 // Reverse an edge
1858 if (cell_info >> (face_start + e) & 1)
1859 {
1861 std::span(v_size_t),
1862 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1863 block_size);
1864 }
1865 dofstart += _edofs[1][e].size();
1866 }
1867 }
1868
1869 if (_cell_tdim == 3)
1870 {
1871 // Permute DOFs on faces
1872 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1873 {
1874 auto& trans = _etrans_inv.at(_cell_subentity_types[2][f]);
1875
1876 // Rotate a face
1877 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1878 {
1879 auto& m = trans[0];
1880 auto& v_size_t = std::get<0>(m);
1881 auto& matrix = std::get<1>(m);
1883 std::span(v_size_t),
1884 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1885 block_size);
1886 }
1887
1888 // Reflect a face
1889 if (cell_info >> (3 * f) & 1)
1890 {
1891 auto& m = trans[1];
1892 auto& v_size_t = std::get<0>(m);
1893 auto& matrix = std::get<1>(m);
1895 std::span(v_size_t),
1896 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1897 block_size);
1898 }
1899
1900 dofstart += _edofs[2][f].size();
1901 }
1902 }
1903 }
1904}
1905//-----------------------------------------------------------------------------
1906
1907} // namespace basix
A finite element.
Definition: finite-element.h:144
std::map< cell::type, std::pair< std::vector< double >, std::array< std::size_t, 3 > > > entity_transformations() const
Definition: finite-element.cpp:1418
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:1042
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:1837
cell::type cell_type() const
Definition: finite-element.cpp:1116
int highest_complete_degree() const
Definition: finite-element.cpp:1122
FiniteElement & operator=(const FiniteElement &element)=default
Assignment operator.
element::lagrange_variant lagrange_variant() const
Definition: finite-element.cpp:1460
FiniteElement(FiniteElement &&element)=default
Move constructor.
bool operator==(const FiniteElement &e) const
Definition: finite-element.cpp:994
const std::array< std::vector< std::pair< std::vector< double >, std::array< std::size_t, 2 > > >, 4 > & x() const
Definition: finite-element.cpp:1437
void apply_dof_transformation(const std::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1352
bool dof_transformations_are_permutations() const
Definition: finite-element.cpp:1140
maps::type map_type() const
Definition: finite-element.cpp:1136
std::vector< std::tuple< std::vector< FiniteElement >, std::vector< int > > > get_tensor_product_representation() const
Definition: finite-element.cpp:1478
const std::vector< std::vector< std::vector< int > > > & entity_dofs() const
Definition: finite-element.cpp:1157
const std::pair< std::vector< double >, std::array< std::size_t, 2 > > & dual_matrix() const
Definition: finite-element.cpp:1424
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:1699
int interpolation_nderivs() const
The number of derivatives needed when interpolating.
Definition: finite-element.cpp:1472
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, 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:603
FiniteElement & operator=(FiniteElement &&element)=default
Move assignment operator.
const std::vector< std::vector< std::vector< int > > > & entity_closure_dofs() const
Definition: finite-element.cpp:1163
const std::pair< std::vector< double >, std::array< std::size_t, 2 > > & points() const
Definition: finite-element.cpp:1237
int degree() const
Definition: finite-element.cpp:1118
void apply_transpose_dof_transformation(const std::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1423
const std::vector< std::size_t > & value_shape() const
Definition: finite-element.cpp:1127
const std::pair< std::vector< double >, std::array< std::size_t, 2 > > & coefficient_matrix() const
Definition: finite-element.cpp:1450
bool has_tensor_product_factorisation() const
Definition: finite-element.cpp:1455
~FiniteElement()=default
Destructor.
const std::pair< std::vector< double >, std::array< std::size_t, 2 > > & wcoeffs() const
Definition: finite-element.cpp:1430
int dim() const
Definition: finite-element.cpp:1132
int highest_degree() const
Definition: finite-element.cpp:1120
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:1243
std::array< std::size_t, 4 > tabulate_shape(std::size_t nd, std::size_t num_points) const
Definition: finite-element.cpp:1028
element::dpc_variant dpc_variant() const
Definition: finite-element.cpp:1465
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:1768
void apply_dof_transformation_to_transpose(const std::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1629
element::family family() const
Definition: finite-element.cpp:1134
void permute_dofs(const std::span< std::int32_t > &dofs, std::uint32_t cell_info) const
Definition: finite-element.cpp:1310
bool dof_transformations_are_identity() const
Definition: finite-element.cpp:1145
void apply_inverse_dof_transformation(const std::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1559
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:1278
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:1151
bool interpolation_is_identity() const
Definition: finite-element.cpp:1467
std::pair< std::vector< double >, std::array< std::size_t, 3 > > base_transformations() const
Get the base transformations.
Definition: finite-element.cpp:1169
const std::array< std::vector< std::pair< std::vector< double >, std::array< std::size_t, 4 > > >, 4 > & M() const
Definition: finite-element.cpp:1444
std::function< void(O &, const P &, const Q &, double, const R &)> map_fn() const
Definition: finite-element.h:630
bool discontinuous() const
Definition: finite-element.cpp:1138
void apply_inverse_transpose_dof_transformation(const std::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1492
void unpermute_dofs(const std::span< std::int32_t > &dofs, std::uint32_t cell_info) const
Definition: finite-element.cpp:1364
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:111
impl::cmdspan2_t cmdspan2_t
Typedef for mdspan.
Definition: finite-element.h:109
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
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:1485
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, bool discontinuous, int highest_complete_degree, int highest_degree)
Definition: finite-element.cpp:464