8 #include "element-families.h"
10 #include "precompute.h"
16 #include <xtensor/xtensor.hpp>
17 #include <xtensor/xview.hpp>
18 #include <xtl/xspan.hpp>
39 std::tuple<std::array<std::vector<xt::xtensor<double, 2>>, 4>,
40 std::array<std::vector<xt::xtensor<double, 3>>, 4>>
42 const std::array<std::vector<xt::xtensor<double, 3>>, 4>& M,
43 int tdim,
int value_size);
222 const xt::xtensor<double, 2>&
wcoeffs,
223 const std::array<std::vector<xt::xtensor<double, 2>>, 4>&
x,
224 const std::array<std::vector<xt::xtensor<double, 3>>, 4>&
M,
226 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
263 std::size_t num_points)
const;
286 xt::xtensor<double, 4>
tabulate(
int nd,
const xt::xarray<double>&
x)
const;
313 void tabulate(
int nd,
const xt::xarray<double>&
x,
314 xt::xtensor<double, 4>& basis)
const;
376 xt::xtensor<double, 3>
push_forward(
const xt::xtensor<double, 3>& U,
377 const xt::xtensor<double, 3>& J,
378 const xtl::span<const double>& detJ,
379 const xt::xtensor<double, 3>& K)
const;
387 xt::xtensor<double, 3>
pull_back(
const xt::xtensor<double, 3>& u,
388 const xt::xtensor<double, 3>& J,
389 const xtl::span<const double>& detJ,
390 const xt::xtensor<double, 3>& K)
const;
423 template <
typename O,
typename P,
typename Q,
typename R>
424 std::function<void(O&,
const P&,
const Q&,
double,
const R&)>
map_fn()
const
428 case maps::type::identity:
429 return [](O& u,
const P& U,
const Q&, double,
const R&) { u.assign(U); };
430 case maps::type::L2Piola:
431 return [](O& u,
const P& U,
const Q& J,
double detJ,
const R& K) {
434 case maps::type::covariantPiola:
435 return [](O& u,
const P& U,
const Q& J,
double detJ,
const R& K) {
438 case maps::type::contravariantPiola:
439 return [](O& u,
const P& U,
const Q& J,
double detJ,
const R& K) {
442 case maps::type::doubleCovariantPiola:
443 return [](O& u,
const P& U,
const Q& J,
double detJ,
const R& K) {
446 case maps::type::doubleContravariantPiola:
447 return [](O& u,
const P& U,
const Q& J,
double detJ,
const R& K) {
451 throw std::runtime_error(
"Map not implemented");
482 const std::vector<std::vector<std::vector<int>>>&
entity_dofs()
const;
585 std::uint32_t cell_info)
const;
595 std::uint32_t cell_info)
const;
605 template <
typename T>
607 std::uint32_t cell_info)
const;
617 template <
typename T>
620 std::uint32_t cell_info)
const;
630 template <
typename T>
632 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const;
642 template <
typename T>
645 std::uint32_t cell_info)
const;
655 template <
typename T>
658 std::uint32_t cell_info)
const;
668 template <
typename T>
670 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const;
680 template <
typename T>
682 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const;
692 template <
typename T>
694 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const;
700 const xt::xtensor<double, 2>&
points()
const;
791 const xt::xtensor<double, 2>&
wcoeffs()
const;
796 const std::array<std::vector<xt::xtensor<double, 2>>, 4>&
x()
const;
830 const std::array<std::vector<xt::xtensor<double, 3>>, 4>&
M()
const;
861 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
873 std::size_t _cell_tdim;
876 std::vector<std::vector<cell::type>> _cell_subentity_types;
891 std::vector<int> _value_shape;
901 xt::xtensor<double, 2> _coeffs;
909 std::vector<std::vector<int>> _num_edofs;
913 std::vector<std::vector<int>> _num_e_closure_dofs;
916 std::vector<std::vector<std::vector<int>>> _edofs;
919 std::vector<std::vector<std::vector<int>>> _e_closure_dofs;
922 std::map<cell::type, xt::xtensor<double, 3>> _entity_transformations;
929 xt::xtensor<double, 2> _points;
933 std::array<std::vector<xt::xtensor<double, 2>>, 4> _x;
936 xt::xtensor<double, 2> _matM;
940 bool _dof_transformations_are_permutations;
943 bool _dof_transformations_are_identity;
948 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm;
953 std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_rev;
957 std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
958 xt::xtensor<double, 2>>>>
963 std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
964 xt::xtensor<double, 2>>>>
969 std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
970 xt::xtensor<double, 2>>>>
975 std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
976 xt::xtensor<double, 2>>>>
984 xt::xtensor<double, 2> _dual_matrix;
990 std::array<int, 2> _degree_bounds;
997 std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
1001 bool _interpolation_is_identity;
1005 xt::xtensor<double, 2> _wcoeffs;
1008 std::array<std::vector<xt::xtensor<double, 3>>, 4> _M;
1032 const std::vector<std::size_t>& value_shape,
1033 const xt::xtensor<double, 2>& wcoeffs,
1034 const std::array<std::vector<xt::xtensor<double, 2>>, 4>& x,
1035 const std::array<std::vector<xt::xtensor<double, 3>>, 4>& M,
1036 maps::type map_type,
bool discontinuous,
int highest_complete_degree);
1049 bool discontinuous);
1076 bool discontinuous);
1089 int degree,
bool discontinuous);
1138 template <
typename T>
1141 std::uint32_t cell_info)
const
1143 if (_dof_transformations_are_identity)
1146 if (_cell_tdim >= 2)
1150 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1152 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1155 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1158 if (cell_info >> (face_start + e) & 1)
1160 dofstart, block_size);
1161 dofstart += _num_edofs[1][e];
1164 if (_cell_tdim == 3)
1167 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1170 if (cell_info >> (3 * f) & 1)
1172 data, dofstart, block_size);
1175 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1177 data, dofstart, block_size);
1178 dofstart += _num_edofs[2][f];
1184 template <
typename T>
1186 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1188 if (_dof_transformations_are_identity)
1191 if (_cell_tdim >= 2)
1195 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1197 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1200 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1203 if (cell_info >> (face_start + e) & 1)
1205 dofstart, block_size);
1206 dofstart += _num_edofs[1][e];
1209 if (_cell_tdim == 3)
1212 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1215 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1217 data, dofstart, block_size);
1219 if (cell_info >> (3 * f) & 1)
1221 data, dofstart, block_size);
1222 dofstart += _num_edofs[2][f];
1228 template <
typename T>
1230 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1232 if (_dof_transformations_are_identity)
1235 if (_cell_tdim >= 2)
1239 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1241 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1244 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1247 if (cell_info >> (face_start + e) & 1)
1249 dofstart, block_size);
1250 dofstart += _num_edofs[1][e];
1253 if (_cell_tdim == 3)
1256 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1259 if (cell_info >> (3 * f) & 1)
1261 _etrans_invT.at(_cell_subentity_types[2][f])[1], data, dofstart,
1265 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1267 _etrans_invT.at(_cell_subentity_types[2][f])[0], data, dofstart,
1269 dofstart += _num_edofs[2][f];
1275 template <
typename T>
1277 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1279 if (_dof_transformations_are_identity)
1282 if (_cell_tdim >= 2)
1286 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1288 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1291 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1294 if (cell_info >> (face_start + e) & 1)
1296 dofstart, block_size);
1297 dofstart += _num_edofs[1][e];
1300 if (_cell_tdim == 3)
1303 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1306 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1308 _etrans_inv.at(_cell_subentity_types[2][f])[0], data, dofstart,
1311 if (cell_info >> (3 * f) & 1)
1313 _etrans_inv.at(_cell_subentity_types[2][f])[1], data, dofstart,
1315 dofstart += _num_edofs[2][f];
1321 template <
typename T>
1323 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1325 if (_dof_transformations_are_identity)
1328 if (_cell_tdim >= 2)
1332 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1334 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1337 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1340 if (cell_info >> (face_start + e) & 1)
1342 _etrans.at(cell::type::interval)[0], data, dofstart, block_size);
1343 dofstart += _num_edofs[1][e];
1346 if (_cell_tdim == 3)
1349 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1352 if (cell_info >> (3 * f) & 1)
1354 _etrans.at(_cell_subentity_types[2][f])[1], data, dofstart,
1358 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1360 _etrans.at(_cell_subentity_types[2][f])[0], data, dofstart,
1362 dofstart += _num_edofs[2][f];
1368 template <
typename T>
1370 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1372 if (_dof_transformations_are_identity)
1375 if (_cell_tdim >= 2)
1379 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1381 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1384 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1387 if (cell_info >> (face_start + e) & 1)
1389 _etrans_invT.at(cell::type::interval)[0], data, dofstart,
1391 dofstart += _num_edofs[1][e];
1394 if (_cell_tdim == 3)
1397 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1400 if (cell_info >> (3 * f) & 1)
1402 _etrans_invT.at(_cell_subentity_types[2][f])[1], data, dofstart,
1406 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1408 _etrans_invT.at(_cell_subentity_types[2][f])[0], data, dofstart,
1410 dofstart += _num_edofs[2][f];
1416 template <
typename T>
1418 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1420 if (_dof_transformations_are_identity)
1423 if (_cell_tdim >= 2)
1427 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1429 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1432 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1435 if (cell_info >> (face_start + e) & 1)
1437 _etransT.at(cell::type::interval)[0], data, dofstart, block_size);
1438 dofstart += _num_edofs[1][e];
1441 if (_cell_tdim == 3)
1444 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1447 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1449 _etransT.at(_cell_subentity_types[2][f])[0], data, dofstart,
1453 if (cell_info >> (3 * f) & 1)
1455 _etransT.at(_cell_subentity_types[2][f])[1], data, dofstart,
1458 dofstart += _num_edofs[2][f];
1464 template <
typename T>
1466 const xtl::span<T>& data,
int block_size, std::uint32_t cell_info)
const
1468 if (_dof_transformations_are_identity)
1471 if (_cell_tdim >= 2)
1475 int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1477 = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1480 for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1483 if (cell_info >> (face_start + e) & 1)
1485 _etrans_inv.at(cell::type::interval)[0], data, dofstart,
1487 dofstart += _num_edofs[1][e];
1490 if (_cell_tdim == 3)
1493 for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1496 for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1498 _etrans_inv.at(_cell_subentity_types[2][f])[0], data, dofstart,
1502 if (cell_info >> (3 * f) & 1)
1504 _etrans_inv.at(_cell_subentity_types[2][f])[1], data, dofstart,
1507 dofstart += _num_edofs[2][f];