Basix
Loading...
Searching...
No Matches
finite-element.h
1// Copyright (c) 2020-2024 Chris Richardson, Matthew Scroggs and Garth . Wells
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
512 const std::vector<std::size_t>& value_shape() const { return _value_shape; }
513
518 int dim() const { return _coeffs.second[0]; }
519
522 element::family family() const { return _family; }
523
527 {
528 return _lagrange_variant;
529 }
530
533 element::dpc_variant dpc_variant() const { return _dpc_variant; }
534
537 maps::type map_type() const { return _map_type; }
538
541 sobolev::space sobolev_space() const { return _sobolev_space; }
542
547 bool discontinuous() const { return _discontinuous; }
548
552 {
553 return _dof_transformations_are_permutations;
554 }
555
558 {
559 return _dof_transformations_are_identity;
560 }
561
576 std::pair<std::vector<F>, std::array<std::size_t, 3>>
577 push_forward(impl::mdspan_t<const F, 3> U, impl::mdspan_t<const F, 3> J,
578 std::span<const F> detJ, impl::mdspan_t<const F, 3> K) const;
579
587 std::pair<std::vector<F>, std::array<std::size_t, 3>>
588 pull_back(impl::mdspan_t<const F, 3> u, impl::mdspan_t<const F, 3> J,
589 std::span<const F> detJ, impl::mdspan_t<const F, 3> K) const;
590
622 template <typename O, typename P, typename Q, typename R>
623 std::function<void(O&, const P&, const Q&, F, const R&)> map_fn() const
624 {
625 switch (_map_type)
626 {
627 case maps::type::identity:
628 return [](O& u, const P& U, const Q&, F, const R&)
629 {
630 assert(U.extent(0) == u.extent(0));
631 assert(U.extent(1) == u.extent(1));
632 for (std::size_t i = 0; i < U.extent(0); ++i)
633 for (std::size_t j = 0; j < U.extent(1); ++j)
634 u(i, j) = U(i, j);
635 };
636 case maps::type::covariantPiola:
637 return [](O& u, const P& U, const Q& J, F detJ, const R& K)
638 { maps::covariant_piola(u, U, J, detJ, K); };
639 case maps::type::contravariantPiola:
640 return [](O& u, const P& U, const Q& J, F detJ, const R& K)
642 case maps::type::doubleCovariantPiola:
643 return [](O& u, const P& U, const Q& J, F detJ, const R& K)
645 case maps::type::doubleContravariantPiola:
646 return [](O& u, const P& U, const Q& J, F detJ, const R& K)
648 default:
649 throw std::runtime_error("Map not implemented");
650 }
651 }
652
660 const std::vector<std::vector<std::vector<int>>>& entity_dofs() const
661 {
662 return _edofs;
663 }
664
674 const std::vector<std::vector<std::vector<int>>>& entity_closure_dofs() const
675 {
676 return _e_closure_dofs;
677 }
678
760 std::pair<std::vector<F>, std::array<std::size_t, 3>>
761 base_transformations() const;
762
767 std::map<cell::type, std::pair<std::vector<F>, std::array<std::size_t, 3>>>
769 {
770 return _entity_transformations;
771 }
772
795 void permute(std::span<std::int32_t> d, 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
803 if (_dof_transformations_are_identity)
804 return;
805 else
807 }
808
826 void permute_inv(std::span<std::int32_t> d, std::uint32_t cell_info) const
827 {
828 if (!_dof_transformations_are_permutations)
829 {
830 throw std::runtime_error(
831 "The DOF transformations for this element are not permutations");
832 }
833
834 if (_dof_transformations_are_identity)
835 return;
836 else
838 }
839
855 void permute_subentity_closure(std::span<std::int32_t> d,
856 std::uint32_t cell_info,
858 {
860
861 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
862
863 std::uint32_t entity_info = 0;
864 switch (entity_dim)
865 {
866 case 0:
867 return;
868 case 1:
870 break;
871 case 2:
872 entity_info = cell_info >> (3 * entity_index) & 7;
873 break;
874 default:
875 throw std::runtime_error("Unsupported cell dimension");
876 }
878 }
879
892 void permute_subentity_closure_inv(std::span<std::int32_t> d,
893 std::uint32_t cell_info,
895 int entity_index) const
896 {
898
899 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
900
901 std::uint32_t entity_info;
902 switch (entity_dim)
903 {
904 case 0:
905 return;
906 case 1:
908 break;
909 case 2:
910 entity_info = cell_info >> (3 * entity_index) & 7;
911 break;
912 default:
913 throw std::runtime_error("Unsupported cell dimension");
914 }
916 }
917
932
933 void permute_subentity_closure(std::span<std::int32_t> d,
934 std::uint32_t entity_info,
936 {
937 if (!_dof_transformations_are_permutations)
938 {
939 throw std::runtime_error(
940 "The DOF transformations for this element are not permutations");
941 }
942
944
945 if (entity_dim == 0)
946 return;
947
948 auto& perm = _subentity_closure_perm.at(entity_type);
949 if (entity_dim == 1)
950 {
951 if (entity_info & 1)
952 {
954 }
955 }
956 else if (entity_dim == 2)
957 {
958 // Rotate a face
959 for (std::uint32_t r = 0; r < (entity_info >> 1 & 3); ++r)
960 {
962 }
963
964 // Reflect a face (post rotate)
965 if (entity_info & 1)
966 {
968 }
969 }
970 else
971 {
972 throw std::runtime_error(
973 "Invalid dimension for permute_subentity_closure");
974 }
975 }
976
988 void permute_subentity_closure_inv(std::span<std::int32_t> d,
989 std::uint32_t entity_info,
991 {
992 if (!_dof_transformations_are_permutations)
993 {
994 throw std::runtime_error(
995 "The DOF transformations for this element are not permutations");
996 }
997
999
1000 if (entity_dim == 0)
1001 return;
1002
1003 auto& perm = _subentity_closure_perm_inv.at(entity_type);
1004 if (entity_dim == 1)
1005 {
1006 if (entity_info & 1)
1007 {
1009 }
1010 }
1011 else if (entity_dim == 2)
1012 {
1013 // Reflect a face (pre rotate)
1014 if (entity_info & 1)
1015 {
1017 }
1018
1019 // Rotate a face
1020 for (std::uint32_t r = 0; r < (entity_info >> 1 & 3); ++r)
1021 {
1023 }
1024 }
1025 else
1026 {
1027 throw std::runtime_error(
1028 "Invalid dimension for permute_subentity_closure");
1029 }
1030 }
1031
1062 template <typename T>
1063 void T_apply(std::span<T> u, int n, std::uint32_t cell_info) const;
1064
1075 template <typename T>
1076 void Tt_apply(std::span<T> u, int n, std::uint32_t cell_info) const;
1077
1089 template <typename T>
1090 void Tt_inv_apply(std::span<T> u, int n, std::uint32_t cell_info) const;
1091
1102 template <typename T>
1103 void Tinv_apply(std::span<T> u, int n, std::uint32_t cell_info) const;
1104
1115 template <typename T>
1116 void Tt_apply_right(std::span<T> u, int n, std::uint32_t cell_info) const;
1117
1127 template <typename T>
1128 void T_apply_right(std::span<T> u, int n, std::uint32_t cell_info) const;
1129
1140 template <typename T>
1141 void Tinv_apply_right(std::span<T> u, int n, std::uint32_t cell_info) const;
1142
1153 template <typename T>
1154 void Tt_inv_apply_right(std::span<T> u, int n, std::uint32_t cell_info) const;
1155
1162 const std::pair<std::vector<F>, std::array<std::size_t, 2>>& points() const
1163 {
1164 return _points;
1165 }
1166
1218 const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
1220 {
1221 return _matM;
1222 }
1223
1229 const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
1231 {
1232 return _dual_matrix;
1233 }
1234
1271 const std::pair<std::vector<F>, std::array<std::size_t, 2>>& wcoeffs() const
1272 {
1273 return _wcoeffs;
1274 }
1275
1280 const std::array<
1281 std::vector<std::pair<std::vector<F>, std::array<std::size_t, 2>>>, 4>&
1282 x() const
1283 {
1284 return _x;
1285 }
1286
1323 const std::array<
1324 std::vector<std::pair<std::vector<F>, std::array<std::size_t, 4>>>, 4>&
1325 M() const
1326 {
1327 return _M;
1328 }
1329
1333 const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
1335 {
1336 return _coeffs;
1337 }
1338
1350 {
1351 return !_tensor_factors.empty();
1352 }
1353
1366 std::vector<std::vector<FiniteElement<F>>>
1368 {
1370 throw std::runtime_error("Element has no tensor product representation.");
1371 return _tensor_factors;
1372 }
1373
1378 bool interpolation_is_identity() const { return _interpolation_is_identity; }
1379
1381 int interpolation_nderivs() const { return _interpolation_nderivs; }
1382
1384 const std::vector<int>& dof_ordering() const { return _dof_ordering; }
1385
1386private:
1393 template <typename T, bool post>
1394 void permute_data(
1395 std::span<T> data, int block_size, std::uint32_t cell_info,
1396 const std::map<cell::type, std::vector<std::vector<std::size_t>>>& eperm)
1397 const;
1398
1399 using array2_t = std::pair<std::vector<F>, std::array<std::size_t, 2>>;
1400 using array3_t = std::pair<std::vector<F>, std::array<std::size_t, 3>>;
1401 using trans_data_t
1402 = std::vector<std::pair<std::vector<std::size_t>, array2_t>>;
1403
1410 template <typename T, bool post, typename OP>
1411 void
1412 transform_data(std::span<T> data, int block_size, std::uint32_t cell_info,
1413 const std::map<cell::type, trans_data_t>& etrans, OP op) const;
1414
1415 // Cell type
1416 cell::type _cell_type;
1417
1418 // Polyset type
1419 polyset::type _poly_type;
1420
1421 // Topological dimension of the cell
1422 std::size_t _cell_tdim;
1423
1424 // Topological dimension of the cell
1425 std::vector<std::vector<cell::type>> _cell_subentity_types;
1426
1427 // Finite element family
1428 element::family _family;
1429
1430 // Lagrange variant
1431 element::lagrange_variant _lagrange_variant;
1432
1433 // DPC variant
1434 element::dpc_variant _dpc_variant;
1435
1436 // Degree that was input when creating the element
1437 int _degree;
1438
1439 // Degree
1440 int _interpolation_nderivs;
1441
1442 // Highest degree polynomial in element's polyset
1443 int _embedded_superdegree;
1444
1445 // Highest degree space that is a subspace of element's polyset
1446 int _embedded_subdegree;
1447
1448 // Value shape
1449 std::vector<std::size_t> _value_shape;
1450
1452 maps::type _map_type;
1453
1455 sobolev::space _sobolev_space;
1456
1457 // Shape function coefficient of expansion sets on cell. If shape
1458 // function is given by @f$\psi_i = \sum_{k} \phi_{k}
1459 // \alpha^{i}_{k}@f$, then _coeffs(i, j) = @f$\alpha^i_k@f$. ie
1460 // _coeffs.row(i) are the expansion coefficients for shape function i
1461 // (@f$\psi_{i}@f$).
1462 std::pair<std::vector<F>, std::array<std::size_t, 2>> _coeffs;
1463
1464 // Dofs associated with each cell (sub-)entity
1465 std::vector<std::vector<std::vector<int>>> _edofs;
1466
1467 // Dofs associated with the closdure of each cell (sub-)entity
1468 std::vector<std::vector<std::vector<int>>> _e_closure_dofs;
1469
1470 // Entity transformations
1471 std::map<cell::type, array3_t> _entity_transformations;
1472
1473 // Set of points used for point evaluation
1474 // Experimental - currently used for an implementation of
1475 // "tabulate_dof_coordinates" Most useful for Lagrange. This may change or go
1476 // away. For non-Lagrange elements, these points will be used in combination
1477 // with _interpolation_matrix to perform interpolation
1478 std::pair<std::vector<F>, std::array<std::size_t, 2>> _points;
1479
1480 // Interpolation points on the cell. The shape is (entity_dim, num
1481 // entities of given dimension, num_points, tdim)
1482 std::array<std::vector<std::pair<std::vector<F>, std::array<std::size_t, 2>>>,
1483 4>
1484 _x;
1485
1487 std::pair<std::vector<F>, std::array<std::size_t, 2>> _matM;
1488
1489 // Indicates whether or not the DOF transformations are all
1490 // permutations
1491 bool _dof_transformations_are_permutations;
1492
1493 // Indicates whether or not the DOF transformations are all identity
1494 bool _dof_transformations_are_identity;
1495
1496 // The entity permutations (factorised). This will only be set if
1497 // _dof_transformations_are_permutations is True and
1498 // _dof_transformations_are_identity is False
1499 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm;
1500
1501 // The reverse entity permutations (factorised). This will only be set
1502 // if _dof_transformations_are_permutations is True and
1503 // _dof_transformations_are_identity is False
1504 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_inv;
1505
1506 // The entity transformations in precomputed form
1507 std::map<cell::type, trans_data_t> _etrans;
1508
1509 // The transposed entity transformations in precomputed form
1510 std::map<cell::type, trans_data_t> _etransT;
1511
1512 // The inverse entity transformations in precomputed form
1513 std::map<cell::type, trans_data_t> _etrans_inv;
1514
1515 // The inverse transpose entity transformations in precomputed form
1516 std::map<cell::type, trans_data_t> _etrans_invT;
1517
1518 // The subentity closure permutations (factorised). This will only be set if
1519 // _dof_transformations_are_permutations is True
1520 std::map<cell::type, std::vector<std::vector<std::size_t>>>
1521 _subentity_closure_perm;
1522
1523 // The inverse subentity closure permutations (factorised). This will only be
1524 // set if _dof_transformations_are_permutations is True
1525 std::map<cell::type, std::vector<std::vector<std::size_t>>>
1526 _subentity_closure_perm_inv;
1527
1528 // Indicates whether or not this is the discontinuous version of the
1529 // element
1530 bool _discontinuous;
1531
1532 // The dual matrix
1533 std::pair<std::vector<F>, std::array<std::size_t, 2>> _dual_matrix;
1534
1535 // Dof reordering for different element dof layout compatibility.
1536 // The reference basix layout is ordered by entity, i.e. dofs on
1537 // vertices, followed by edges, faces, then internal dofs.
1538 // _dof_ordering stores the map to the new order required, e.g.
1539 // for a P2 triangle, _dof_ordering=[0 3 5 1 2 4] will place
1540 // dofs 0, 3, 5 on the vertices and 1, 2, 4, on the edges.
1541 std::vector<int> _dof_ordering;
1542
1543 // Tensor product representation
1544 // Entries of tuple are (list of elements on an interval, permutation
1545 // of DOF numbers)
1546 // @todo: For vector-valued elements, a tensor product type and a
1547 // scaling factor may additionally be needed.
1548 std::vector<std::vector<FiniteElement>> _tensor_factors;
1549
1550 // Is the interpolation matrix an identity?
1551 bool _interpolation_is_identity;
1552
1553 // The coefficients that define the polynomial set in terms of the
1554 // orthonormal polynomials
1555 std::pair<std::vector<F>, std::array<std::size_t, 2>> _wcoeffs;
1556
1557 // Interpolation matrices for each entity
1558 using array4_t
1559 = std::vector<std::pair<std::vector<F>, std::array<std::size_t, 4>>>;
1560 std::array<array4_t, 4> _M;
1561};
1562
1588template <std::floating_point T>
1589FiniteElement<T> create_custom_element(
1590 cell::type cell_type, const std::vector<std::size_t>& value_shape,
1591 impl::mdspan_t<const T, 2> wcoeffs,
1592 const std::array<std::vector<impl::mdspan_t<const T, 2>>, 4>& x,
1593 const std::array<std::vector<impl::mdspan_t<const T, 4>>, 4>& M,
1594 int interpolation_nderivs, maps::type map_type,
1595 sobolev::space sobolev_space, bool discontinuous, int embedded_subdegree,
1596 int embedded_superdegree, polyset::type poly_type);
1597
1609template <std::floating_point T>
1610FiniteElement<T> create_element(element::family family, cell::type cell,
1611 int degree, element::lagrange_variant lvariant,
1612 element::dpc_variant dvariant,
1613 bool discontinuous,
1614 std::vector<int> dof_ordering = {});
1615
1627std::vector<int> tp_dof_ordering(element::family family, cell::type cell,
1628 int degree, element::lagrange_variant lvariant,
1629 element::dpc_variant dvariant,
1630 bool discontinuous);
1631
1644template <std::floating_point T>
1645std::vector<std::vector<FiniteElement<T>>>
1646tp_factors(element::family family, cell::type cell, int degree,
1648 bool discontinuous, std::vector<int> dof_ordering);
1649
1661template <std::floating_point T>
1662FiniteElement<T>
1663create_tp_element(element::family family, cell::type cell, int degree,
1665 element::dpc_variant dvariant, bool discontinuous);
1666
1669std::string version();
1670
1671//-----------------------------------------------------------------------------
1672template <std::floating_point F>
1673template <typename T, bool post>
1674void FiniteElement<F>::permute_data(
1675 std::span<T> data, int block_size, std::uint32_t cell_info,
1676 const std::map<cell::type, std::vector<std::vector<std::size_t>>>& eperm)
1677 const
1678{
1679 if (_cell_tdim >= 2)
1680 {
1681 // This assumes 3 bits are used per face. This will need updating if 3D
1682 // cells with faces with more than 4 sides are implemented
1683 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1684
1685 // Permute DOFs on edges
1686 {
1687 auto& trans = eperm.at(cell::type::interval)[0];
1688 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1689 {
1690 // Reverse an edge
1691 if (cell_info >> (face_start + e) & 1)
1692 {
1693 precompute::apply_permutation_mapped(trans, data, _edofs[1][e],
1694 block_size);
1695 }
1696 }
1697 }
1698
1699 if (_cell_tdim == 3)
1700 {
1701 // Permute DOFs on faces
1702 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1703 {
1704 auto& trans = eperm.at(_cell_subentity_types[2][f]);
1705
1706 // Reflect a face (pre rotate)
1707 if (!post and cell_info >> (3 * f) & 1)
1708 {
1709 precompute::apply_permutation_mapped(trans[1], data, _edofs[2][f],
1710 block_size);
1711 }
1712
1713 // Rotate a face
1714 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1715 {
1716 precompute::apply_permutation_mapped(trans[0], data, _edofs[2][f],
1717 block_size);
1718 }
1719
1720 // Reflect a face (post rotate)
1721 if (post and cell_info >> (3 * f) & 1)
1722 {
1723 precompute::apply_permutation_mapped(trans[1], data, _edofs[2][f],
1724 block_size);
1725 }
1726 }
1727 }
1728 }
1729}
1730//-----------------------------------------------------------------------------
1731template <std::floating_point F>
1732template <typename T, bool post, typename OP>
1733void FiniteElement<F>::transform_data(
1734 std::span<T> data, int block_size, std::uint32_t cell_info,
1735 const std::map<cell::type, trans_data_t>& etrans, OP op) const
1736{
1737 if (_cell_tdim >= 2)
1738 {
1739 // This assumes 3 bits are used per face. This will need updating if
1740 // 3D cells with faces with more than 4 sides are implemented
1741 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1742 int dofstart = 0;
1743 for (auto& edofs0 : _edofs[0])
1744 dofstart += edofs0.size();
1745
1746 // Transform DOFs on edges
1747 {
1748 auto& [v_size_t, matrix] = etrans.at(cell::type::interval)[0];
1749 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1750 {
1751 // Reverse an edge
1752 if (cell_info >> (face_start + e) & 1)
1753 {
1754 op(std::span(v_size_t),
1755 mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
1756 dofstart, block_size);
1757 }
1758 dofstart += _edofs[1][e].size();
1759 }
1760 }
1761
1762 if (_cell_tdim == 3)
1763 {
1764 // Permute DOFs on faces
1765 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1766 {
1767 auto& trans = etrans.at(_cell_subentity_types[2][f]);
1768
1769 // Reflect a face (pre rotation)
1770 if (!post and cell_info >> (3 * f) & 1)
1771 {
1772 const auto& m = trans[1];
1773 const auto& v_size_t = std::get<0>(m);
1774 const auto& matrix = std::get<1>(m);
1775 op(std::span(v_size_t),
1776 mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
1777 dofstart, block_size);
1778 }
1779
1780 // Rotate a face
1781 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1782 {
1783 const auto& m = trans[0];
1784 const auto& v_size_t = std::get<0>(m);
1785 const auto& matrix = std::get<1>(m);
1786 op(std::span(v_size_t),
1787 mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
1788 dofstart, block_size);
1789 }
1790
1791 // Reflect a face (post rotation)
1792 if (post and cell_info >> (3 * f) & 1)
1793 {
1794 const auto& m = trans[1];
1795 const auto& v_size_t = std::get<0>(m);
1796 const auto& matrix = std::get<1>(m);
1797 op(std::span(v_size_t),
1798 mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
1799 dofstart, block_size);
1800 }
1801
1802 dofstart += _edofs[2][f].size();
1803 }
1804 }
1805 }
1806}
1807//-----------------------------------------------------------------------------
1808template <std::floating_point F>
1809template <typename T>
1810void FiniteElement<F>::T_apply(std::span<T> u, int n,
1811 std::uint32_t cell_info) const
1812{
1813 if (_dof_transformations_are_identity)
1814 return;
1815
1816 if (_dof_transformations_are_permutations)
1818 else
1819 {
1821 precompute::apply_matrix<F, T>);
1822 }
1823}
1824//-----------------------------------------------------------------------------
1825template <std::floating_point F>
1826template <typename T>
1827void FiniteElement<F>::Tt_apply(std::span<T> u, int n,
1828 std::uint32_t cell_info) const
1829{
1830 if (_dof_transformations_are_identity)
1831 return;
1832 else if (_dof_transformations_are_permutations)
1833 permute_data<T, true>(u, n, cell_info, _eperm_inv);
1834 else
1835 {
1837 precompute::apply_matrix<F, T>);
1838 }
1839}
1840//-----------------------------------------------------------------------------
1841template <std::floating_point F>
1842template <typename T>
1843void FiniteElement<F>::Tt_inv_apply(std::span<T> u, int n,
1844 std::uint32_t cell_info) const
1845{
1846 if (_dof_transformations_are_identity)
1847 return;
1848 else if (_dof_transformations_are_permutations)
1850 else
1851 {
1852 transform_data<T, false>(u, n, cell_info, _etrans_invT,
1853 precompute::apply_matrix<F, T>);
1854 }
1855}
1856//-----------------------------------------------------------------------------
1857template <std::floating_point F>
1858template <typename T>
1859void FiniteElement<F>::Tinv_apply(std::span<T> u, int n,
1860 std::uint32_t cell_info) const
1861{
1862 if (_dof_transformations_are_identity)
1863 return;
1864 else if (_dof_transformations_are_permutations)
1865 permute_data<T, true>(u, n, cell_info, _eperm_inv);
1866 else
1867 {
1868 transform_data<T, true>(u, n, cell_info, _etrans_inv,
1869 precompute::apply_matrix<F, T>);
1870 }
1871}
1872//-----------------------------------------------------------------------------
1873template <std::floating_point F>
1874template <typename T>
1875void FiniteElement<F>::Tt_apply_right(std::span<T> u, int n,
1876 std::uint32_t cell_info) const
1877{
1878 if (_dof_transformations_are_identity)
1879 return;
1880 else if (_dof_transformations_are_permutations)
1881 {
1882 assert(u.size() % n == 0);
1883 const int step = u.size() / n;
1884 for (int i = 0; i < n; ++i)
1885 {
1886 std::span<T> dblock(u.data() + i * step, step);
1888 }
1889 }
1890 else
1891 {
1893 precompute::apply_tranpose_matrix_right<F, T>);
1894 }
1895}
1896//-----------------------------------------------------------------------------
1897template <std::floating_point F>
1898template <typename T>
1900 std::uint32_t cell_info) const
1901{
1902 if (_dof_transformations_are_identity)
1903 return;
1904 else if (_dof_transformations_are_permutations)
1905 {
1906 assert(u.size() % n == 0);
1907 const int step = u.size() / n;
1908 for (int i = 0; i < n; ++i)
1909 {
1910 std::span<T> dblock(u.data() + i * step, step);
1912 }
1913 }
1914 else
1915 {
1916 transform_data<T, false>(u, n, cell_info, _etrans_invT,
1917 precompute::apply_tranpose_matrix_right<F, T>);
1918 }
1919}
1920//-----------------------------------------------------------------------------
1921template <std::floating_point F>
1922template <typename T>
1923void FiniteElement<F>::T_apply_right(std::span<T> u, int n,
1924 std::uint32_t cell_info) const
1925{
1926 if (_dof_transformations_are_identity)
1927 return;
1928 else if (_dof_transformations_are_permutations)
1929 {
1930 assert(u.size() % n == 0);
1931 const int step = u.size() / n;
1932 for (int i = 0; i < n; ++i)
1933 {
1934 std::span<T> dblock(u.data() + i * step, step);
1935 permute_data<T, true>(dblock, 1, cell_info, _eperm_inv);
1936 }
1937 }
1938 else
1939 {
1941 precompute::apply_tranpose_matrix_right<F, T>);
1942 }
1943}
1944//-----------------------------------------------------------------------------
1945template <std::floating_point F>
1946template <typename T>
1948 std::uint32_t cell_info) const
1949{
1950 if (_dof_transformations_are_identity)
1951 return;
1952 else if (_dof_transformations_are_permutations)
1953 {
1954 assert(u.size() % n == 0);
1955 const int step = u.size() / n;
1956 for (int i = 0; i < n; ++i)
1957 {
1958 std::span<T> dblock(u.data() + i * step, step);
1959 permute_data<T, true>(dblock, 1, cell_info, _eperm_inv);
1960 }
1961 }
1962 else
1963 {
1964 transform_data<T, true>(u, n, cell_info, _etrans_inv,
1965 precompute::apply_tranpose_matrix_right<F, T>);
1966 }
1967}
1968//-----------------------------------------------------------------------------
1969
1970} // namespace basix
A finite element.
Definition finite-element.h:139
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.
const std::array< std::vector< std::pair< std::vector< F >, std::array< std::size_t, 4 > > >, 4 > & M() const
Get the interpolation matrices for each subentity.
Definition finite-element.h:1325
void Tinv_apply_right(std::span< T > u, int n, std::uint32_t cell_info) const
Right(post)-apply the inverse of the operator applied by T_apply().
Definition finite-element.h:1899
std::pair< std::vector< F >, std::array< std::size_t, 3 > > base_transformations() const
Get the base transformations.
Definition finite-element.cpp:1632
void T_apply(std::span< T > u, int n, std::uint32_t cell_info) const
Transform basis functions from the reference element ordering and orientation to the globally consist...
Definition finite-element.h:1810
void Tt_apply(std::span< T > u, int n, std::uint32_t cell_info) const
Apply the transpose of the operator applied by T_apply().
Definition finite-element.h:1827
void permute_subentity_closure_inv(std::span< std::int32_t > d, std::uint32_t cell_info, cell::type entity_type, int entity_index) const
Perform the inverse of the operation applied by permute_subentity_closure().
Definition finite-element.h:892
bool dof_transformations_are_identity() const
Indicates is the dof transformations are all the identity.
Definition finite-element.h:557
bool operator==(const FiniteElement &e) const
Check if two elements are the same.
Definition finite-element.cpp:1451
std::map< cell::type, std::pair< std::vector< F >, std::array< std::size_t, 3 > > > entity_transformations() const
Return the entity dof transformation matrices.
Definition finite-element.h:768
const std::vector< int > & dof_ordering() const
Get dof layout.
Definition finite-element.h:1384
void permute_subentity_closure_inv(std::span< std::int32_t > d, std::uint32_t entity_info, cell::type entity_type) const
Perform the inverse of the operation applied by permute_subentity_closure().
Definition finite-element.h:988
int embedded_subdegree() const
Highest degree n such that a Lagrange (or vector Lagrange) element of degree n is a subspace of this ...
Definition finite-element.h:505
const std::vector< std::vector< std::vector< int > > > & entity_dofs() const
Get the dofs on each topological entity: (vertices, edges, faces, cell) in that order.
Definition finite-element.h:660
void T_apply_right(std::span< T > u, int n, std::uint32_t cell_info) const
Right(post)-apply the operator applied by T_apply().
Definition finite-element.h:1923
FiniteElement(FiniteElement &&element)=default
Move constructor.
bool has_tensor_product_factorisation() const
Indicates whether or not this element can be represented as a product of elements defined on lower-di...
Definition finite-element.h:1349
int interpolation_nderivs() const
The number of derivatives needed when interpolating.
Definition finite-element.h:1381
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:1540
int embedded_superdegree() const
Lowest degree n such that the highest degree polynomial in this element is contained in a Lagrange (o...
Definition finite-element.h:500
int dim() const
Dimension of the finite element space.
Definition finite-element.h:518
void permute(std::span< std::int32_t > d, std::uint32_t cell_info) const
Permute indices associated with degree-of-freedoms on the reference element ordering to the globally ...
Definition finite-element.h:795
std::size_t hash() const
Get a unique hash of this element.
Definition finite-element.cpp:1490
void permute_inv(std::span< std::int32_t > d, std::uint32_t cell_info) const
Perform the inverse of the operation applied by permute().
Definition finite-element.h:826
void Tt_apply_right(std::span< T > u, int n, std::uint32_t cell_info) const
Right(post)-apply the transpose of the operator applied by T_apply().
Definition finite-element.h:1875
const std::vector< std::size_t > & value_shape() const
Element value tensor shape.
Definition finite-element.h:512
void Tt_inv_apply(std::span< T > u, int n, std::uint32_t cell_info) const
Apply the inverse transpose of the operator applied by T_apply().
Definition finite-element.h:1843
void permute_subentity_closure(std::span< std::int32_t > d, std::uint32_t cell_info, cell::type entity_type, int entity_index) const
Permute indices associated with degree-of-freedoms on the closure of a sub-entity of the reference el...
Definition finite-element.h:855
element::lagrange_variant lagrange_variant() const
Lagrange variant of the element.
Definition finite-element.h:526
void Tt_inv_apply_right(std::span< T > u, int n, std::uint32_t cell_info) const
Right(post)-apply the transpose inverse of the operator applied by T_apply().
Definition finite-element.h:1947
const std::array< std::vector< std::pair< std::vector< F >, std::array< std::size_t, 2 > > >, 4 > & x() const
Get the interpolation points for each subentity.
Definition finite-element.h:1282
int degree() const
Get the element polynomial degree.
Definition finite-element.h:494
void permute_subentity_closure(std::span< std::int32_t > d, std::uint32_t entity_info, cell::type entity_type) const
Permute indices associated with degree-of-freedoms on the closure of a sub-entity of the reference el...
Definition finite-element.h:933
void Tinv_apply(std::span< T > u, int n, std::uint32_t cell_info) const
Apply the inverse of the operator applied by T_apply().
Definition finite-element.h:1859
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & points() const
Return the interpolation points.
Definition finite-element.h:1162
sobolev::space sobolev_space() const
Underlying Sobolev space for this element.
Definition finite-element.h:541
const std::vector< std::vector< std::vector< int > > > & entity_closure_dofs() const
Get the dofs on the closure of each topological entity: (vertices, edges, faces, cell) in that order.
Definition finite-element.h:674
FiniteElement(const FiniteElement &element)=default
Copy constructor.
std::array< std::size_t, 4 > tabulate_shape(std::size_t nd, std::size_t num_points) const
Array shape for tabulate basis values and derivatives at set of points.
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:1219
polyset::type polyset_type() const
Get the element polyset type.
Definition finite-element.h:490
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & wcoeffs() const
Get the coefficients that define the polynomial set in terms of the orthonormal polynomials.
Definition finite-element.h:1271
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
Map function values from a physical cell to the reference.
Definition finite-element.cpp:1741
bool dof_transformations_are_permutations() const
Indicates if the degree-of-freedom transformations are all permutations.
Definition finite-element.h:551
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
Map function values from the reference to a physical cell.
Definition finite-element.cpp:1701
cell::type cell_type() const
Get the element cell type.
Definition finite-element.h:486
bool interpolation_is_identity() const
Indicates whether or not the interpolation matrix for this element is an identity matrix.
Definition finite-element.h:1378
bool discontinuous() const
Indicates whether this element is the discontinuous variant.
Definition finite-element.h:547
std::vector< std::vector< FiniteElement< F > > > get_tensor_product_representation() const
Get the tensor product representation of this element.
Definition finite-element.h:1367
maps::type map_type() const
Map type for the element.
Definition finite-element.h:537
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & dual_matrix() const
Get the dual matrix.
Definition finite-element.h:1230
element::dpc_variant dpc_variant() const
DPC variant of the element.
Definition finite-element.h:533
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & coefficient_matrix() const
Get the matrix of coefficients.
Definition finite-element.h:1334
FiniteElement & operator=(const FiniteElement &element)=default
Assignment operator.
element::family family() const
The finite element family.
Definition finite-element.h:522
~FiniteElement()=default
Destructor.
std::function< void(O &, const P &, const Q &, F, const R &)> map_fn() const
Return a function that performs the appropriate push-forward/pull-back for the element type.
Definition finite-element.h:623
type
Cell type.
Definition cell.h:21
int topological_dimension(cell::type celltype)
Definition cell.cpp:295
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:521
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:61
void contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Contravariant Piola map.
Definition maps.h:81
void double_contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Double contravariant Piola map.
Definition maps.h:135
void double_covariant_piola(O &&r, const P &U, const Q &J, double, const R &K)
Double covariant Piola map.
Definition maps.h:103
type
Map type.
Definition maps.h:39
type
Cell type.
Definition polyset.h:136
void apply_permutation(std::span< const std::size_t > perm, std::span< E > data, std::size_t offset=0, std::size_t n=1)
Apply a (precomputed) permutation .
Definition precompute.h:144
void apply_permutation_mapped(std::span< const std::size_t > perm, std::span< E > data, std::span< const int > emap, std::size_t n=1)
Permutation of mapped data.
Definition precompute.h:154
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:194
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:400
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:350
std::string version()
Definition finite-element.cpp:1779
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:610
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:331