deal.II version 9.8.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
reference_cell.h
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2020 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13#ifndef dealii_tria_reference_cell_h
14#define dealii_tria_reference_cell_h
15
16#include <deal.II/base/config.h>
17
21#include <deal.II/base/point.h>
22#include <deal.II/base/tensor.h>
23#include <deal.II/base/types.h>
25
27
28#include <boost/container/small_vector.hpp>
29
30#include <iosfwd>
31#include <string>
32
34
35// Forward declarations
36#ifndef DOXYGEN
37template <int dim, int spacedim>
38class Mapping;
39
40template <int dim>
41class Quadrature;
42
43template <int dim>
44class ReferenceCell;
45#endif
46
47
48namespace internal
49{
62 template <int dim>
63 constexpr ReferenceCell<dim>
64 make_reference_cell_from_int(const std::uint8_t kind);
65} // namespace internal
66
67
68
75enum class IsotropicRefinementChoice : std::uint8_t
76{
93};
94
95
96
127template <int dim>
129{
130public:
140 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
141 "Use ReferenceCells::n_vertices_to_reference_cell() instead")
142 static ReferenceCell
143 n_vertices_to_type(const int dim_, const unsigned int n_vertices);
144
153 constexpr ReferenceCell();
154
159
165 constexpr bool
167
171 constexpr bool
172 is_simplex() const;
173
179 constexpr unsigned int
181
185
190
201 double
202 d_linear_shape_function(const Point<dim> &xi, const unsigned int i) const;
203
210 const unsigned int i) const;
211
221 template <int spacedim = dim>
222 std::unique_ptr<Mapping<dim, spacedim>>
223 get_default_mapping(const unsigned int degree) const;
224
236 template <int spacedim = dim>
239
249 get_gauss_type_quadrature(const unsigned n_points_1d) const;
250
262
276 const Quadrature<dim> &
278
282
287
293 constexpr unsigned int
294 n_vertices() const;
295
300 std_cxx20::ranges::iota_view<unsigned int, unsigned int>
302
311 vertex(const unsigned int v) const;
312
318 constexpr unsigned int
319 n_lines() const;
320
325 std_cxx20::ranges::iota_view<unsigned int, unsigned int>
327
333 constexpr unsigned int
334 n_faces() const;
335
340 std_cxx20::ranges::iota_view<unsigned int, unsigned int>
342
348 std_cxx20::ranges::iota_view<unsigned int, unsigned int>
349 face_indices_by_type(const ReferenceCell<dim - 1> &face_ref_type) const;
350
354 constexpr unsigned int
356
361 get_isotropic_refinement_choice(const unsigned int ref_choice) const;
362
374 constexpr unsigned int
376
381 unsigned int
384
389 std_cxx20::ranges::iota_view<unsigned int, unsigned int>
391
411 std::pair<unsigned int, RefinementCase<dim - 1>>
413 const types::geometric_orientation combined_face_orientation,
414 const internal::SubfaceCase<dim> subface_case,
415 const unsigned int subface_no) const;
416
426
439 constexpr ReferenceCell<dim - 1>
440 face_reference_cell(const unsigned int face_index) const;
441
445
450
458 "Use numbers::default_geometric_orientation instead.")
459 static constexpr types::geometric_orientation
461
471 "Use numbers::reverse_line_orientation instead.")
472 static constexpr types::geometric_orientation
474
483 "Use the version of this function which takes the orientation as the third "
484 "parameter instead.")
485 unsigned int
486 child_cell_on_face(const unsigned int face, const unsigned int subface) const;
487
528 unsigned int
530 const unsigned int face,
531 const unsigned int subface,
532 const types::geometric_orientation combined_orientation) const;
533
543 std::array<unsigned int, 2>
545
555 std::array<unsigned int, 2>
556 standard_line_to_face_and_line_index(const unsigned int line) const;
557
566 unsigned int
567 line_to_cell_vertices(const unsigned int line,
568 const unsigned int vertex) const;
569
588 unsigned int
589 face_to_cell_lines(const unsigned int face,
590 const unsigned int line,
591 const types::geometric_orientation face_orientation) const;
592
609 unsigned int
611 const unsigned int face,
612 const unsigned int vertex,
613 const types::geometric_orientation face_orientation) const;
614
635 unsigned int
637 const unsigned int face,
638 const unsigned int line,
639 const unsigned int vertex,
640 const types::geometric_orientation face_orientation) const;
641
661 Point<dim>
662 face_vertex_location(const unsigned int face,
663 const unsigned int vertex) const;
664
678 Point<dim>
680 const unsigned int face_no,
681 const unsigned int subface_no,
682 const unsigned int subface_vertex_no,
683 const RefinementCase<dim - 1> face_refinement_case) const;
684
689 unsigned int
691 const unsigned int vertex,
692 const unsigned int face,
693 const types::geometric_orientation face_orientation) const;
694
699 unsigned int
701 const unsigned int line,
702 const unsigned int face,
703 const types::geometric_orientation face_orientation) const;
704
713 "Use face_to_cell_line_orientation() instead.")
714 types::geometric_orientation
716 const unsigned int line,
717 const unsigned int face,
718 const types::geometric_orientation face_orientation,
719 const types::geometric_orientation line_orientation) const;
720
743 types::geometric_orientation
745 const unsigned int face_line_no,
746 const unsigned int face_no,
747 const types::geometric_orientation face_orientation,
748 const types::geometric_orientation line_orientation) const;
749
775 unsigned int
776 opposite_face_index(const unsigned int face_no) const;
777
781
787
801 double
802 volume() const;
803
819 double
820 face_measure(const unsigned int face_no) const;
821
834 Point<dim>
835 barycenter() const;
836
856 bool
857 contains_point(const Point<dim> &p, const double tolerance = 0) const;
858
863 Point<dim>
864 closest_point(const Point<dim> &p) const;
865
873 Tensor<1, dim>
874 face_tangent_vector(const unsigned int face_no, const unsigned int i) const;
875
885 DEAL_II_DEPRECATED_WITH_COMMENT("Use face_tangent_vector() instead.")
886 Tensor<1, dim>
887 unit_tangential_vectors(const unsigned int face_no,
888 const unsigned int i) const;
889
893 Tensor<1, dim>
894 face_normal_vector(const unsigned int face_no) const;
895
901 DEAL_II_DEPRECATED_WITH_COMMENT("Use face_normal_vector() instead.")
902 Tensor<1, dim>
903 unit_normal_vectors(const unsigned int face_no) const;
904
910 unsigned int
911 n_face_orientations(const unsigned int face_no) const;
912
929 template <typename T, std::size_t N>
930 DEAL_II_DEPRECATED types::geometric_orientation
931 compute_orientation(const std::array<T, N> &vertices_0,
932 const std::array<T, N> &vertices_1) const;
933
958 template <typename T>
959 types::geometric_orientation
960 get_combined_orientation(const ArrayView<const T> &vertices_0,
961 const ArrayView<const T> &vertices_1) const;
962
963
977 template <typename T, std::size_t N>
978 DEAL_II_DEPRECATED std::array<T, N>
979 permute_according_orientation(const std::array<T, N> &vertices,
980 const unsigned int orientation) const;
981
993 template <typename T>
994 boost::container::small_vector<T, 8>
996 const ArrayView<const T> &vertices,
997 const types::geometric_orientation orientation) const;
998
1004 types::geometric_orientation
1006 const types::geometric_orientation orientation) const;
1007
1011 ArrayView<const unsigned int>
1012 faces_for_given_vertex(const unsigned int vertex_index) const;
1013
1025 constexpr ::ndarray<unsigned int, 13, 4>
1026 new_isotropic_child_face_lines(const unsigned int refinement_choice) const;
1027
1041 constexpr ::ndarray<unsigned int, 13, 4, 2>
1043 const unsigned int refinement_choice) const;
1044
1056 constexpr ::ndarray<unsigned int, 10, 6>
1057 new_isotropic_child_cell_faces(const unsigned int refinement_choice) const;
1058
1070 constexpr ::ndarray<unsigned int, 10, 8>
1071 new_isotropic_child_cell_vertices(const unsigned int refinement_choice) const;
1072
1076
1081
1085 unsigned int
1086 exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const;
1087
1091 unsigned int
1092 exodusii_face_to_deal_face(const unsigned int face_n) const;
1093
1097 unsigned int
1098 unv_vertex_to_deal_vertex(const unsigned int vertex_n) const;
1099
1103 unsigned int
1104 vtk_linear_type() const;
1105
1110 unsigned int
1111 vtk_quadratic_type() const;
1112
1117 unsigned int
1118 vtk_lagrange_type() const;
1119
1132 unsigned int
1134 const std::array<unsigned, dim> &node_indices,
1135 const std::array<unsigned, dim> &nodes_per_direction,
1136 const bool legacy_format) const;
1137
1141 unsigned int
1142 vtk_vertex_to_deal_vertex(const unsigned int vertex_index) const;
1143
1147 unsigned int
1148 gmsh_element_type() const;
1149
1153
1158
1163 std::string
1164 to_string() const;
1165
1169 constexpr operator std::uint8_t() const;
1170
1174 constexpr bool
1175 operator==(const ReferenceCell<dim> &type) const;
1176
1180 constexpr bool
1181 operator!=(const ReferenceCell<dim> &type) const;
1182
1188 template <class Archive>
1189 void
1190 serialize(Archive &archive, const unsigned int /*version*/);
1191
1195 static constexpr std::size_t
1197
1201
1202private:
1207
1214 constexpr ReferenceCell(const std::uint8_t kind);
1215
1221 static constexpr ndarray<unsigned int, 1, 1> vertex_vertex_permutations = {
1222 {{{0}}}};
1223
1230
1241
1246 static constexpr ndarray<unsigned int, 8, 4>
1248 {{0, 1, 2, 3}},
1249 {{0, 2, 1, 3}},
1250 {{2, 0, 3, 1}},
1251 {{2, 3, 0, 1}},
1252 {{3, 2, 1, 0}},
1253 {{3, 1, 2, 0}},
1254 {{1, 3, 0, 2}},
1255 {{1, 0, 3, 2}},
1256 }};
1257
1262 {{{0, 1, 2}},
1263 {{2, 1, 0}},
1264 {{2, 0, 1}},
1265 {{1, 0, 2}},
1266 {{1, 2, 0}},
1267 {{0, 2, 1}}}};
1268
1274 {{
1275 {{0, 1, 2, 3}},
1276 {{2, 3, 0, 1}},
1277 {{3, 2, 0, 1}},
1278 {{0, 1, 3, 2}},
1279 {{1, 0, 3, 2}},
1280 {{3, 2, 1, 0}},
1281 {{2, 3, 1, 0}},
1282 {{1, 0, 2, 3}},
1283 }};
1284
1289 template <int dim_>
1290 friend constexpr ReferenceCell<dim_>
1292
1293 template <int dim_>
1294 friend std::ostream &
1295 operator<<(std::ostream &out, const ReferenceCell<dim_> &reference_cell);
1296
1297 template <int dim_>
1298 friend std::istream &
1299 operator>>(std::istream &in, ReferenceCell<dim_> &reference_cell);
1300};
1301
1302
1311template <int dim>
1312std::ostream &
1313operator<<(std::ostream &out, const ReferenceCell<dim> &reference_cell);
1314
1322template <int dim>
1323std::istream &
1324operator>>(std::istream &in, ReferenceCell<dim> &reference_cell);
1325
1326
1327
1328template <int dim>
1329inline constexpr ReferenceCell<dim>::ReferenceCell(const std::uint8_t kind)
1330 : kind(kind)
1331{}
1332
1333
1334
1335template <int dim>
1336inline constexpr ReferenceCell<dim>::operator std::uint8_t() const
1337{
1338 return kind;
1339}
1340
1341
1342
1343template <int dim>
1344inline constexpr bool
1346{
1347 return kind == type.kind;
1348}
1349
1350
1351
1352template <int dim>
1353inline constexpr bool
1355{
1356 return kind != type.kind;
1357}
1358
1359
1360
1361namespace internal
1362{
1363 template <int dim>
1364 inline constexpr ReferenceCell<dim>
1365 make_reference_cell_from_int(const std::uint8_t kind)
1366 {
1367#ifndef DEAL_II_CXX14_CONSTEXPR_BUG
1368 // Make sure these are the only indices from which objects can be
1369 // created.
1370 Assert((kind == std::numeric_limits<std::uint8_t>::max()) || (kind < 8),
1372#endif
1373
1374 // Call the private constructor, which we can from here because this
1375 // function is a 'friend'.
1376 return {kind};
1377 }
1378} // namespace internal
1379
1380
1381
1390{
1407
1408 template <int dim>
1411 std::numeric_limits<std::uint8_t>::max());
1412
1417 template <int dim>
1418 constexpr const std::vector<ReferenceCell<dim>>
1420
1426 template <int dim>
1427 constexpr const ReferenceCell<dim> &
1428 get_simplex();
1429
1435 template <int dim>
1436 constexpr const ReferenceCell<dim> &
1437 get_hypercube();
1438
1443 template <int dim>
1444 constexpr unsigned int
1446
1459 template <int dim>
1461 n_vertices_to_reference_cell(const unsigned int n_vertices);
1462
1470 template <int structdim>
1471 constexpr unsigned int
1473
1481 template <int structdim>
1482 constexpr unsigned int
1483 max_n_lines();
1484
1500 template <int structdim>
1501 constexpr unsigned int
1502 max_n_faces();
1503
1513 template <int structdim>
1514 constexpr unsigned int
1516} // namespace ReferenceCells
1517
1518
1519
1520template <int dim>
1522 : ReferenceCell(ReferenceCells::Invalid<dim>)
1523{}
1524
1525
1526
1527template <int dim>
1528template <class Archive>
1529inline void
1530ReferenceCell<dim>::serialize(Archive &archive, const unsigned int /*version*/)
1531{
1532 archive &kind;
1533}
1534
1535
1536
1537template <int dim>
1538inline constexpr std::size_t
1543
1544
1545
1546template <int dim>
1549{
1551 switch (this->kind)
1552 {
1554 return {&GeometryInfo<2>::vertex_to_face[vertex][0], 1};
1556 return {&GeometryInfo<2>::vertex_to_face[vertex][0], 2};
1558 {
1559 static constexpr ndarray<unsigned int, 3, 2> table = {
1560 {{{0, 2}}, {{0, 1}}, {{1, 2}}}};
1561 return table[vertex];
1562 }
1564 {
1565 static constexpr ndarray<unsigned int, 4, 3> table = {
1566 {{{0, 1, 2}}, {{0, 1, 3}}, {{0, 2, 3}}, {{1, 2, 3}}}};
1567
1568 return table[vertex];
1569 }
1571 {
1572 static constexpr unsigned int X = numbers::invalid_unsigned_int;
1573 static constexpr ndarray<unsigned int, 5, 4> table = {
1574 {{{0, 1, 3, X}},
1575 {{0, 2, 3, X}},
1576 {{0, 1, 4, X}},
1577 {{0, 2, 4, X}},
1578 {{1, 2, 3, 4}}}};
1579
1580 return {&table[vertex][0], vertex == 4 ? 4u : 3u};
1581 }
1583 {
1585 static constexpr ndarray<unsigned int, 6, 3> table = {{{{0, 2, 4}},
1586 {{0, 2, 3}},
1587 {{0, 3, 4}},
1588 {{1, 2, 4}},
1589 {{1, 2, 3}},
1590 {{1, 3, 4}}}};
1591
1592 return table[vertex];
1593 }
1595 return {&GeometryInfo<3>::vertex_to_face[vertex][0], 3};
1596 default:
1598 }
1599
1600 return {};
1601}
1602
1603
1604
1605template <int dim>
1606constexpr ::ndarray<unsigned int, 13, 4>
1608 const unsigned int refinement_choice) const
1609{
1611
1612 if constexpr (dim == 3)
1613 {
1614 const unsigned int X = numbers::invalid_unsigned_int;
1615
1616 switch (this->kind)
1617 {
1619 {
1620 switch (refinement_choice)
1621 {
1622 case 0:
1623 // new line is (6,8)
1624 return {{{{2, 3, 8, X}},
1625 {{0, 9, 5, X}},
1626 {{1, 6, 11, X}},
1627 {{4, 10, 7, X}},
1628 {{2, 12, 5, X}},
1629 {{1, 9, 12, X}},
1630 {{4, 8, 12, X}},
1631 {{6, 12, 10, X}},
1632 {{X, X, X, X}},
1633 {{X, X, X, X}},
1634 {{X, X, X, X}},
1635 {{X, X, X, X}},
1636 {{X, X, X, X}}}};
1637 case 1:
1638 // new line is (5,7)
1639 return {{{{2, 3, 8, X}},
1640 {{0, 9, 5, X}},
1641 {{1, 6, 11, X}},
1642 {{4, 10, 7, X}},
1643 {{0, 3, 12, X}},
1644 {{1, 12, 8, X}},
1645 {{4, 12, 9, X}},
1646 {{7, 11, 12, X}},
1647 {{X, X, X, X}},
1648 {{X, X, X, X}},
1649 {{X, X, X, X}},
1650 {{X, X, X, X}},
1651 {{X, X, X, X}}}};
1652 case 2:
1653 // new line is (4,9)
1654 return {{{{2, 3, 8, X}},
1655 {{0, 9, 5, X}},
1656 {{1, 6, 11, X}},
1657 {{4, 10, 7, X}},
1658 {{0, 12, 11, X}},
1659 {{2, 6, 12, X}},
1660 {{3, 12, 7, X}},
1661 {{5, 10, 12, X}},
1662 {{X, X, X, X}},
1663 {{X, X, X, X}},
1664 {{X, X, X, X}},
1665 {{X, X, X, X}},
1666 {{X, X, X, X}}}};
1667 default:
1669 }
1670 }
1671
1673 {
1674 constexpr ::ndarray<unsigned int, 13, 4>
1675 new_quad_lines_pyramid = {
1676 {{{0, 10, 16, X}}, // child 0
1677 {{2, 16, 6, X}},
1678 {{0, 17, 12, X}}, // child 1
1679 {{3, 7, 17, X}},
1680 {{1, 18, 15, X}}, // child 2
1681 {{2, 4, 18, X}},
1682 {{1, 13, 19, X}}, // child 3
1683 {{3, 19, 9, X}},
1684 {{5, 8, 11, 14}}, // child 9 (top)
1685 {{11, 17, 16, X}}, // child 8 (upside down)
1686 {{14, 18, 19, X}},
1687 {{5, 16, 18, X}},
1688 {{8, 19, 17, X}}}};
1689 return new_quad_lines_pyramid;
1690 }
1691
1693 {
1694 constexpr ::ndarray<unsigned int, 13, 4>
1695 new_quad_lines_wedge = {{{{2, 20, 11, X}}, // mid tri
1696 {{3, 6, 18, X}},
1697 {{19, 7, 10, X}},
1698 {{18, 19, 20, X}},
1699 {{4, 0, 14, 18}}, // lower
1700 {{0, 8, 12, 20}},
1701 {{8, 4, 13, 19}},
1702 {{5, 1, 18, 15}}, // upper
1703 {{1, 9, 20, 17}},
1704 {{9, 5, 19, 16}},
1705 {{X, X, X, X}},
1706 {{X, X, X, X}},
1707 {{X, X, X, X}}}};
1708 return new_quad_lines_wedge;
1709 }
1710
1712 {
1713 constexpr ::ndarray<unsigned int, 13, 4>
1714 new_quad_lines_hex = {{{{10, 28, 16, 24}},
1715 {{28, 14, 17, 25}},
1716 {{11, 29, 24, 20}},
1717 {{29, 15, 25, 21}},
1718 {{18, 26, 0, 28}},
1719 {{26, 22, 1, 29}},
1720 {{19, 27, 28, 4}},
1721 {{27, 23, 29, 5}},
1722 {{2, 24, 8, 26}},
1723 {{24, 6, 9, 27}},
1724 {{3, 25, 26, 12}},
1725 {{25, 7, 27, 13}},
1726 {{X, X, X, X}}}};
1727 return new_quad_lines_hex;
1728 }
1729
1730 default:
1732 }
1733 }
1734 else
1735 // We should never get here except in the 3d case.
1737
1738 return {};
1739}
1740
1741
1742
1743template <int dim>
1744constexpr ::ndarray<unsigned int, 13, 4, 2>
1746 const unsigned int refinement_choice) const
1747{
1749
1750 if constexpr (dim == 3)
1751 {
1752 const unsigned int X = numbers::invalid_unsigned_int;
1753
1754 switch (this->kind)
1755 {
1757 {
1758 switch (refinement_choice)
1759 {
1760 case 0:
1761 // new line is (6, 8)
1762 return {{{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
1763 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
1764 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
1765 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
1766 {{{{4, 6}}, {{6, 8}}, {{8, 4}}, {{X, X}}}},
1767 {{{{6, 5}}, {{5, 8}}, {{8, 6}}, {{X, X}}}},
1768 {{{{8, 7}}, {{7, 6}}, {{6, 8}}, {{X, X}}}},
1769 {{{{9, 6}}, {{6, 8}}, {{8, 9}}, {{X, X}}}},
1770 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1771 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1772 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1773 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1774 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}};
1775 case 1:
1776 // new line is (5, 7)
1777 return {{{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
1778 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
1779 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
1780 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
1781 {{{{5, 4}}, {{4, 7}}, {{7, 5}}, {{X, X}}}},
1782 {{{{6, 5}}, {{5, 7}}, {{7, 6}}, {{X, X}}}},
1783 {{{{8, 7}}, {{7, 5}}, {{5, 8}}, {{X, X}}}},
1784 {{{{7, 9}}, {{9, 5}}, {{5, 7}}, {{X, X}}}},
1785 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1786 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1787 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1788 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1789 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}};
1790 case 2:
1791 // new line is (4, 9)
1792 return {{{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
1793 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
1794 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
1795 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
1796 {{{{5, 4}}, {{4, 9}}, {{9, 5}}, {{X, X}}}},
1797 {{{{4, 6}}, {{6, 9}}, {{9, 4}}, {{X, X}}}},
1798 {{{{7, 4}}, {{4, 9}}, {{9, 7}}, {{X, X}}}},
1799 {{{{4, 8}}, {{8, 9}}, {{9, 4}}, {{X, X}}}},
1800 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1801 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1802 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1803 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1804 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}};
1805 default:
1807 }
1808 }
1809
1811 {
1812 constexpr ::ndarray<unsigned int, 13, 4, 2>
1813 quad_lines_vertices_pyramid = {
1814 {// child 0
1815 {{{{13, 7}}, {{7, 9}}, {{9, 13}}, {{X, X}}}},
1816 {{{{5, 13}}, {{13, 9}}, {{9, 5}}, {{X, X}}}},
1817 // child 1
1818 {{{{7, 13}}, {{13, 10}}, {{10, 7}}, {{X, X}}}},
1819 {{{{13, 6}}, {{6, 10}}, {{10, 13}}, {{X, X}}}},
1820 // child 2
1821 {{{{8, 13}}, {{13, 11}}, {{11, 8}}, {{X, X}}}},
1822 {{{{13, 5}}, {{5, 11}}, {{11, 13}}, {{X, X}}}},
1823 // child 3
1824 {{{{13, 8}}, {{8, 12}}, {{12, 13}}, {{X, X}}}},
1825 {{{{6, 13}}, {{13, 12}}, {{12, 6}}, {{X, X}}}},
1826 // child 9 (top)
1827 {{{{9, 11}}, {{10, 12}}, {{9, 10}}, {{11, 12}}}},
1828 // child 8 (upside down)
1829 {{{{9, 10}}, {{10, 13}}, {{13, 9}}, {{X, X}}}},
1830 {{{{12, 11}}, {{11, 13}}, {{13, 12}}, {{X, X}}}},
1831 {{{{11, 9}}, {{9, 13}}, {{13, 11}}, {{X, X}}}},
1832 {{{{10, 12}}, {{12, 13}}, {{13, 10}}, {{X, X}}}}}};
1833 return quad_lines_vertices_pyramid;
1834 }
1835
1837 {
1838 constexpr ::ndarray<unsigned int, 13, 4, 2>
1839 quad_lines_vertices_wedge = {
1840 {// mid tri
1841 {{{{12, 15}}, {{15, 17}}, {{17, 12}}, {{X, X}}}},
1842 {{{{15, 13}}, {{13, 16}}, {{16, 15}}, {{X, X}}}},
1843 {{{{17, 16}}, {{16, 14}}, {{14, 17}}, {{X, X}}}},
1844 {{{{15, 16}}, {{16, 17}}, {{17, 15}}, {{X, X}}}},
1845 // lower
1846 {{{{7, 16}}, {{6, 15}}, {{7, 6}}, {{16, 15}}}},
1847 {{{{6, 15}}, {{8, 17}}, {{6, 8}}, {{15, 17}}}},
1848 {{{{8, 17}}, {{7, 16}}, {{8, 7}}, {{17, 16}}}},
1849 // upper
1850 {{{{16, 10}}, {{15, 9}}, {{16, 15}}, {{10, 9}}}},
1851 {{{{15, 9}}, {{17, 11}}, {{15, 17}}, {{9, 11}}}},
1852 {{{{17, 11}}, {{16, 10}}, {{17, 16}}, {{11, 10}}}},
1853 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1854 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1855 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}};
1856 return quad_lines_vertices_wedge;
1857 }
1858
1860 {
1861 constexpr ::ndarray<unsigned int, 13, 4, 2>
1862 quad_lines_vertices_hex = {
1863 {{{{{10, 22}}, {{24, 26}}, {{10, 24}}, {{22, 26}}}},
1864 {{{{24, 26}}, {{11, 23}}, {{24, 11}}, {{26, 23}}}},
1865 {{{{22, 14}}, {{26, 25}}, {{22, 26}}, {{14, 25}}}},
1866 {{{{26, 25}}, {{23, 15}}, {{26, 23}}, {{25, 15}}}},
1867 {{{{8, 24}}, {{20, 26}}, {{8, 20}}, {{24, 26}}}},
1868 {{{{20, 26}}, {{12, 25}}, {{20, 12}}, {{26, 25}}}},
1869 {{{{24, 9}}, {{26, 21}}, {{24, 26}}, {{9, 21}}}},
1870 {{{{26, 21}}, {{25, 13}}, {{26, 25}}, {{21, 13}}}},
1871 {{{{16, 20}}, {{22, 26}}, {{16, 22}}, {{20, 26}}}},
1872 {{{{22, 26}}, {{17, 21}}, {{22, 17}}, {{26, 21}}}},
1873 {{{{20, 18}}, {{26, 23}}, {{20, 26}}, {{18, 23}}}},
1874 {{{{26, 23}}, {{21, 19}}, {{26, 21}}, {{23, 19}}}},
1875 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}};
1876 return quad_lines_vertices_hex;
1877 }
1878 default:
1880 }
1881 }
1882 else
1883 // We should never get here except in the 3d case.
1885
1886 return {};
1887}
1888
1889
1890
1891template <int dim>
1892constexpr ::ndarray<unsigned int, 10, 6>
1894 const unsigned int refinement_choice) const
1895{
1897
1898 if constexpr (dim == 3)
1899 {
1900 const unsigned int X = numbers::invalid_unsigned_int;
1901
1902 switch (this->kind)
1903 {
1905 {
1906 switch (refinement_choice)
1907 {
1908 case 0:
1909 // new line is (6, 8)
1910 return {{
1911 {{8, 13, 16, 0, X, X}},
1912 {{9, 12, 1, 21, X, X}},
1913 {{10, 2, 17, 20, X, X}},
1914 {{3, 14, 18, 22, X, X}},
1915 {{11, 1, 4, 5, X, X}},
1916 {{15, 0, 4, 6, X, X}},
1917 {{19, 7, 6, 3, X, X}},
1918 {{23, 5, 2, 7, X, X}},
1919 {{X, X, X, X, X, X}},
1920 {{X, X, X, X, X, X}},
1921 }};
1922 case 1:
1923 // new line is (5, 7)
1924 return {{
1925 {{8, 13, 16, 0, X, X}},
1926 {{9, 12, 1, 21, X, X}},
1927 {{10, 2, 17, 20, X, X}},
1928 {{3, 14, 18, 22, X, X}},
1929 {{11, 4, 0, 5, X, X}},
1930 {{15, 4, 1, 6, X, X}},
1931 {{19, 2, 5, 7, X, X}},
1932 {{23, 6, 7, 3, X, X}},
1933 {{X, X, X, X, X, X}},
1934 {{X, X, X, X, X, X}},
1935 }};
1936 case 2:
1937 // new line is (4, 9)
1938 return {{
1939 {{8, 13, 16, 0, X, X}},
1940 {{9, 12, 1, 21, X, X}},
1941 {{10, 2, 17, 20, X, X}},
1942 {{3, 14, 18, 22, X, X}},
1943 {{11, 4, 5, 2, X, X}},
1944 {{15, 6, 7, 3, X, X}},
1945 {{19, 5, 0, 6, X, X}},
1946 {{23, 1, 4, 7, X, X}},
1947 {{X, X, X, X, X, X}},
1948 {{X, X, X, X, X, X}},
1949 }};
1950 default:
1952 }
1953 }
1955 {
1956 constexpr ::ndarray<unsigned int, 10, 6>
1957 cell_quads_pyramid = {{
1958 {{13, 17, 0, 26, 1, X}}, // bottom pyramids
1959 {{14, 2, 22, 25, 3, X}}, //
1960 {{15, 18, 4, 5, 29, X}}, //
1961 {{16, 6, 21, 7, 30, X}}, //
1962 {{5, 1, 20, 11, X, X}}, // bottom wedges
1963 {{7, 3, 12, 24, X, X}}, //
1964 {{0, 2, 28, 9, X, X}}, //
1965 {{4, 6, 10, 32, X, X}}, //
1966 {{8, 9, 10, 11, 12, X}}, // upside down pyramid
1967 {{8, 19, 23, 27, 31, X}} // top pyramid
1968 }};
1969 return cell_quads_pyramid;
1970 }
1971
1973 {
1974 // Considerations that went into creating this table
1975 // - Keep aspect ratio of 1.0 if parent had 1.0 (this is why the
1976 // centre wedge is rotated by 180° compared to the parent)
1977 // - Order of children is similar to the order of a reference
1978 // triangle in the x-y-plane (i.e. like the top triangle of the
1979 // parent)
1980 constexpr ::ndarray<unsigned int, 10, 6> cell_quads_wedge =
1981 {{
1982 {{23, 0, 10, 5, 19, X}}, // bottom children
1983 {{22, 1, 11, 14, 4, X}}, //
1984 {{24, 2, 6, 15, 18, X}}, //
1985 {{25, 3, 6, 5, 4, X}}, //
1986 {{0, 26, 12, 8, 21, X}}, // top children
1987 {{1, 27, 13, 16, 7, X}}, //
1988 {{2, 28, 9, 17, 20, X}}, //
1989 {{3, 29, 9, 8, 7, X}}, //
1990 {{X, X, X, X, X, X}},
1991 {{X, X, X, X, X, X}},
1992 }};
1993 return cell_quads_wedge;
1994 }
1995
1997 {
1998 constexpr ::ndarray<unsigned int, 10, 6> cell_quads_hex = {{
1999 {{12, 0, 20, 4, 28, 8}}, // bottom children
2000 {{0, 16, 22, 6, 29, 9}}, //
2001 {{13, 1, 4, 24, 30, 10}}, //
2002 {{1, 17, 6, 26, 31, 11}}, //
2003 {{14, 2, 21, 5, 8, 32}}, // top children
2004 {{2, 18, 23, 7, 9, 33}}, //
2005 {{15, 3, 5, 25, 10, 34}}, //
2006 {{3, 19, 7, 27, 11, 35}}, //
2007 {{X, X, X, X, X, X}},
2008 {{X, X, X, X, X, X}},
2009 }};
2010 return cell_quads_hex;
2011 }
2012 default:
2014 }
2015 }
2016 else
2017 // We should never get here except in the 3d case.
2019
2020 return {};
2021}
2022
2023
2024
2025template <int dim>
2026constexpr ::ndarray<unsigned int, 10, 8>
2028 const unsigned int refinement_choice) const
2029{
2031
2032 if constexpr (dim == 3)
2033 {
2034 constexpr unsigned int X = numbers::invalid_unsigned_int;
2035
2036 switch (this->kind)
2037 {
2039 {
2040 switch (refinement_choice)
2041 {
2042 case 0:
2043 // new line is (6,8)
2044 return {{
2045 {{0, 4, 6, 7, X, X, X, X}},
2046 {{4, 1, 5, 8, X, X, X, X}},
2047 {{6, 5, 2, 9, X, X, X, X}},
2048 {{7, 8, 9, 3, X, X, X, X}},
2049 {{4, 5, 6, 8, X, X, X, X}},
2050 {{4, 7, 8, 6, X, X, X, X}},
2051 {{6, 9, 7, 8, X, X, X, X}},
2052 {{5, 8, 9, 6, X, X, X, X}},
2053 {{X, X, X, X, X, X, X, X}},
2054 {{X, X, X, X, X, X, X, X}},
2055 }};
2056 case 1:
2057 // new line is (5,7)
2058 return {{
2059 {{0, 4, 6, 7, X, X, X, X}},
2060 {{4, 1, 5, 8, X, X, X, X}},
2061 {{6, 5, 2, 9, X, X, X, X}},
2062 {{7, 8, 9, 3, X, X, X, X}},
2063 {{4, 5, 6, 7, X, X, X, X}},
2064 {{4, 7, 8, 5, X, X, X, X}},
2065 {{6, 9, 7, 5, X, X, X, X}},
2066 {{5, 8, 9, 7, X, X, X, X}},
2067 {{X, X, X, X, X, X, X, X}},
2068 {{X, X, X, X, X, X, X, X}},
2069 }};
2070 case 2:
2071 // new line is (4,9)
2072 return {{
2073 {{0, 4, 6, 7, X, X, X, X}},
2074 {{4, 1, 5, 8, X, X, X, X}},
2075 {{6, 5, 2, 9, X, X, X, X}},
2076 {{7, 8, 9, 3, X, X, X, X}},
2077 {{4, 5, 6, 9, X, X, X, X}},
2078 {{4, 7, 8, 9, X, X, X, X}},
2079 {{6, 9, 7, 4, X, X, X, X}},
2080 {{5, 8, 9, 4, X, X, X, X}},
2081 {{X, X, X, X, X, X, X, X}},
2082 {{X, X, X, X, X, X, X, X}},
2083 }};
2084
2085 default:
2087 }
2088 }
2090 {
2091 constexpr ::ndarray<unsigned int, 10, 8>
2092 cell_vertices_pyramid = {{
2093 {{0, 7, 5, 13, 9, X, X, X}}, // bottom pyramid
2094 {{7, 1, 13, 6, 10, X, X, X}}, //
2095 {{5, 13, 2, 8, 11, X, X, X}}, //
2096 {{13, 6, 8, 3, 12, X, X, X}}, //
2097 {{5, 13, 11, 9, X, X, X, X}}, // bottom wedges
2098 {{13, 6, 12, 10, X, X, X, X}}, //
2099 {{7, 13, 9, 10, X, X, X, X}}, //
2100 {{13, 8, 11, 12, X, X, X, X}}, //
2101 {{9, 11, 10, 12, 13, X, X, X}}, // upside down pyramid
2102 {{9, 10, 11, 12, 4, X, X, X}}, // top pyramid
2103 }};
2104 return cell_vertices_pyramid;
2105 }
2106
2108 {
2109 constexpr ::ndarray<unsigned int, 10, 8>
2110 cell_vertices_wedge = {{
2111 {{0, 6, 8, 12, 15, 17, X, X}}, // bottom children
2112 {{6, 1, 7, 15, 13, 16, X, X}}, //
2113 {{8, 7, 2, 17, 16, 14, X, X}}, //
2114 {{7, 8, 6, 16, 17, 15, X, X}}, //
2115 {{12, 15, 17, 3, 9, 11, X, X}}, // top children
2116 {{15, 13, 16, 9, 4, 10, X, X}}, //
2117 {{17, 16, 14, 11, 10, 5, X, X}}, //
2118 {{16, 17, 15, 10, 11, 9, X, X}}, //
2119 {{X, X, X, X, X, X, X, X}},
2120 {{X, X, X, X, X, X, X, X}},
2121 }};
2122 return cell_vertices_wedge;
2123 }
2124
2126 {
2127 constexpr ::ndarray<unsigned int, 10, 8> cell_vertices_hex =
2128 {{
2129 {{0, 10, 8, 24, 16, 22, 20, 26}}, // bottom children
2130 {{10, 1, 24, 9, 22, 17, 26, 21}}, //
2131 {{8, 24, 2, 11, 20, 26, 18, 23}}, //
2132 {{24, 9, 11, 3, 26, 21, 23, 19}}, //
2133 {{16, 22, 20, 26, 4, 14, 12, 25}}, // top children
2134 {{22, 17, 26, 21, 14, 5, 25, 13}}, //
2135 {{20, 26, 18, 23, 12, 25, 6, 15}}, //
2136 {{26, 21, 23, 19, 25, 13, 15, 7}}, //
2137 {{X, X, X, X, X, X, X, X}},
2138 {{X, X, X, X, X, X, X, X}},
2139 }};
2140 return cell_vertices_hex;
2141 }
2142
2143 default:
2145 }
2146 }
2147 else
2148 // We should never get here except in the 3d case.
2150
2151 return {};
2152}
2153
2154
2155
2156template <int dim>
2157inline constexpr bool
2159{
2160 return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
2163}
2164
2165
2166
2167template <int dim>
2168inline constexpr bool
2170{
2171 return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
2172 *this == ReferenceCells::Triangle ||
2174}
2175
2176
2177
2178template <int dim>
2179inline constexpr unsigned int
2181{
2182 return dim;
2183}
2184
2185
2186
2187template <int dim>
2190{
2191 return Quadrature<dim>(std::vector<Point<dim>>({barycenter()}),
2192 std::vector<double>({volume()}));
2193}
2194
2195
2196
2197template <int dim>
2198inline constexpr unsigned int
2200{
2201 switch (this->kind)
2202 {
2204 return 1;
2206 return 2;
2208 return 3;
2210 return 4;
2212 return 4;
2214 return 5;
2216 return 6;
2218 return 8;
2219 default:
2221 }
2222
2224}
2225
2226
2227
2228template <int dim>
2229inline constexpr unsigned int
2231{
2232 switch (this->kind)
2233 {
2235 return 0;
2237 return 1;
2239 return 3;
2241 return 4;
2243 return 6;
2245 return 8;
2247 return 9;
2249 return 12;
2250 default:
2252 }
2253
2255}
2256
2257
2258
2259template <int dim>
2261ReferenceCell<dim>::vertex(const unsigned int v) const
2262{
2264
2265 switch (dim)
2266 {
2267 case 0:
2268 {
2269 if (*this == ReferenceCells::Vertex)
2270 return Point<dim>();
2271 break;
2272 }
2273 case 1:
2274 {
2275 static const Point<dim> vertices[2] = {
2276 Point<dim>(), // the origin
2277 Point<dim>::unit_vector(0) // unit point along x-axis
2278 };
2279 if (*this == ReferenceCells::Line)
2280 return vertices[v];
2281 break;
2282 }
2283 case 2:
2284 {
2285 switch (this->kind)
2286 {
2288 {
2289 static const Point<dim> vertices[3] = {
2290 Point<dim>(), // the origin
2291 Point<dim>::unit_vector(0), // unit point along x-axis
2292 Point<dim>::unit_vector(1) // unit point along y-axis
2293 };
2294 return vertices[v];
2295 }
2297 {
2298 static const Point<dim> vertices[4] = {
2299 // First the two points on the x-axis
2300 Point<dim>(),
2302 // Then these two points shifted in the y-direction
2305 return vertices[v];
2306 }
2307 }
2308 break;
2309 }
2310 case 3:
2311 {
2312 switch (this->kind)
2313 {
2315 {
2316 static const Point<dim> vertices[4] = {
2317 Point<dim>(), // the origin
2318 Point<dim>::unit_vector(0), // unit point along x-axis
2319 Point<dim>::unit_vector(1), // unit point along y-axis
2320 Point<dim>::unit_vector(2) // unit point along z-axis
2321 };
2322 return vertices[v];
2323 }
2325 {
2326 static const Point<dim> vertices[5] = {
2327 Point<dim>{-1.0, -1.0, 0.0},
2328 Point<dim>{+1.0, -1.0, 0.0},
2329 Point<dim>{-1.0, +1.0, 0.0},
2330 Point<dim>{+1.0, +1.0, 0.0},
2331 Point<dim>{+0.0, +0.0, 1.0}};
2332 return vertices[v];
2333 }
2335 {
2336 static const Point<dim> vertices[6] = {
2337 // First the three points on the triangular base of the
2338 // wedge:
2339 Point<dim>(),
2342 // And now everything shifted in the z-direction again
2346 return vertices[v];
2347 }
2349 {
2350 static const Point<dim> vertices[8] = {
2351 // First the two points on the x-axis
2352 Point<dim>(),
2354 // Then these two points shifted in the y-direction
2357 // And now all four points shifted in the z-direction
2364 return vertices[v];
2365 }
2366 }
2367 break;
2368 }
2369 default:
2371 }
2372
2374 return Point<dim>();
2375}
2376
2377
2378
2379template <int dim>
2380inline constexpr unsigned int
2382{
2383 switch (this->kind)
2384 {
2386 return 0;
2388 return 2;
2390 return 3;
2392 return 4;
2394 return 4;
2396 return 5;
2398 return 5;
2400 return 6;
2401 default:
2403 }
2404
2406}
2407
2408
2409
2410template <int dim>
2411inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
2413{
2414 return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(0U,
2415 n_faces());
2416}
2417
2418
2419
2420template <int dim>
2421inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
2423 const ReferenceCell<dim - 1> &face_ref_type) const
2424{
2425 switch (this->kind)
2426 {
2428 return this->face_indices(); // no faces
2430 if (face_ref_type == ReferenceCells::Vertex)
2431 return this->face_indices();
2432 else
2433 return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(0U,
2434 0U);
2437 if (face_ref_type == ReferenceCells::Line)
2438 return this->face_indices();
2439 else
2440 return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(0U,
2441 0U);
2443 if (face_ref_type == ReferenceCells::Triangle)
2444 return this->face_indices();
2445 else
2446 return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(0U,
2447 0U);
2449 switch (face_ref_type)
2450 {
2452 return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(
2453 1U, 5U);
2455 return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(
2456 0U, 1U);
2457 default:
2458 return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(
2459 0U, 0U);
2460 }
2462 switch (face_ref_type)
2463 {
2465 return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(
2466 0U, 2U);
2468 return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(
2469 2U, 5U);
2470 default:
2471 return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(
2472 0U, 0U);
2473 }
2475 if (face_ref_type == ReferenceCells::Quadrilateral)
2476 return this->face_indices();
2477 else
2478 return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(0U,
2479 0U);
2480 default:
2482 }
2483
2484 return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(0U, 0U);
2485}
2486
2487
2488
2489template <int dim>
2490inline constexpr unsigned int
2492{
2493 switch (this->kind)
2494 {
2496 return 1;
2498 return 1;
2500 return 1;
2502 return 1;
2504 return 3;
2506 return 1;
2508 return 1;
2510 return 1;
2511 default:
2513 }
2514
2516}
2517
2518
2519
2520template <int dim>
2523 const unsigned int ref_choice) const
2524{
2526 switch (this->kind)
2527 {
2529 {
2531 isotropic_ref_choices = {{IsotropicRefinementChoice::cut_tet_68,
2534 return isotropic_ref_choices[ref_choice];
2535 }
2536 default:
2537 {
2539 }
2540 }
2541
2543}
2544
2545
2546
2547template <int dim>
2548inline constexpr unsigned int
2550{
2551 switch (this->kind)
2552 {
2554 return 1;
2556 return 2;
2558 return 4;
2560 return 4;
2562 return 8;
2564 return 10;
2566 return 8;
2568 return 8;
2569 default:
2571 }
2572
2574}
2575
2576
2577
2578template <int dim>
2579unsigned int
2581{
2582 // Use GeometryInfo here to keep it the single source of truth
2583 if (this->is_hyper_cube())
2584 return GeometryInfo<dim>::n_children(ref_case);
2585
2586 // TODO: Rewrite once refinement cases are cleaned up
2587 // For all other cells, we assume that only isotropic refinement is
2588 // implemented. This is a strong assumption that should be removed as soon as
2589 // the issue regarding the uniqueness of refinement cases is resolved.
2590 //
2591 // This is a short outline of the issue: Different refinement options are
2592 // stored as an enum, primarily focusing on hypercubes. They therefore are (to
2593 // some extent) unsuitable for other cell types. For example, there are
2594 // multiple ways to isotropically refine a tetrahedron. Distinguishing between
2595 // these required the introduction of the additional
2596 // `IsotropicRefinementChoice` enum. However, this leads to a conflict where
2597 // the same integer value is assigned to multiple different refinement options
2598 // across the existing refinement cases. The most critical conflict is likely
2599 // `IsotropicRefinementChoice::isotropic_refinement` shadowing
2600 // `RefinementPossibilities<~>::no_refinement`.
2601 // The following code attempts to circumvent this mess by predicting whether
2602 // the given refinement case is isotropic and then returning either the number
2603 // of isotropic children or `0` (see below for why `0` is used).
2604
2605 // Detect isotropic refinement.
2606 // Only in 3D and only for Tetrahedrons we have ambiguity. But first assume
2607 // that we don't have a tet.
2608 bool is_isotropic =
2610
2611 // check if assumption of line above was correct, else apply fix for tet
2612 if constexpr (dim == 3) // check for dim first (only 3D relevant for tet)
2613 if (*this == ReferenceCells::Tetrahedron) // tets
2614 is_isotropic =
2615 ((ref_case == //
2617 // We must not check
2618 // `IsotropicRefinementChoice::isotropic_refinement`
2619 // because it coincides with
2620 // `RefinementPossibilities<dim>::no_refinement`
2621 // (both are 0).
2622 (ref_case ==
2623 static_cast<typename RefinementPossibilities<dim>::Possibilities>(
2625 (ref_case ==
2626 static_cast<typename RefinementPossibilities<dim>::Possibilities>(
2628 (ref_case ==
2629 static_cast<typename RefinementPossibilities<dim>::Possibilities>(
2631
2632
2633 // return the according value
2634 if (is_isotropic)
2635 return n_isotropic_children();
2636
2637 // Unsupported case
2638
2639 // DEAL_II_NOT_IMPLEMENTED();
2640 // We must return some data. Many parts of the code aren't aware
2641 // which refinement cases are actually supported, calling this
2642 // function with impossible refinement cases and expecting some
2643 // output. One example is the constructor of `FiniteElement`.
2644 return 0;
2645}
2646
2647
2648
2649template <int dim>
2650inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
2652{
2653 return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(
2654 0U, n_isotropic_children());
2655}
2656
2657
2658
2659template <int dim>
2660inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
2662{
2663 return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(0U,
2664 n_vertices());
2665}
2666
2667
2668
2669template <int dim>
2670inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
2672{
2673 return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(0U,
2674 n_lines());
2675}
2676
2677
2678
2679template <int dim>
2682{
2684
2685 if constexpr (dim == 0)
2686 {
2687 // Like equivalent_refinement_case(), to better enable generic
2688 // programming, return a value implying vertices can be refined since 1d
2689 // Triangulations can be refined.
2690 static constexpr std::array<RefinementCase<0>, 1> possibilities{
2692 return make_array_view(possibilities);
2693 }
2694 else if constexpr (dim == 1)
2695 {
2696 static constexpr std::array<RefinementCase<1>, 1> possibilities{
2698 return make_array_view(possibilities);
2699 }
2700 else if constexpr (dim == 2)
2701 {
2702 switch (this->kind)
2703 {
2705 {
2706 static constexpr std::array<RefinementCase<2>, 1> possibilities{
2708 return make_array_view(possibilities);
2709 }
2711 {
2712 static constexpr std::array<RefinementCase<2>, 3> possibilities{
2716 return make_array_view(possibilities);
2717 }
2718 default:
2720 }
2721 }
2722 else if constexpr (dim == 3)
2723 {
2724 switch (this->kind)
2725 {
2729 {
2730 static constexpr std::array<RefinementCase<3>, 1> possibilities{
2732 return make_array_view(possibilities);
2733 }
2735 {
2736 static constexpr std::array<RefinementCase<3>, 7> possibilities{
2744 return make_array_view(possibilities);
2745 }
2746 default:
2748 }
2749 }
2750
2753}
2754
2755
2756
2757template <int dim>
2758inline constexpr ReferenceCell<dim - 1>
2759ReferenceCell<dim>::face_reference_cell(const unsigned int face_no) const
2760{
2761 AssertIndexRange(face_no, n_faces());
2762
2763 if constexpr (dim == 0)
2765 else if constexpr (dim == 1)
2767 else if constexpr (dim == 2)
2768 return ReferenceCells::Line;
2769 else if constexpr (dim == 3)
2770 {
2771 switch (this->kind)
2772 {
2776 if (face_no == 0)
2778 else
2781 if (face_no > 1)
2783 else
2787 default:
2789 }
2790 }
2791 else
2793}
2794
2795
2796
2797template <int dim>
2798inline constexpr types::geometric_orientation
2803
2804
2805
2806template <int dim>
2807inline constexpr types::geometric_orientation
2809{
2810 // For a reversed line 'orientation' is false and neither flip nor rotate
2811 // are defined.
2813}
2814
2815
2816
2817template <int dim>
2818inline unsigned int
2820 const unsigned int subface) const
2821{
2822 return child_cell_on_face(face,
2823 subface,
2825}
2826
2827
2828
2829template <int dim>
2830inline unsigned int
2832 const unsigned int face,
2833 const unsigned int subface,
2834 const types::geometric_orientation combined_orientation) const
2835{
2836 AssertIndexRange(face, n_faces());
2838
2839 switch (this->kind)
2840 {
2843 {
2845 break;
2846 }
2848 {
2849 static constexpr ndarray<unsigned int, 3, 2> subcells = {
2850 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
2851
2852 Assert(combined_orientation ==
2854 combined_orientation == numbers::reverse_line_orientation,
2856 return subcells[face][combined_orientation ==
2858 subface :
2859 1 - subface];
2860 }
2862 {
2863 const auto [face_orientation, face_rotation, face_flip] =
2864 internal::split_face_orientation(combined_orientation);
2865
2868 face,
2869 subface,
2870 face_orientation,
2871 face_flip,
2872 face_rotation);
2873 }
2877 {
2879 break;
2880 }
2882 {
2883 const auto [face_orientation, face_rotation, face_flip] =
2884 internal::split_face_orientation(combined_orientation);
2885
2888 face,
2889 subface,
2890 face_orientation,
2891 face_flip,
2892 face_rotation);
2893 }
2894 default:
2896 }
2897
2899}
2900
2901
2902
2903template <int dim>
2904inline std::array<unsigned int, 2>
2906 const unsigned int vertex) const
2907{
2909 // Work around a GCC warning at higher optimization levels by making all of
2910 // these tables the same size
2911 constexpr unsigned int X = numbers::invalid_unsigned_int;
2912
2913 switch (this->kind)
2914 {
2917 break;
2919 {
2920 static constexpr ndarray<unsigned int, 2, 2> table = {
2921 {{{0, 0}}, {{1, 0}}}};
2922
2923 return table[vertex];
2924 }
2926 {
2927 static constexpr ndarray<unsigned int, 6, 2> table = {
2928 {{{0, 0}}, {{0, 1}}, {{1, 1}}, {{X, X}}, {{X, X}}, {{X, X}}}};
2929
2930 return table[vertex];
2931 }
2933 {
2935 vertex);
2936 }
2938 {
2939 static constexpr ndarray<unsigned int, 6, 2> table = {
2940 {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 2}}, {{X, X}}, {{X, X}}}};
2941
2942 return table[vertex];
2943 }
2945 {
2946 static constexpr ndarray<unsigned int, 6, 2> table = {
2947 {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{0, 3}}, {{1, 2}}, {{X, X}}}};
2948
2949 return table[vertex];
2950 }
2952 {
2953 static constexpr ndarray<unsigned int, 6, 2> table = {
2954 {{{0, 1}}, {{0, 0}}, {{0, 2}}, {{1, 0}}, {{1, 1}}, {{1, 2}}}};
2955
2956 return table[vertex];
2957 }
2959 {
2961 vertex);
2962 }
2963 default:
2965 }
2966
2967 return {};
2968}
2969
2970
2971
2972template <int dim>
2973inline std::array<unsigned int, 2>
2975 const unsigned int line) const
2976{
2977 AssertIndexRange(line, n_lines());
2978
2979 switch (this->kind)
2980 {
2985 {
2987 break;
2988 }
2990 {
2991 static const std::array<unsigned int, 2> table[6] = {
2992 {{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}};
2993
2994 return table[line];
2995 }
2997 {
2998 static const std::array<unsigned int, 2> table[8] = {{{0, 0}},
2999 {{0, 1}},
3000 {{0, 2}},
3001 {{0, 3}},
3002 {{1, 2}},
3003 {{2, 1}},
3004 {{1, 1}},
3005 {{2, 2}}};
3006
3007 return table[line];
3008 }
3010 {
3011 static const std::array<unsigned int, 2> table[9] = {{{0, 0}},
3012 {{0, 2}},
3013 {{0, 1}},
3014 {{1, 0}},
3015 {{1, 1}},
3016 {{1, 2}},
3017 {{2, 0}},
3018 {{2, 1}},
3019 {{3, 1}}};
3020
3021 return table[line];
3022 }
3024 {
3026 }
3027 default:
3029 }
3030
3031 return {};
3032}
3033
3034
3035
3036template <int dim>
3037inline unsigned int
3039 const unsigned int vertex) const
3040{
3042 AssertIndexRange(line, n_lines());
3043
3044 switch (this->kind)
3045 {
3048 return vertex;
3050 {
3051 static constexpr ndarray<unsigned int, 3, 2> table = {
3052 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
3053 return table[line][vertex];
3054 }
3056 {
3057 static constexpr ndarray<unsigned int, 4, 2> table = {
3058 {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
3059 return table[line][vertex];
3060 }
3062 {
3063 static constexpr ndarray<unsigned int, 6, 2> table = {
3064 {{{0, 1}}, {{1, 2}}, {{2, 0}}, {{0, 3}}, {{1, 3}}, {{2, 3}}}};
3065 return table[line][vertex];
3066 }
3068 {
3069 static constexpr ndarray<unsigned int, 8, 2> table = {{{{0, 2}},
3070 {{1, 3}},
3071 {{0, 1}},
3072 {{2, 3}},
3073 {{0, 4}},
3074 {{1, 4}},
3075 {{2, 4}},
3076 {{3, 4}}}};
3077 return table[line][vertex];
3078 }
3080 {
3081 static constexpr ndarray<unsigned int, 9, 2> table = {{{{0, 1}},
3082 {{1, 2}},
3083 {{2, 0}},
3084 {{3, 4}},
3085 {{4, 5}},
3086 {{5, 3}},
3087 {{0, 3}},
3088 {{1, 4}},
3089 {{2, 5}}}};
3090 return table[line][vertex];
3091 }
3093 {
3094 // first four lines comprise the bottom face, next four are the top,
3095 // and the last four are 'bottom to top'
3096 static constexpr ndarray<unsigned int, 12, 2> table = {{{{0, 2}},
3097 {{1, 3}},
3098 {{0, 1}},
3099 {{2, 3}},
3100 {{4, 6}},
3101 {{5, 7}},
3102 {{4, 5}},
3103 {{6, 7}},
3104 {{0, 4}},
3105 {{1, 5}},
3106 {{2, 6}},
3107 {{3, 7}}}};
3108 return table[line][vertex];
3109 }
3110
3111 default:
3113 }
3114
3116}
3117
3118
3119
3120template <int dim>
3121inline unsigned int
3123 const unsigned int face,
3124 const unsigned int line,
3125 const types::geometric_orientation combined_face_orientation) const
3126{
3127 AssertIndexRange(face, n_faces());
3129
3130 static constexpr unsigned int X = numbers::invalid_unsigned_int;
3131
3132 switch (this->kind)
3133 {
3135 {
3136 // We can't get here based on the assertions above: vertices
3137 // have no faces.
3139 break;
3140 }
3142 {
3143 // We can't get here based on the assertions above: lines
3144 // have vertices as faces, which have no bounding lines.
3146 break;
3147 }
3150 {
3151 // 2d cells have only one line per face, so the line we are
3152 // looking for is actually the face we were given:
3153 return face;
3154 }
3156 {
3157 static constexpr ndarray<unsigned int, 4, 3> table = {
3158 {{{0, 1, 2}}, {{0, 3, 4}}, {{2, 5, 3}}, {{1, 4, 5}}}};
3159
3160 return table[face][standard_to_real_face_line(
3161 line, face, combined_face_orientation)];
3162 }
3164 {
3165 static constexpr ndarray<unsigned int, 5, 4> table = {
3166 {{{0, 1, 2, 3}},
3167 {{0, 6, 4, X}},
3168 {{1, 5, 7, X}},
3169 {{2, 4, 5, X}},
3170 {{3, 7, 6, X}}}};
3171
3172 return table[face][standard_to_real_face_line(
3173 line, face, combined_face_orientation)];
3174 }
3176 {
3177 static constexpr ndarray<unsigned int, 5, 4> table = {
3178 {{{0, 2, 1, X}},
3179 {{3, 4, 5, X}},
3180 {{6, 7, 0, 3}},
3181 {{7, 8, 1, 4}},
3182 {{8, 6, 2, 5}}}};
3183
3184 return table[face][standard_to_real_face_line(
3185 line, face, combined_face_orientation)];
3186 }
3188 {
3189 const auto [face_orientation, face_rotation, face_flip] =
3190 internal::split_face_orientation(combined_face_orientation);
3191
3193 face, line, face_orientation, face_flip, face_rotation);
3194 }
3195 default:
3197 }
3198
3200}
3201
3202
3203
3204template <int dim>
3205inline unsigned int
3207 const unsigned int face,
3208 const unsigned int vertex,
3209 const types::geometric_orientation combined_face_orientation) const
3210{
3211 AssertIndexRange(face, n_faces());
3213 AssertIndexRange(combined_face_orientation, n_face_orientations(face));
3214
3215 switch (this->kind)
3216 {
3218 {
3219 // We can't get here based on the assertions above: vertices
3220 // have no faces.
3222 break;
3223 }
3225 {
3226 // A line has only one vertex per face, so the vertex we are
3227 // looking for is actually the face we were given:
3228 return face;
3229 }
3231 {
3232 static constexpr ndarray<unsigned int, 3, 2> table = {
3233 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
3234
3235 return table[face][combined_face_orientation ==
3237 vertex :
3238 (1 - vertex)];
3239 }
3241 {
3242 const auto [face_orientation, face_rotation, face_flip] =
3243 internal::split_face_orientation(combined_face_orientation);
3244
3246 face, vertex, face_orientation, face_flip, face_rotation);
3247 }
3249 {
3250 static constexpr ndarray<unsigned int, 4, 3> table = {
3251 {{{0, 1, 2}}, {{1, 0, 3}}, {{0, 2, 3}}, {{2, 1, 3}}}};
3252
3253 return table[face][standard_to_real_face_vertex(
3254 vertex, face, combined_face_orientation)];
3255 }
3257 {
3258 constexpr auto X = numbers::invalid_unsigned_int;
3259 static constexpr ndarray<unsigned int, 5, 4> table = {
3260 {{{0, 1, 2, 3}},
3261 {{0, 2, 4, X}},
3262 {{3, 1, 4, X}},
3263 {{1, 0, 4, X}},
3264 {{2, 3, 4, X}}}};
3265
3266 return table[face][standard_to_real_face_vertex(
3267 vertex, face, combined_face_orientation)];
3268 }
3270 {
3271 constexpr auto X = numbers::invalid_unsigned_int;
3272 static constexpr ndarray<unsigned int, 6, 4> table = {
3273 {{{1, 0, 2, X}},
3274 {{3, 4, 5, X}},
3275 {{0, 1, 3, 4}},
3276 {{1, 2, 4, 5}},
3277 {{2, 0, 5, 3}}}};
3278
3279 return table[face][standard_to_real_face_vertex(
3280 vertex, face, combined_face_orientation)];
3281 }
3283 {
3284 const auto [face_orientation, face_rotation, face_flip] =
3285 internal::split_face_orientation(combined_face_orientation);
3286
3288 face, vertex, face_orientation, face_flip, face_rotation);
3289 }
3290 default:
3292 }
3293
3295}
3296
3297
3298
3299template <int dim>
3300inline unsigned int
3302 const unsigned int face,
3303 const unsigned int line,
3304 const unsigned int vertex,
3305 const types::geometric_orientation combined_face_orientation) const
3306{
3307 AssertIndexRange(face, n_faces());
3308 AssertIndexRange(line, this->face_reference_cell(face).n_lines());
3310 AssertIndexRange(combined_face_orientation, n_face_orientations(face));
3311
3312 switch (this->kind)
3313 {
3315 {
3316 // We can't get here based on the assertions above: vertices
3317 // have no faces.
3319 break;
3320 }
3322 {
3323 // We can't get here based on the assertions above: Lines'
3324 // faces are vertices that have no lines as sub-objects.
3326 break;
3327 }
3330 {
3331 // For 2d objects, each face has (and is identical to) exactly
3332 // one line, so we can ask the other function:
3333 return face_to_cell_vertices(face, vertex, combined_face_orientation);
3334 }
3336 {
3337 Assert(combined_face_orientation ==
3340
3341 const static ::ndarray<unsigned int, 4, 3, 2> table = {
3342 {{{{{0, 1}}, {{1, 2}}, {{2, 0}}}},
3343 {{{{1, 0}}, {{0, 3}}, {{3, 1}}}},
3344 {{{{0, 2}}, {{2, 3}}, {{3, 0}}}},
3345 {{{{2, 1}}, {{1, 3}}, {{3, 2}}}}}};
3346
3347 return table[face][line][vertex];
3348 }
3350 {
3351 Assert(combined_face_orientation ==
3354
3355 static const unsigned int X = static_cast<unsigned int>(-1);
3356 const static ::ndarray<unsigned int, 5, 4, 2> table = {
3357 {{{{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}},
3358 {{{{0, 2}}, {{2, 4}}, {{4, 0}}, {{X, X}}}},
3359 {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}},
3360 {{{{1, 0}}, {{0, 4}}, {{4, 1}}, {{X, X}}}},
3361 {{{{2, 3}}, {{3, 4}}, {{4, 2}}, {{X, X}}}}}};
3362
3363 return table[face][line][vertex];
3364 }
3366 {
3367 Assert(combined_face_orientation ==
3370
3371 static const unsigned int X = static_cast<unsigned int>(-1);
3372 const static ::ndarray<unsigned int, 5, 4, 2> table = {
3373 {{{{{1, 0}}, {{0, 2}}, {{2, 1}}, {{X, X}}}},
3374 {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}},
3375 {{{{0, 3}}, {{1, 4}}, {{0, 1}}, {{3, 4}}}},
3376 {{{{1, 4}}, {{2, 5}}, {{1, 2}}, {{4, 5}}}},
3377 {{{{2, 5}}, {{0, 3}}, {{2, 0}}, {{5, 3}}}}}};
3378
3379 return table[face][line][vertex];
3380 }
3382 {
3383 Assert(combined_face_orientation ==
3386
3387 const static ::ndarray<unsigned int, 6, 4, 2> table = {
3388 {{{{{0, 4}}, {{2, 6}}, {{0, 2}}, {{4, 6}}}},
3389 {{{{1, 5}}, {{3, 7}}, {{1, 3}}, {{5, 7}}}},
3390 {{{{0, 1}}, {{4, 5}}, {{0, 4}}, {{1, 5}}}},
3391 {{{{2, 3}}, {{6, 7}}, {{2, 6}}, {{3, 7}}}},
3392 {{{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}},
3393 {{{{4, 6}}, {{5, 7}}, {{4, 5}}, {{6, 7}}}}}};
3394
3395 return table[face][line][vertex];
3396 }
3397 default:
3399 }
3400
3402}
3403
3404
3405
3406template <int dim>
3409 const unsigned int vertex) const
3410{
3411 return this->vertex(face_to_cell_vertices(
3413}
3414
3415
3416
3417template <int dim>
3418inline unsigned int
3420 const unsigned int vertex,
3421 const unsigned int face,
3422 const types::geometric_orientation face_orientation) const
3423{
3424 AssertIndexRange(face, n_faces());
3426 AssertIndexRange(face_orientation, n_face_orientations(face));
3427
3428 switch (face_reference_cell(face))
3429 {
3431 // test to ensure that face_orientation is
3432 // default_geometric_orientation already done with
3433 // AssertIndexRange(face_orientation, ...) above.
3434 return vertex;
3436 return line_vertex_permutations[face_orientation][vertex];
3438 return triangle_vertex_permutations[face_orientation][vertex];
3440 return quadrilateral_vertex_permutations[face_orientation][vertex];
3441 default:
3443 }
3444
3447}
3448
3449
3450
3451template <int dim>
3452inline unsigned int
3454 const unsigned int line,
3455 const unsigned int face,
3456 const types::geometric_orientation face_orientation) const
3457{
3458 AssertIndexRange(face, n_faces());
3460 AssertIndexRange(face_orientation, n_face_orientations(face));
3461
3462 switch (face_reference_cell(face))
3463 {
3465 return triangle_line_permutations[face_orientation][line];
3467 return quadrilateral_line_permutations[face_orientation][line];
3468 // case ReferenceCells::Vertex:
3469 // case ReferenceCells::Line:
3470 // case ReferenceCells::Invalid:
3471 default:
3473 }
3474
3477}
3478
3479
3480
3481namespace ReferenceCells
3482{
3483 template <int dim>
3484 constexpr const std::vector<ReferenceCell<dim>>
3486 {
3487 if constexpr (dim == 0)
3488 return {{ReferenceCells::Vertex}};
3489 else if constexpr (dim == 1)
3490 return {{ReferenceCells::Line}};
3491 else if constexpr (dim == 2)
3493 else if constexpr (dim == 3)
3498 else
3500
3501 return {};
3502 }
3503
3504
3505
3506 template <int dim>
3507 inline constexpr const ReferenceCell<dim> &
3509 {
3510 if constexpr (dim == 0)
3512 else if constexpr (dim == 1)
3513 return ReferenceCells::Line;
3514 else if constexpr (dim == 2)
3516 else if constexpr (dim == 3)
3518 else
3520
3522 }
3523
3524
3525
3526 template <int dim>
3527 inline constexpr const ReferenceCell<dim> &
3529 {
3530 if constexpr (dim == 0)
3532 else if constexpr (dim == 1)
3533 return ReferenceCells::Line;
3534 else if constexpr (dim == 2)
3536 else if constexpr (dim == 3)
3538 else
3541 }
3542
3543
3544
3545 template <int dim>
3546 inline constexpr unsigned int
3548 {
3549 if constexpr (dim == 0 || dim == 1)
3550 return 1; // Vertex or Line
3551 else if (dim == 2)
3552 return 2; // Quad and Tri
3553 else if (dim == 3)
3554 return 4; // Tet, Pyramid, Wedge and Hex
3555 else
3557 return 0;
3558 }
3559
3560
3561
3562 template <int dim>
3563 inline ReferenceCell<dim>
3564 n_vertices_to_reference_cell(const unsigned int n_vertices)
3565 {
3566 if constexpr (dim == 0)
3567 if (n_vertices == 1)
3569
3570 if constexpr (dim == 1)
3571 {
3572 if (n_vertices == 2)
3573 return ReferenceCells::Line;
3574 }
3575
3576 if constexpr (dim == 2)
3577 {
3578 if (n_vertices == 3)
3580 else if (n_vertices == 4)
3582 }
3583
3584 if constexpr (dim == 3)
3585 {
3586 if (n_vertices == 4)
3588 else if (n_vertices == 5)
3590 if (n_vertices == 6)
3591 return ReferenceCells::Wedge;
3592 else if (n_vertices == 8)
3594 }
3595
3596 Assert(false,
3597 ExcMessage("The combination of dim = " + std::to_string(dim) +
3598 " and n_vertices = " + std::to_string(n_vertices) +
3599 " does not correspond to a known reference cell type."));
3601 }
3602
3603
3604
3605 template <int structdim>
3606 inline constexpr unsigned int
3611
3612
3613
3614 template <int structdim>
3615 inline constexpr unsigned int
3620
3621
3622
3623 template <int structdim>
3624 inline constexpr unsigned int
3629
3630
3631
3632 template <int structdim>
3633 inline constexpr unsigned int
3635 {
3636 // Pyramids (3D) have 10 children.
3637 return structdim == 3 ? 10 : GeometryInfo<structdim>::max_children_per_cell;
3638 }
3639} // namespace ReferenceCells
3640
3641
3642
3643template <int dim>
3644inline ReferenceCell<dim>
3646 const unsigned int n_vertices)
3647{
3648 AssertDimension(dim, dim_);
3650}
3651
3652
3653
3654template <int dim>
3655inline double
3657 const unsigned int i) const
3658{
3660
3661 switch (this->kind)
3662 {
3664 return 1.0;
3668 if constexpr (dim > 0)
3671 // see also BarycentricPolynomials<2>::compute_value
3673 {
3674 switch (i)
3675 {
3676 case 0:
3677 return 1.0 - xi[std::min(0, dim - 1)] -
3678 xi[std::min(1, dim - 1)];
3679 case 1:
3680 return xi[std::min(0, dim - 1)];
3681 case 2:
3682 return xi[std::min(1, dim - 1)];
3683 default:
3685 }
3686 }
3687 // see also BarycentricPolynomials<3>::compute_value
3689 {
3690 switch (i)
3691 {
3692 case 0:
3693 return 1.0 - xi[std::min(0, dim - 1)] -
3694 xi[std::min(1, dim - 1)] - xi[std::min(2, dim - 1)];
3695 case 1:
3696 return xi[std::min(0, dim - 1)];
3697 case 2:
3698 return xi[std::min(1, dim - 1)];
3699 case 3:
3700 return xi[std::min(2, dim - 1)];
3701 default:
3703 }
3704 }
3705 // see also ScalarLagrangePolynomialPyramid::compute_value()
3707 {
3708 const double Q14 = 0.25;
3709
3710 const double r = xi[std::min(0, dim - 1)];
3711 const double s = xi[std::min(1, dim - 1)];
3712 const double t = xi[std::min(2, dim - 1)];
3713
3714 const double ratio =
3715 (std::fabs(t - 1.0) > 1.0e-14 ? (r * s * t) / (1.0 - t) : 0.0);
3716
3717 if (i == 0)
3718 return Q14 * ((1.0 - r) * (1.0 - s) - t + ratio);
3719 if (i == 1)
3720 return Q14 * ((1.0 + r) * (1.0 - s) - t - ratio);
3721 if (i == 2)
3722 return Q14 * ((1.0 - r) * (1.0 + s) - t - ratio);
3723 if (i == 3)
3724 return Q14 * ((1.0 + r) * (1.0 + s) - t + ratio);
3725 else
3726 return t;
3727 }
3728 // see also ScalarLagrangePolynomialWedge::compute_value()
3730 return ReferenceCells::Triangle.d_linear_shape_function(
3731 Point<2>(xi[std::min(0, dim - 1)], xi[std::min(1, dim - 1)]),
3732 i % 3) *
3733 ReferenceCells::Line.d_linear_shape_function(
3734 Point<1>(xi[std::min(2, dim - 1)]), i / 3);
3735 default:
3737 }
3738
3740 return 0.0;
3741}
3742
3743
3744
3745template <int dim>
3746inline Tensor<1, dim>
3748 const unsigned int i) const
3749{
3750 switch (this->kind)
3751 {
3757 // see also BarycentricPolynomials<2>::compute_grad()
3759 switch (i)
3760 {
3761 case 0:
3762 return Point<dim>(-1.0, -1.0);
3763 case 1:
3764 return Point<dim>(+1.0, +0.0);
3765 case 2:
3766 return Point<dim>(+0.0, +1.0);
3767 default:
3769 }
3770 default:
3772 }
3773
3774 return Point<dim>(+0.0, +0.0, +0.0);
3775}
3776
3777
3778
3779template <int dim>
3780inline double
3782{
3783 switch (this->kind)
3784 {
3786 return 1;
3788 return 1;
3790 return 1. / 2.;
3792 return 1;
3794 return 1. / 6.;
3796 return 4. / 3.;
3798 return 1. / 2.;
3800 return 1;
3801 default:
3803 }
3804
3805 return 0.0;
3806}
3807
3808
3809
3810template <int dim>
3811inline double
3812ReferenceCell<dim>::face_measure(const unsigned int face_no) const
3813{
3814 AssertIndexRange(face_no, n_faces());
3815 constexpr double X = std::numeric_limits<double>::signaling_NaN();
3816
3817 static const std::array<double, 5> tri_faces{
3818 {1.0, std::sqrt(2.0), 1.0, X, X}};
3819 static const std::array<double, 5> tet_faces{
3820 {0.5, 0.5, 0.5, std::sqrt(3.0) / 2.0, X}};
3821 static const std::array<double, 5> pyramid_faces{
3822 {4, std::sqrt(2.0), std::sqrt(2.0), std::sqrt(2.0), std::sqrt(2.0)}};
3823 static const std::array<double, 5> wedge_faces{
3824 {0.5, 0.5, 1.0, std::sqrt(2.0), 1.0}};
3825
3826 switch (this->kind)
3827 {
3832 return 1.0;
3834 return tri_faces[face_no];
3836 return tet_faces[face_no];
3838 return pyramid_faces[face_no];
3840 return wedge_faces[face_no];
3841 }
3843 return X;
3844}
3845
3846
3847
3848template <int dim>
3849inline Point<dim>
3851{
3852 switch (this->kind)
3853 {
3855 return Point<dim>();
3857 return Point<dim>(1. / 2.);
3859 return Point<dim>(1. / 3., 1. / 3.);
3861 return Point<dim>(1. / 2., 1. / 2.);
3863 return Point<dim>(1. / 4., 1. / 4., 1. / 4.);
3865 return Point<dim>(0, 0, 1. / 4.);
3867 return Point<dim>(1. / 3, 1. / 3, 1. / 2.);
3869 return Point<dim>(1. / 2., 1. / 2., 1. / 2.);
3870 default:
3872 }
3873
3874 return Point<dim>();
3875}
3876
3877
3878
3879template <int dim>
3880inline bool
3882 const double tolerance) const
3883{
3884 // Introduce abbreviations to silence compiler warnings about invalid
3885 // array accesses (that can't happen, of course, but the compiler
3886 // doesn't know that).
3887 constexpr unsigned int x_coordinate = 0;
3888 constexpr unsigned int y_coordinate = (dim >= 2 ? 1 : 0);
3889 constexpr unsigned int z_coordinate = (dim >= 3 ? 2 : 0);
3890
3891 switch (this->kind)
3892 {
3894 {
3895 // Vertices are special cases in that they do not actually
3896 // have coordinates. Error out if this function is called
3897 // with a vertex:
3898 Assert(false,
3899 ExcMessage("Vertices are zero-dimensional objects and "
3900 "as a consequence have no coordinates. You "
3901 "cannot meaningfully ask whether a point is "
3902 "inside a vertex (within a certain tolerance) "
3903 "without coordinate values."));
3904 return false;
3905 }
3909 {
3910 for (unsigned int d = 0; d < dim; ++d)
3911 if ((p[d] < -tolerance) || (p[d] > 1 + tolerance))
3912 return false;
3913 return true;
3914 }
3917 {
3918 // First make sure that we are in the first quadrant or octant
3919 for (unsigned int d = 0; d < dim; ++d)
3920 if (p[d] < -tolerance)
3921 return false;
3922
3923 // Now we also need to make sure that we are below the diagonal line
3924 // or plane that delineates the simplex. This diagonal is given by
3925 // sum(p[d])<=1, and a diagonal a distance eps away is given by
3926 // sum(p[d])<=1+eps*sqrt(d). (For example, the point at (1,1) is a
3927 // distance of 1/sqrt(2) away from the diagonal. That is, its
3928 // sum satisfies
3929 // sum(p[d]) = 2 <= 1 + (1/sqrt(2)) * sqrt(2)
3930 // in other words, it satisfies the predicate with eps=1/sqrt(2).)
3931 double sum = 0;
3932 for (unsigned int d = 0; d < dim; ++d)
3933 sum += p[d];
3934 return (sum <= 1 + tolerance * std::sqrt(1. * dim));
3935 }
3937 {
3938 // A pyramid only lives in the upper half-space:
3939 if (p[z_coordinate] < -tolerance)
3940 return false;
3941
3942 // It also only lives in the space below z=1:
3943 if (p[z_coordinate] > 1 + tolerance)
3944 return false;
3945
3946 // Within what's left of the space, a pyramid is a cone that tapers
3947 // towards the top. First compute the distance of the point to the
3948 // axis in the max norm (this is the right norm because the vertices
3949 // of the pyramid are at points +/-1, +/-1):
3950 const double distance_from_axis =
3951 std::max(std::fabs(p[x_coordinate]), std::fabs(p[y_coordinate]));
3952
3953 // We are inside the pyramid if the distance from the axis is less
3954 // than (1-z)
3955 return (distance_from_axis <= 1 + tolerance - p[z_coordinate]);
3956 }
3958 {
3959 // The wedge we use is a triangle extruded into the third
3960 // dimension by one unit. So we can use the same logic as for
3961 // triangles above (i.e., for the simplex above, using dim==2)
3962 // and then check the third dimension separately.
3963
3964 if ((p[x_coordinate] < -tolerance) || (p[y_coordinate] < -tolerance))
3965 return false;
3966
3967 const double sum = p[x_coordinate] + p[y_coordinate];
3968 if (sum > 1 + tolerance * std::sqrt(2.0))
3969 return false;
3970
3971 if (p[z_coordinate] < -tolerance)
3972 return false;
3973 if (p[z_coordinate] > 1 + tolerance)
3974 return false;
3975
3976 return true;
3977 }
3978 default:
3980 }
3981
3982 return false;
3983}
3984
3985
3986
3987template <int dim>
3988inline Tensor<1, dim>
3990 const unsigned int i) const
3991{
3992 return face_tangent_vector(face_no, i);
3993}
3994
3995
3996
3997template <int dim>
3998inline Tensor<1, dim>
4000 const unsigned int i) const
4001{
4002 AssertIndexRange(face_no, n_faces());
4003
4004 // Simplify the 1d case by setting tangents (which are used for boundary
4005 // forms and thus Jacobians) equal to the unit vector
4006 if constexpr (dim == 1)
4007 {
4008 return Tensor<1, dim>{{1.0}};
4009 }
4010
4011 AssertIndexRange(i, dim - 1);
4012
4013 switch (this->kind)
4014 {
4022 if constexpr (dim == 2)
4023 {
4024 AssertIndexRange(face_no, 3);
4025 constexpr std::array<Tensor<1, dim>, 3> table = {
4026 {Point<dim>(1, 0),
4028 Point<dim>(0, -1)}};
4029
4030 return table[face_no];
4031 }
4032 else
4035 if constexpr (dim == 3)
4036 {
4037 AssertIndexRange(face_no, 4);
4038 // We need std::pow(1.0/3.0, 0.25) in a constexpr context, but that
4039 // function isn't constexpr yet so hard-code the value and assert
4040 // that it is close
4041 constexpr double third_r4 = 0.7598356856515925473311877;
4042 Assert(std::abs((third_r4 * third_r4) * (third_r4 * third_r4) -
4043 1.0 / 3.0) < 1e-14,
4045 constexpr ndarray<Tensor<1, dim>, 4, 2> table = {
4046 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
4047 {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
4048 {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
4049 {{Point<dim>(-third_r4, +third_r4, 0),
4050 Point<dim>(-third_r4, 0, +third_r4)}}}};
4051
4052 return table[face_no][i];
4053 }
4054 else
4057 if constexpr (dim == 3)
4058 {
4059 AssertIndexRange(face_no, 5);
4060 constexpr ndarray<Tensor<1, dim>, 5, 2> table = {
4061 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
4063 Point<dim>(0, 1, 0)}},
4065 Point<dim>(0, 1, 0)}},
4066 {{Point<dim>(1, 0, 0),
4068 {{Point<dim>(1, 0, 0),
4070
4071 return table[face_no][i];
4072 }
4073 else
4076 if constexpr (dim == 3)
4077 {
4078 AssertIndexRange(face_no, 5);
4079 constexpr ndarray<Tensor<1, dim>, 5, 2> table = {
4080 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
4081 {{Point<dim>(1, 0, 0), Point<dim>(0, 1, 0)}},
4082 {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
4084 Point<dim>(0, 0, 1)}},
4085 {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}}}};
4086
4087 return table[face_no][i];
4088 }
4089 else
4091 default:
4093 }
4094
4095
4096 return {};
4097}
4098
4099
4100
4101template <int dim>
4102inline Tensor<1, dim>
4103ReferenceCell<dim>::unit_normal_vectors(const unsigned int face_no) const
4104{
4105 return face_normal_vector(face_no);
4106}
4107
4108
4109
4110template <int dim>
4111inline Tensor<1, dim>
4112ReferenceCell<dim>::face_normal_vector(const unsigned int face_no) const
4113{
4114 if (is_hyper_cube())
4115 {
4118 }
4119 else if (dim == 2)
4120 {
4122
4123 // Return the rotated vector
4124 return cross_product_2d(face_tangent_vector(face_no, 0));
4125 }
4126 else if (dim == 3)
4127 {
4128 return cross_product_3d(face_tangent_vector(face_no, 0),
4129 face_tangent_vector(face_no, 1));
4130 }
4131
4133
4134 return {};
4135}
4136
4137
4138
4139template <int dim>
4140inline unsigned int
4141ReferenceCell<dim>::n_face_orientations(const unsigned int face_no) const
4142{
4143 AssertIndexRange(face_no, n_faces());
4144
4145 switch (face_reference_cell(face_no))
4146 {
4148 return 1;
4150 return 2;
4152 return 6;
4154 return 8;
4155 default:
4157 }
4158
4161}
4162
4163
4164
4165template <int dim>
4168 const unsigned int line,
4169 const unsigned int face,
4170 const types::geometric_orientation combined_face_orientation,
4171 const types::geometric_orientation line_orientation) const
4172{
4174 face,
4175 combined_face_orientation,
4176 line_orientation);
4177}
4178
4179
4180
4181template <int dim>
4184 const unsigned int face_line_no,
4185 const unsigned int face_no,
4186 const types::geometric_orientation combined_face_orientation,
4187 const types::geometric_orientation line_orientation) const
4188{
4189 constexpr auto D = numbers::default_geometric_orientation;
4190 constexpr auto R = numbers::reverse_line_orientation;
4191 static constexpr unsigned int X = numbers::invalid_unsigned_int;
4192
4194 {
4195 static constexpr ::ndarray<types::geometric_orientation, 2, 8>
4196 orientation_quad{
4197 {{{D, D, D, R, R, R, R, D}}, {{D, D, R, D, R, R, D, R}}}};
4198 // We use face_line_no / 2 here since lines i and i + 1 are parallel
4199 // and, on a given face, have the same relative orientations.
4200 const bool match =
4201 line_orientation ==
4202 orientation_quad[face_line_no / 2][combined_face_orientation];
4203
4206 }
4207 else if (this->face_reference_cell(face_no) == ReferenceCells::Triangle)
4208 {
4209 Assert(face_line_no < 3,
4210 ExcMessage("Invalid face_line_no number for a triangular face"));
4211
4212 // if the line in the face and the line on the reference cell have the
4213 // same orientation, then switching vertices 1 and 2 reverses the
4214 // orientation (i.e. in the first array the 2nd, the 4th and the last
4215 // entry) but rotating the triangle does not switch the orientation (the
4216 // remaining entries)
4217 // if the line in the face and the line on the reference cell do not have
4218 // the same orientation, then it is the other way around
4219 static constexpr ::ndarray<types::geometric_orientation, 2, 6>
4220 orientation_triangle{{{{D, R, D, R, D, R}}, {{R, D, R, D, R, D}}}};
4221
4222 unsigned int combined_line;
4223
4224 if (*this == ReferenceCells::Tetrahedron)
4225 {
4226 static constexpr ::ndarray<unsigned int, 4, 3> combined_lines{
4227 {{{0, 0, 0}}, {{X, 0, 1}}, {{X, 0, X}}, {{X, X, X}}}};
4228
4229 combined_line = combined_lines[face_no][face_line_no];
4230
4231 Assert(combined_line != X,
4232 ExcMessage(
4233 "This function can only be called for following face-line "
4234 "combinations: (0,0), (0,1), (0,2), (1,1), (1,2), (2,1)"));
4235 }
4236 else if (*this == ReferenceCells::Pyramid)
4237 {
4238 // face_line_no 0 already done above
4239 Assert(face_line_no == 1 || face_line_no == 2,
4240 ExcMessage(
4241 "Invalid face line number for the triangular faces"));
4242
4243 // face 0 with lines 0-3 already done above
4244 // lines 6 (face 1, face_line_no 1) and 5 (face 2, face_line_no 1) are
4245 // oriented the same way in the cell and on the faces
4246 // lines 4 (face 1, face_line 2) and 7 (face 2, face_line 2) are not
4247 // oriented the same way in the cell and on the faces
4248 static constexpr ::ndarray<unsigned int, 3, 3> combined_lines{
4249 {{{X, X, X}}, {{X, 0, 1}}, {{X, 0, 1}}}};
4250
4251 combined_line = combined_lines[face_no][face_line_no];
4252 }
4253 else if (*this == ReferenceCells::Wedge)
4254 {
4255 // lines 0,1,2 (face 0, face_line_no 0,1,2) are not
4256 // oriented the same way in the cell and on the face
4257 // lines 4,5,6 (face 1, face_line 0,1,2) are
4258 // oriented the same way in the cell and on the faces
4259 static constexpr ::ndarray<unsigned int, 2, 3> combined_lines{
4260 {{{1, 1, 1}}, {{0, 0, 0}}}};
4261
4262 combined_line = combined_lines[face_no][face_line_no];
4263 }
4264 else
4266
4267 const bool match =
4268 line_orientation ==
4269 orientation_triangle[combined_line][combined_face_orientation];
4270
4273 }
4274 else
4277}
4278
4279
4280
4281template <int dim>
4282inline unsigned int
4283ReferenceCell<dim>::opposite_face_index(const unsigned int face_no) const
4284{
4285 AssertIndexRange(face_no, n_faces());
4286 switch (this->kind)
4287 {
4291 {
4292 return GeometryInfo<2>::opposite_face[face_no];
4293 }
4295 {
4296 return GeometryInfo<3>::opposite_face[face_no];
4297 }
4301 {
4302 return (face_no + 1) % (this->get_dimension() + 1);
4303 }
4305 {
4306 switch (face_no)
4307 {
4308 case 0:
4309 return 0;
4310 case 1:
4311 return 2;
4312 case 2:
4313 return 1;
4314 case 3:
4315 return 4;
4316 case 4:
4317 return 3;
4318 default:
4320 }
4321 }
4323 {
4324 switch (face_no)
4325 {
4326 case 0:
4327 return 1;
4328 case 1:
4329 return 0;
4330 case 2:
4331 return 3;
4332 case 3:
4333 return 4;
4334 case 4:
4335 return 2;
4336 default:
4338 }
4339 }
4340 default:
4342 }
4344}
4345
4346namespace internal
4347{
4348 template <int dim, typename T>
4350 {
4351 public:
4365
4369 virtual ~NoPermutation() noexcept override = default;
4370
4374 virtual void
4375 print_info(std::ostream &out) const override
4376 {
4377 out << '[';
4378
4379 const unsigned int n_vertices = entity_type.n_vertices();
4380
4381 for (unsigned int i = 0; i < n_vertices; ++i)
4382 {
4383 out << vertices_0[i];
4384 if (i + 1 != n_vertices)
4385 out << ',';
4386 }
4387
4388 out << "] is not a valid permutation of [";
4389
4390 for (unsigned int i = 0; i < n_vertices; ++i)
4391 {
4392 out << vertices_1[i];
4393 if (i + 1 != n_vertices)
4394 out << ',';
4395 }
4396
4397 out << "]." << std::endl;
4398 }
4399
4404
4409
4414 };
4415
4416} // namespace internal
4417
4418
4419
4420template <int dim>
4421template <typename T, std::size_t N>
4424 const std::array<T, N> &vertices_0,
4425 const std::array<T, N> &vertices_1) const
4426{
4427 Assert(N >= n_vertices(),
4428 ExcMessage("The number of array elements must be equal to or "
4429 "greater than the number of vertices of the cell "
4430 "referenced by this object."));
4431
4432 // Call the non-deprecated function, taking care of calling it only with
4433 // those array elements that we actually care about (see the note
4434 // in the documentation about the arguments potentially being
4435 // larger arrays than necessary).
4437 make_array_view(vertices_0.begin(), vertices_0.begin() + n_vertices()),
4438 make_array_view(vertices_1.begin(), vertices_1.begin() + n_vertices()));
4439}
4440
4441
4442
4443template <int dim>
4444template <typename T>
4447 const ArrayView<const T> &vertices_0,
4448 const ArrayView<const T> &vertices_1) const
4449{
4450 Assert(vertices_0.size() == n_vertices(),
4451 ExcMessage("The number of array elements must be equal to "
4452 "the number of vertices of the cell "
4453 "referenced by this object."));
4454 Assert(vertices_1.size() == n_vertices(),
4455 ExcMessage("The number of array elements must be equal to "
4456 "the number of vertices of the cell "
4457 "referenced by this object."));
4458
4459 auto compute_orientation = [&](const auto &table,
4460 const unsigned int n_vertices) {
4461 for (types::geometric_orientation o = 0; o < table.size(); ++o)
4462 {
4463 bool match = true;
4464 for (unsigned int j = 0; j < n_vertices; ++j)
4465 match = (match && vertices_0[j] == vertices_1[table[o][j]]);
4466
4467 if (match)
4468 return o;
4469 }
4470
4471 // Do not use `this` in Assert because nvcc when using C++20 assumes that
4472 // `this` is an integer and we get the following error: invalid type
4473 // argument of unary '*' (have 'int')
4474 [[maybe_unused]] const auto &ref_cell = *this;
4475 Assert(false,
4476 (internal::NoPermutation<dim, T>(ref_cell, vertices_0, vertices_1)));
4477 return std::numeric_limits<types::geometric_orientation>::max();
4478 };
4479
4480 switch (this->kind)
4481 {
4490 default:
4492 }
4493
4494 Assert(false,
4495 (internal::NoPermutation<dim, T>(*this, vertices_0, vertices_1)));
4496 return std::numeric_limits<types::geometric_orientation>::max();
4497}
4498
4499
4500
4501template <int dim>
4502template <typename T, std::size_t N>
4503inline std::array<T, N>
4505 const std::array<T, N> &vertices,
4506 const unsigned int orientation) const
4507{
4508 Assert(N >= n_vertices(),
4509 ExcMessage("The number of array elements must be equal to or "
4510 "greater than the number of vertices of the cell "
4511 "referenced by this object."));
4512
4513 // Call the non-deprecated function, taking care of calling it only with
4514 // those array elements that we actually care about (see the note
4515 // in the documentation about the arguments potentially being
4516 // larger arrays than necessary).
4517 const auto permutation = permute_by_combined_orientation(
4518 make_array_view(vertices.begin(), vertices.begin() + n_vertices()),
4519 orientation);
4520
4521 std::array<T, N> temp;
4522 std::copy(permutation.begin(), permutation.end(), temp.begin());
4523
4524 return temp;
4525}
4526
4527
4528
4529template <int dim>
4530template <typename T>
4531boost::container::small_vector<T, 8>
4533 const ArrayView<const T> &vertices,
4534 const types::geometric_orientation orientation) const
4535{
4536 AssertDimension(vertices.size(), n_vertices());
4537 boost::container::small_vector<T, 8> result(n_vertices());
4538
4539 auto permute = [&](const auto &table, unsigned int n_vertices) {
4540 AssertIndexRange(orientation, table.size());
4541 for (unsigned int j = 0; j < n_vertices; ++j)
4542 result[j] = vertices[table[orientation][j]];
4543 };
4544
4545 switch (this->kind)
4546 {
4548 permute(vertex_vertex_permutations, 1);
4549 break;
4551 permute(line_vertex_permutations, 2);
4552 break;
4554 permute(triangle_vertex_permutations, 3);
4555 break;
4558 break;
4559 default:
4561 }
4562
4563 return result;
4564}
4565
4566
4567
4568template <int dim>
4571 const types::geometric_orientation orientation) const
4572{
4573 switch (this->kind)
4574 {
4576 // Things are always default-oriented in 1D
4577 return orientation;
4578
4580 // the 1d orientations are the identity and a flip: i.e., the identity
4581 // and an involutory mapping
4582 return orientation;
4583
4585 {
4586 AssertIndexRange(orientation, 6);
4587 constexpr std::array<types::geometric_orientation, 6> inverses{
4588 {0, 1, 4, 3, 2, 5}};
4589 return inverses[orientation];
4590 }
4592 {
4593 AssertIndexRange(orientation, 8);
4594 constexpr std::array<types::geometric_orientation, 8> inverses{
4595 {0, 1, 6, 3, 4, 5, 2, 7}};
4596 return inverses[orientation];
4597 }
4598 default:
4600 }
4601
4602 return std::numeric_limits<types::geometric_orientation>::max();
4603}
4604
4605
4607
4608#endif
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
std::size_t size() const
Definition array_view.h:737
Abstract base class for mapping classes.
Definition mapping.h:318
Definition point.h:111
static constexpr Point< dim, Number > unit_vector(const unsigned int i)
unsigned int unv_vertex_to_deal_vertex(const unsigned int vertex_n) const
const Quadrature< dim > & get_nodal_type_quadrature() const
unsigned int face_to_cell_lines(const unsigned int face, const unsigned int line, const types::geometric_orientation face_orientation) const
Tensor< 1, dim > unit_normal_vectors(const unsigned int face_no) const
static constexpr ndarray< unsigned int, 8, 4 > quadrilateral_vertex_permutations
constexpr bool is_simplex() const
double volume() const
constexpr unsigned int n_vertices() const
types::geometric_orientation compute_orientation(const std::array< T, N > &vertices_0, const std::array< T, N > &vertices_1) const
IsotropicRefinementChoice get_isotropic_refinement_choice(const unsigned int ref_choice) const
types::geometric_orientation face_to_cell_line_orientation(const unsigned int face_line_no, const unsigned int face_no, const types::geometric_orientation face_orientation, const types::geometric_orientation line_orientation) const
void serialize(Archive &archive, const unsigned int)
double face_measure(const unsigned int face_no) const
unsigned int child_cell_on_face(const unsigned int face, const unsigned int subface) const
constexpr bool is_hyper_cube() const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > line_indices() const
Tensor< 1, dim > unit_tangential_vectors(const unsigned int face_no, const unsigned int i) const
ArrayView< const unsigned int > faces_for_given_vertex(const unsigned int vertex_index) const
unsigned int vtk_quadratic_type() const
Point< dim > closest_point(const Point< dim > &p) const
Quadrature< dim > get_midpoint_quadrature() const
double d_linear_shape_function(const Point< dim > &xi, const unsigned int i) const
constexpr unsigned int get_dimension() const
static ReferenceCell n_vertices_to_type(const int dim_, const unsigned int n_vertices)
unsigned int vtk_lagrange_type() const
types::geometric_orientation get_combined_orientation(const ArrayView< const T > &vertices_0, const ArrayView< const T > &vertices_1) const
static constexpr types::geometric_orientation default_combined_face_orientation()
std::array< unsigned int, 2 > standard_vertex_to_face_and_vertex_index(const unsigned int vertex) const
ArrayView< const RefinementCase< dim > > refinement_cases() const
Quadrature< dim > get_gauss_type_quadrature(const unsigned n_points_1d) const
unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex) const
std::uint8_t kind
std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices() const
constexpr ReferenceCell< dim - 1 > face_reference_cell(const unsigned int face_index) const
friend std::istream & operator>>(std::istream &in, ReferenceCell< dim_ > &reference_cell)
unsigned int n_children(const RefinementCase< dim > ref_case=RefinementCase< dim >::isotropic_refinement) const
constexpr ::ndarray< unsigned int, 10, 8 > new_isotropic_child_cell_vertices(const unsigned int refinement_choice) const
constexpr unsigned int n_faces() const
constexpr ReferenceCell()
unsigned int standard_to_real_face_line(const unsigned int line, const unsigned int face, const types::geometric_orientation face_orientation) const
Point< dim > barycenter() const
constexpr unsigned int n_isotropic_refinement_choices() const
static constexpr ndarray< unsigned int, 8, 4 > quadrilateral_line_permutations
unsigned int vtk_vertex_to_deal_vertex(const unsigned int vertex_index) const
static constexpr ndarray< unsigned int, 3, 4 > line_vertex_permutations
constexpr unsigned int n_lines() const
Point< dim > vertex(const unsigned int v) const
unsigned int standard_to_real_face_vertex(const unsigned int vertex, const unsigned int face, const types::geometric_orientation face_orientation) const
std::pair< unsigned int, RefinementCase< dim - 1 > > equivalent_refinement_case(const types::geometric_orientation combined_face_orientation, const internal::SubfaceCase< dim > subface_case, const unsigned int subface_no) const
static constexpr ndarray< unsigned int, 6, 3 > triangle_line_permutations
unsigned int opposite_face_index(const unsigned int face_no) const
std::unique_ptr< Mapping< dim, spacedim > > get_default_mapping(const unsigned int degree) const
static constexpr ndarray< unsigned int, 1, 1 > vertex_vertex_permutations
Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i) const
constexpr operator std::uint8_t() const
std::string to_string() const
unsigned int n_face_orientations(const unsigned int face_no) const
unsigned int vtk_lexicographic_to_node_index(const std::array< unsigned, dim > &node_indices, const std::array< unsigned, dim > &nodes_per_direction, const bool legacy_format) const
static constexpr ndarray< unsigned int, 6, 4 > triangle_vertex_permutations
unsigned int gmsh_element_type() const
constexpr bool operator!=(const ReferenceCell< dim > &type) const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > isotropic_child_indices() const
constexpr ::ndarray< unsigned int, 13, 4 > new_isotropic_child_face_lines(const unsigned int refinement_choice) const
friend std::ostream & operator<<(std::ostream &out, const ReferenceCell< dim_ > &reference_cell)
types::geometric_orientation standard_vs_true_line_orientation(const unsigned int line, const unsigned int face, const types::geometric_orientation face_orientation, const types::geometric_orientation line_orientation) const
std::array< unsigned int, 2 > standard_line_to_face_and_line_index(const unsigned int line) const
Point< dim > face_vertex_location(const unsigned int face, const unsigned int vertex) const
constexpr ::ndarray< unsigned int, 13, 4, 2 > new_isotropic_child_face_line_vertices(const unsigned int refinement_choice) const
Tensor< 1, dim > face_tangent_vector(const unsigned int face_no, const unsigned int i) const
unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const types::geometric_orientation face_orientation) const
unsigned int exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices_by_type(const ReferenceCell< dim - 1 > &face_ref_type) const
const Mapping< dim, spacedim > & get_default_linear_mapping() const
static constexpr types::geometric_orientation reversed_combined_line_orientation()
boost::container::small_vector< T, 8 > permute_by_combined_orientation(const ArrayView< const T > &vertices, const types::geometric_orientation orientation) const
types::geometric_orientation get_inverse_combined_orientation(const types::geometric_orientation orientation) const
static constexpr std::size_t memory_consumption()
Tensor< 1, dim > face_normal_vector(const unsigned int face_no) const
constexpr ::ndarray< unsigned int, 10, 6 > new_isotropic_child_cell_faces(const unsigned int refinement_choice) const
constexpr unsigned int n_isotropic_children() const
Point< dim > subface_vertex_location(const unsigned int face_no, const unsigned int subface_no, const unsigned int subface_vertex_no, const RefinementCase< dim - 1 > face_refinement_case) const
bool contains_point(const Point< dim > &p, const double tolerance=0) const
unsigned int exodusii_face_to_deal_face(const unsigned int face_n) const
unsigned int face_and_line_to_cell_vertices(const unsigned int face, const unsigned int line, const unsigned int vertex, const types::geometric_orientation face_orientation) const
unsigned int vtk_linear_type() const
std::array< T, N > permute_according_orientation(const std::array< T, N > &vertices, const unsigned int orientation) const
constexpr bool operator==(const ReferenceCell< dim > &type) const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices() const
const ArrayView< const T > vertices_1
virtual void print_info(std::ostream &out) const override
const ReferenceCell< dim > entity_type
NoPermutation(const ReferenceCell< dim > &entity_type, const ArrayView< const T > &vertices_0, const ArrayView< const T > &vertices_1)
virtual ~NoPermutation() noexcept override=default
const ArrayView< const T > vertices_0
#define DEAL_II_DEPRECATED
Definition config.h:286
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
#define DEAL_II_DEPRECATED_WITH_COMMENT(comment)
Definition config.h:287
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
constexpr unsigned int max_n_children()
constexpr ReferenceCell< 3 > Hexahedron
constexpr ReferenceCell< 2 > Quadrilateral
ReferenceCell< dim > n_vertices_to_reference_cell(const unsigned int n_vertices)
constexpr ReferenceCell< 1 > Line
constexpr unsigned int max_n_lines()
constexpr const std::vector< ReferenceCell< dim > > get_reference_cells_in_dim()
constexpr ReferenceCell< 2 > Triangle
constexpr unsigned int max_n_vertices()
constexpr ReferenceCell< 3 > Tetrahedron
constexpr const ReferenceCell< dim > & get_simplex()
constexpr unsigned int max_n_faces()
constexpr ReferenceCell< dim > Invalid
constexpr ReferenceCell< 3 > Pyramid
constexpr const ReferenceCell< dim > & get_hypercube()
constexpr unsigned int n_reference_cells_in_dim()
constexpr ReferenceCell< 3 > Wedge
constexpr ReferenceCell< 0 > Vertex
std::tuple< bool, bool, bool > split_face_orientation(const types::geometric_orientation combined_orientation)
constexpr ReferenceCell< dim > make_reference_cell_from_int(const std::uint8_t kind)
constexpr types::geometric_orientation invalid_geometric_orientation
Definition types.h:377
constexpr unsigned int invalid_unsigned_int
Definition types.h:236
constexpr double SQRT1_2
Definition numbers.h:260
constexpr types::geometric_orientation reverse_line_orientation
Definition types.h:363
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:350
STL namespace.
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition types.h:30
std::uint8_t geometric_orientation
Definition types.h:38
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:105
IsotropicRefinementChoice
std::istream & operator>>(std::istream &in, ReferenceCell< dim > &reference_cell)
std::ostream & operator<<(std::ostream &out, const ReferenceCell< dim > &reference_cell)
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)
static constexpr ndarray< Tensor< 1, dim >, faces_per_cell, dim - 1 > unit_tangential_vectors
static unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static constexpr unsigned int vertices_per_cell
static constexpr ndarray< unsigned int, vertices_per_cell, dim > vertex_to_face
static constexpr unsigned int lines_per_cell
static std::array< unsigned int, 2 > standard_hex_line_to_quad_line_index(const unsigned int line)
static constexpr unsigned int faces_per_cell
static std::array< unsigned int, 2 > standard_hex_vertex_to_quad_vertex_index(const unsigned int vertex)
static std::array< unsigned int, 2 > standard_quad_vertex_to_line_vertex_index(const unsigned int vertex)
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
static constexpr std::array< unsigned int, faces_per_cell > opposite_face
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
static unsigned int face_to_cell_lines(const unsigned int face, const unsigned int line, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static constexpr std::array< Tensor< 1, dim >, faces_per_cell > unit_normal_vector
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
static constexpr unsigned int max_children_per_cell
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition tensor.h:2639
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition tensor.h:2664