8 #include "element-families.h"
11 #include "precompute.h"
28 = std::experimental::mdspan<double,
29 std::experimental::dextents<std::size_t, 2>>;
31 = std::experimental::mdspan<double,
32 std::experimental::dextents<std::size_t, 4>>;
34 = std::experimental::mdspan<
const double,
35 std::experimental::dextents<std::size_t, 2>>;
37 = std::experimental::mdspan<
const double,
38 std::experimental::dextents<std::size_t, 3>>;
40 = std::experimental::mdspan<
const double,
41 std::experimental::dextents<std::size_t, 4>>;
44 = std::experimental::mdarray<double,
45 std::experimental::dextents<std::size_t, 2>>;
47 = std::experimental::mdarray<double,
48 std::experimental::dextents<std::size_t, 4>>;
52 inline std::array<std::vector<cmdspan2_t>, 4>
53 to_mdspan(std::array<std::vector<mdarray2_t>, 4>& x)
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());
65 inline std::array<std::vector<cmdspan4_t>, 4>
66 to_mdspan(std::array<std::vector<mdarray4_t>, 4>& M)
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());
78 inline std::array<std::vector<cmdspan2_t>, 4>
79 to_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)
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]));
92 inline std::array<std::vector<cmdspan4_t>, 4>
93 to_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)
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]));
128 std::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>>
133 const std::array<std::vector<cmdspan4_t>, 4>& M,
134 std::size_t tdim, std::size_t value_size);
146 = std::experimental::mdspan<
const double,
147 std::experimental::dextents<std::size_t, 2>>;
149 = std::experimental::mdspan<
const double,
150 std::experimental::dextents<std::size_t, 4>>;
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>>>
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>>>
351 const std::array<std::vector<cmdspan2_t>, 4>&
x,
352 const std::array<std::vector<cmdspan4_t>, 4>&
M,
356 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
364 const std::array<std::vector<cmdspan2_t>, 4>&
x,
365 const std::array<std::vector<cmdspan4_t>, 4>&
M,
368 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
403 std::size_t num_points)
const;
426 std::pair<std::vector<double>, std::array<std::size_t, 4>>
427 tabulate(
int nd, impl::cmdspan2_t
x)
const;
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;
481 void tabulate(
int nd, impl::cmdspan2_t
x, impl::mdspan4_t basis)
const;
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;
534 const std::vector<std::size_t>&
value_shape()
const;
583 std::pair<std::vector<double>, std::array<std::size_t, 3>>
585 std::span<const double> detJ, impl::cmdspan3_t K)
const;
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;
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
634 case maps::type::identity:
635 return [](O& u,
const P& U,
const Q&, double,
const R&)
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)
643 case maps::type::covariantPiola:
644 return [](O& u,
const P& U,
const Q& J,
double detJ,
const R& K)
646 case maps::type::contravariantPiola:
647 return [](O& u,
const P& U,
const Q& J,
double detJ,
const R& K)
649 case maps::type::doubleCovariantPiola:
650 return [](O& u,
const P& U,
const Q& J,
double detJ,
const R& K)
652 case maps::type::doubleContravariantPiola:
653 return [](O& u,
const P& U,
const Q& J,
double detJ,
const R& K)
656 throw std::runtime_error(
"Map not implemented");
666 const std::vector<std::vector<std::vector<int>>>&
entity_dofs()
const;
758 std::pair<std::vector<double>, std::array<std::size_t, 3>>
766 std::pair<std::vector<double>, std::array<std::size_t, 3>>>
777 std::uint32_t cell_info)
const;
787 std::uint32_t cell_info)
const;
797 template <
typename T>
799 std::uint32_t cell_info)
const;
809 template <
typename T>
812 std::uint32_t cell_info)
const;
822 template <
typename T>
824 const std::span<T>& data,
int block_size, std::uint32_t cell_info)
const;
834 template <
typename T>
837 std::uint32_t cell_info)
const;
847 template <
typename T>
850 std::uint32_t cell_info)
const;
860 template <
typename T>
862 const std::span<T>& data,
int block_size, std::uint32_t cell_info)
const;
873 template <
typename T>
875 const std::span<T>& data,
int block_size, std::uint32_t cell_info)
const;
885 template <
typename T>
887 const std::span<T>& data,
int block_size, std::uint32_t cell_info)
const;
893 const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
948 const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
956 const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
993 const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
1001 std::vector<std::pair<std::vector<double>, std::array<std::size_t, 2>>>,
1042 std::vector<std::pair<std::vector<double>, std::array<std::size_t, 4>>>,
1051 const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
1074 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
1089 std::size_t _cell_tdim;
1092 std::vector<std::vector<cell::type>> _cell_subentity_types;
1107 int _interpolation_nderivs;
1110 int _highest_degree;
1113 int _highest_complete_degree;
1116 std::vector<std::size_t> _value_shape;
1126 std::pair<std::vector<double>, std::array<std::size_t, 2>> _coeffs;
1129 std::vector<std::vector<std::vector<int>>> _edofs;
1132 std::vector<std::vector<std::vector<int>>> _e_closure_dofs;
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>>;
1138 std::map<cell::type, array3_t> _entity_transformations;
1145 std::pair<std::vector<double>, std::array<std::size_t, 2>> _points;
1150 std::vector<std::pair<std::vector<double>, std::array<std::size_t, 2>>>,
1155 std::pair<std::vector<double>, std::array<std::size_t, 2>> _matM;
1159 bool _dof_transformations_are_permutations;
1162 bool _dof_transformations_are_identity;
1167 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm;
1172 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_rev;
1175 = std::vector<std::pair<std::vector<std::size_t>, array2_t>>;
1178 std::map<cell::type, trans_data_t> _etrans;
1181 std::map<cell::type, trans_data_t> _etransT;
1184 std::map<cell::type, trans_data_t> _etrans_inv;
1187 std::map<cell::type, trans_data_t> _etrans_invT;
1191 bool _discontinuous;
1194 std::pair<std::vector<double>, std::array<std::size_t, 2>> _dual_matrix;
1201 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
1205 bool _interpolation_is_identity;
1209 std::pair<std::vector<double>, std::array<std::size_t, 2>> _wcoeffs;
1213 = std::vector<std::pair<std::vector<double>, std::array<std::size_t, 4>>>;
1214 std::array<array4_t, 4> _M;
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);
1262 bool discontinuous);
1289 bool discontinuous);
1302 int degree,
bool discontinuous);
1351 template <
typename T>
1354 std::uint32_t cell_info)
const
1356 if (_dof_transformations_are_identity)
1359 if (_cell_tdim >= 2)
1363 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1365 for (
auto& edofs0 : _edofs[0])
1366 dofstart += edofs0.size();
1370 auto& [v_size_t, matrix] = _etrans.at(cell::type::interval)[0];
1371 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1374 if (cell_info >> (face_start + e) & 1)
1377 std::span(v_size_t),
1378 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1381 dofstart += _edofs[1][e].size();
1385 if (_cell_tdim == 3)
1388 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1390 auto& trans = _etrans.at(_cell_subentity_types[2][f]);
1393 if (cell_info >> (3 * f) & 1)
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,
1405 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
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,
1416 dofstart += _edofs[2][f].size();
1422 template <
typename T>
1424 const std::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1426 if (_dof_transformations_are_identity)
1429 if (_cell_tdim >= 2)
1433 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1435 for (
auto& edofs0 : _edofs[0])
1436 dofstart += edofs0.size();
1440 auto& [v_size_t, matrix] = _etransT.at(cell::type::interval)[0];
1441 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1444 if (cell_info >> (face_start + e) & 1)
1447 std::span(v_size_t),
1448 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1451 dofstart += _edofs[1][e].size();
1455 if (_cell_tdim == 3)
1458 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1460 auto& trans = _etransT.at(_cell_subentity_types[2][f]);
1463 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
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,
1475 if (cell_info >> (3 * f) & 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,
1485 dofstart += _edofs[2][f].size();
1491 template <
typename T>
1493 const std::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1495 if (_dof_transformations_are_identity)
1498 if (_cell_tdim >= 2)
1502 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1504 for (
auto& edofs0 : _edofs[0])
1505 dofstart += edofs0.size();
1508 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1511 if (cell_info >> (face_start + e) & 1)
1513 auto& [v_size_t, matrix] = _etrans_invT.at(cell::type::interval)[0];
1515 cmdspan2_t(matrix.first.data(), matrix.second),
1516 data, dofstart, block_size);
1518 dofstart += _edofs[1][e].size();
1521 if (_cell_tdim == 3)
1524 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1526 auto& trans = _etrans_invT.at(_cell_subentity_types[2][f]);
1529 if (cell_info >> (3 * f) & 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,
1541 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
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,
1552 dofstart += _edofs[2][f].size();
1558 template <
typename T>
1560 const std::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1562 if (_dof_transformations_are_identity)
1565 if (_cell_tdim >= 2)
1569 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1571 for (
auto& edofs0 : _edofs[0])
1572 dofstart += edofs0.size();
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)
1580 if (cell_info >> (face_start + e) & 1)
1583 std::span(v_size_t),
1584 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1587 dofstart += _edofs[1][e].size();
1591 if (_cell_tdim == 3)
1594 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1596 auto& trans = _etrans_inv.at(_cell_subentity_types[2][f]);
1599 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
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,
1611 if (cell_info >> (3 * f) & 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,
1622 dofstart += _edofs[2][f].size();
1628 template <
typename T>
1630 const std::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1632 if (_dof_transformations_are_identity)
1635 if (_cell_tdim >= 2)
1639 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1641 for (
auto& edofs0 : _edofs[0])
1642 dofstart += edofs0.size();
1646 auto& [v_size_t, matrix] = _etrans.at(cell::type::interval)[0];
1647 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1650 if (cell_info >> (face_start + e) & 1)
1653 std::span(v_size_t),
1654 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1658 dofstart += _edofs[1][e].size();
1662 if (_cell_tdim == 3)
1665 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1667 auto& trans = _etrans.at(_cell_subentity_types[2][f]);
1670 if (cell_info >> (3 * f) & 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,
1682 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
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,
1692 dofstart += _edofs[2][f].size();
1698 template <
typename T>
1700 const std::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1702 if (_dof_transformations_are_identity)
1705 if (_cell_tdim >= 2)
1709 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1711 for (
auto& edofs0 : _edofs[0])
1712 dofstart += edofs0.size();
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)
1720 if (cell_info >> (face_start + e) & 1)
1723 std::span(v_size_t),
1724 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1727 dofstart += _edofs[1][e].size();
1731 if (_cell_tdim == 3)
1734 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1736 auto& trans = _etrans_invT.at(_cell_subentity_types[2][f]);
1739 if (cell_info >> (3 * f) & 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,
1751 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
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,
1761 dofstart += _edofs[2][f].size();
1767 template <
typename T>
1769 const std::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1771 if (_dof_transformations_are_identity)
1774 if (_cell_tdim >= 2)
1778 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1780 for (
auto& edofs0 : _edofs[0])
1781 dofstart += edofs0.size();
1785 auto& [v_size_t, matrix] = _etransT.at(cell::type::interval)[0];
1786 for (std::size_t e = 0; e < _edofs[1].size(); ++e)
1789 if (cell_info >> (face_start + e) & 1)
1792 std::span(v_size_t),
1793 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1796 dofstart += _edofs[1][e].size();
1800 if (_cell_tdim == 3)
1803 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1805 auto& trans = _etransT.at(_cell_subentity_types[2][f]);
1808 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
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,
1820 if (cell_info >> (3 * f) & 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,
1830 dofstart += _edofs[2][f].size();
1836 template <
typename T>
1838 const std::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1840 if (_dof_transformations_are_identity)
1843 if (_cell_tdim >= 2)
1847 int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
1849 for (
auto& edofs0 : _edofs[0])
1850 dofstart += edofs0.size();
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)
1858 if (cell_info >> (face_start + e) & 1)
1861 std::span(v_size_t),
1862 cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
1865 dofstart += _edofs[1][e].size();
1869 if (_cell_tdim == 3)
1872 for (std::size_t f = 0; f < _edofs[2].size(); ++f)
1874 auto& trans = _etrans_inv.at(_cell_subentity_types[2][f]);
1877 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
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,
1889 if (cell_info >> (3 * f) & 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,
1900 dofstart += _edofs[2][f].size();