Reference documentation for deal.II version GIT relicensing-480-geae235577e 2024-04-25 01:00:02+00:00
\(\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: LGPL-2.1-or-later
4// Copyright (C) 2020 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_tria_reference_cell_h
16#define dealii_tria_reference_cell_h
17
18#include <deal.II/base/config.h>
19
23#include <deal.II/base/point.h>
24#include <deal.II/base/tensor.h>
26
28
29#include <boost/container/small_vector.hpp>
30
31#include <iosfwd>
32#include <string>
33
35
36// Forward declarations
37#ifndef DOXYGEN
38template <int dim, int spacedim>
39class Mapping;
40
41template <int dim>
42class Quadrature;
43
44class ReferenceCell;
45#endif
46
47
48namespace internal
49{
62 constexpr ReferenceCell
63 make_reference_cell_from_int(const std::uint8_t kind);
64} // namespace internal
65
66
67
99{
100public:
108 static ReferenceCell
109 n_vertices_to_type(const int dim, const unsigned int n_vertices);
110
119 constexpr ReferenceCell();
120
131 bool
132 is_hyper_cube() const;
133
137 bool
138 is_simplex() const;
139
144 unsigned int
145 get_dimension() const;
146
160 template <int dim>
161 double
162 d_linear_shape_function(const Point<dim> &xi, const unsigned int i) const;
163
168 template <int dim>
171 const unsigned int i) const;
172
182 template <int dim, int spacedim = dim>
183 std::unique_ptr<Mapping<dim, spacedim>>
184 get_default_mapping(const unsigned int degree) const;
185
197 template <int dim, int spacedim = dim>
200
209 template <int dim>
211 get_gauss_type_quadrature(const unsigned n_points_1d) const;
212
222 template <int dim>
225
239 template <int dim>
240 const Quadrature<dim> &
242
257 unsigned int
258 n_vertices() const;
259
265 vertex_indices() const;
266
274 template <int dim>
276 vertex(const unsigned int v) const;
277
283 unsigned int
284 n_lines() const;
285
291 line_indices() const;
292
298 unsigned int
299 n_faces() const;
300
306 face_indices() const;
307
319 unsigned int
320 n_isotropic_children() const;
321
328
342 face_reference_cell(const unsigned int face_no) const;
343
357 static constexpr unsigned char
359
370 static constexpr unsigned char
372
414 unsigned int
415 child_cell_on_face(const unsigned int face,
416 const unsigned int subface,
417 const unsigned char face_orientation =
419
429 std::array<unsigned int, 2>
430 standard_vertex_to_face_and_vertex_index(const unsigned int vertex) const;
431
441 std::array<unsigned int, 2>
442 standard_line_to_face_and_line_index(const unsigned int line) const;
443
452 unsigned int
453 line_to_cell_vertices(const unsigned int line,
454 const unsigned int vertex) const;
455
459 unsigned int
460 face_to_cell_lines(const unsigned int face,
461 const unsigned int line,
462 const unsigned char face_orientation) const;
463
467 unsigned int
468 face_to_cell_vertices(const unsigned int face,
469 const unsigned int vertex,
470 const unsigned char face_orientation) const;
471
491 template <int dim>
493 face_vertex_location(const unsigned int face,
494 const unsigned int vertex) const;
495
499 unsigned int
500 standard_to_real_face_vertex(const unsigned int vertex,
501 const unsigned int face,
502 const unsigned char face_orientation) const;
503
507 unsigned int
508 standard_to_real_face_line(const unsigned int line,
509 const unsigned int face,
510 const unsigned char face_orientation) const;
511
522 bool
523 standard_vs_true_line_orientation(const unsigned int line,
524 const unsigned int face,
525 const unsigned char face_orientation,
526 const bool line_orientation) const;
527
551 double
552 volume() const;
553
566 template <int dim>
568 barycenter() const;
569
589 template <int dim>
590 bool
591 contains_point(const Point<dim> &p, const double tolerance = 0) const;
592
597 template <int dim>
599 closest_point(const Point<dim> &p) const;
600
608 template <int dim>
610 unit_tangential_vectors(const unsigned int face_no,
611 const unsigned int i) const;
612
616 template <int dim>
618 unit_normal_vectors(const unsigned int face_no) const;
619
625 unsigned int
626 n_face_orientations(const unsigned int face_no) const;
627
644 template <typename T, std::size_t N>
645 DEAL_II_DEPRECATED unsigned char
646 compute_orientation(const std::array<T, N> &vertices_0,
647 const std::array<T, N> &vertices_1) const;
648
674 template <typename T>
675 unsigned char
677 const ArrayView<const T> &vertices_1) const;
678
679
693 template <typename T, std::size_t N>
694 DEAL_II_DEPRECATED std::array<T, N>
695 permute_according_orientation(const std::array<T, N> &vertices,
696 const unsigned int orientation) const;
697
709 template <typename T>
710 boost::container::small_vector<T, 8>
712 const unsigned char orientation) const;
713
719 unsigned char
720 get_inverse_combined_orientation(const unsigned char orientation) const;
721
726 faces_for_given_vertex(const unsigned int vertex_index) const;
727
740 unsigned int
741 exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const;
742
746 unsigned int
747 exodusii_face_to_deal_face(const unsigned int face_n) const;
748
752 unsigned int
753 unv_vertex_to_deal_vertex(const unsigned int vertex_n) const;
754
758 unsigned int
759 vtk_linear_type() const;
760
765 unsigned int
766 vtk_quadratic_type() const;
767
772 unsigned int
773 vtk_lagrange_type() const;
774
787 template <int dim>
788 unsigned int
790 const std::array<unsigned, dim> &node_indices,
791 const std::array<unsigned, dim> &nodes_per_direction,
792 const bool legacy_format) const;
793
797 unsigned int
798 vtk_vertex_to_deal_vertex(const unsigned int vertex_index) const;
799
803 unsigned int
804 gmsh_element_type() const;
805
819 std::string
820 to_string() const;
821
825 constexpr operator std::uint8_t() const;
826
830 constexpr bool
831 operator==(const ReferenceCell &type) const;
832
836 constexpr bool
837 operator!=(const ReferenceCell &type) const;
838
844 template <class Archive>
845 void
846 serialize(Archive &archive, const unsigned int /*version*/);
847
851 static constexpr std::size_t
853
858private:
862 std::uint8_t kind;
863
870 constexpr ReferenceCell(const std::uint8_t kind);
871
878 {{{0}}}};
879
884 {{{1, 0}}, {{0, 1}}}};
885
890 {{{0, 2, 1}},
891 {{0, 1, 2}},
892 {{2, 1, 0}},
893 {{2, 0, 1}},
894 {{1, 0, 2}},
895 {{1, 2, 0}}}};
896
900 static constexpr ndarray<unsigned int, 8, 4>
902 {{0, 1, 2, 3}},
903 {{2, 3, 0, 1}},
904 {{2, 0, 3, 1}},
905 {{3, 1, 2, 0}},
906 {{3, 2, 1, 0}},
907 {{1, 0, 3, 2}},
908 {{1, 3, 0, 2}}}};
909
914 friend constexpr ReferenceCell
916
917 friend std::ostream &
918 operator<<(std::ostream &out, const ReferenceCell &reference_cell);
919
920 friend std::istream &
921 operator>>(std::istream &in, ReferenceCell &reference_cell);
922};
923
924
932std::ostream &
933operator<<(std::ostream &out, const ReferenceCell &reference_cell);
934
942std::istream &
943operator>>(std::istream &in, ReferenceCell &reference_cell);
944
945
946inline constexpr ReferenceCell::ReferenceCell(const std::uint8_t kind)
947 : kind(kind)
948{}
949
950
951
952inline constexpr ReferenceCell::operator std::uint8_t() const
953{
954 return kind;
955}
956
957
958
959inline constexpr bool
961{
962 return kind == type.kind;
963}
964
965
966
967inline constexpr bool
969{
970 return kind != type.kind;
971}
972
973
974
975namespace internal
976{
977 inline constexpr ReferenceCell
978 make_reference_cell_from_int(const std::uint8_t kind)
979 {
980#ifndef DEAL_II_CXX14_CONSTEXPR_BUG
981 // Make sure these are the only indices from which objects can be
982 // created.
983 Assert((kind == std::numeric_limits<std::uint8_t>::max()) || (kind < 8),
985#endif
986
987 // Call the private constructor, which we can from here because this
988 // function is a 'friend'.
989 return {kind};
990 }
991} // namespace internal
992
993
994
1003{
1004 constexpr const ReferenceCell Vertex =
1006 constexpr const ReferenceCell Line =
1008 constexpr const ReferenceCell Triangle =
1012 constexpr const ReferenceCell Tetrahedron =
1014 constexpr const ReferenceCell Pyramid =
1016 constexpr const ReferenceCell Wedge =
1018 constexpr const ReferenceCell Hexahedron =
1020 constexpr const ReferenceCell Invalid =
1022 std::numeric_limits<std::uint8_t>::max());
1023
1029 template <int dim>
1030 constexpr const ReferenceCell &
1031 get_simplex();
1032
1038 template <int dim>
1039 constexpr const ReferenceCell &
1040 get_hypercube();
1041} // namespace ReferenceCells
1042
1043
1044
1046 : ReferenceCell(ReferenceCells::Invalid)
1047{}
1048
1049
1050
1051template <class Archive>
1052inline void
1053ReferenceCell::serialize(Archive &archive, const unsigned int /*version*/)
1054{
1055 archive &kind;
1056}
1057
1058
1059
1060inline constexpr std::size_t
1065
1066
1067
1069ReferenceCell::faces_for_given_vertex(const unsigned int vertex) const
1070{
1072 switch (this->kind)
1073 {
1075 return {&GeometryInfo<2>::vertex_to_face[vertex][0], 1};
1077 return {&GeometryInfo<2>::vertex_to_face[vertex][0], 2};
1079 {
1080 static constexpr ndarray<unsigned int, 3, 2> table = {
1081 {{{0, 2}}, {{0, 1}}, {{1, 2}}}};
1082 return table[vertex];
1083 }
1085 {
1086 static constexpr ndarray<unsigned int, 4, 3> table = {
1087 {{{0, 1, 2}}, {{0, 1, 3}}, {{0, 2, 3}}, {{1, 2, 3}}}};
1088
1089 return table[vertex];
1090 }
1092 {
1093 static constexpr unsigned int X = numbers::invalid_unsigned_int;
1094 static constexpr ndarray<unsigned int, 5, 4> table = {
1095 {{{0, 1, 3, X}},
1096 {{0, 2, 3, X}},
1097 {{0, 1, 4, X}},
1098 {{0, 2, 4, X}},
1099 {{1, 2, 3, 4}}}};
1100
1101 return {&table[vertex][0], vertex == 4 ? 4u : 3u};
1102 }
1104 {
1106 static constexpr ndarray<unsigned int, 6, 3> table = {{{{0, 2, 4}},
1107 {{0, 2, 3}},
1108 {{0, 3, 4}},
1109 {{1, 2, 4}},
1110 {{1, 2, 3}},
1111 {{1, 3, 4}}}};
1112
1113 return table[vertex];
1114 }
1116 return {&GeometryInfo<3>::vertex_to_face[vertex][0], 3};
1117 default:
1119 }
1120
1121 return {};
1122}
1123
1124
1125
1126inline bool
1128{
1129 return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
1132}
1133
1134
1135
1136inline bool
1138{
1139 return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
1140 *this == ReferenceCells::Triangle ||
1142}
1143
1144
1145
1146inline unsigned int
1148{
1149 switch (this->kind)
1150 {
1152 return 0;
1154 return 1;
1157 return 2;
1162 return 3;
1163 default:
1165 }
1166
1168}
1169
1170
1171
1172template <int dim>
1175{
1176 return Quadrature<dim>(std::vector<Point<dim>>({barycenter<dim>()}),
1177 std::vector<double>({volume()}));
1178}
1179
1180
1181
1182inline unsigned int
1184{
1185 switch (this->kind)
1186 {
1188 return 1;
1190 return 2;
1192 return 3;
1194 return 4;
1196 return 4;
1198 return 5;
1200 return 6;
1202 return 8;
1203 default:
1205 }
1206
1208}
1209
1210
1211
1212inline unsigned int
1214{
1215 switch (this->kind)
1216 {
1218 return 0;
1220 return 1;
1222 return 3;
1224 return 4;
1226 return 6;
1228 return 8;
1230 return 9;
1232 return 12;
1233 default:
1235 }
1236
1238}
1239
1240
1241
1242template <int dim>
1244ReferenceCell::vertex(const unsigned int v) const
1245{
1248
1249 switch (dim)
1250 {
1251 case 0:
1252 {
1253 if (*this == ReferenceCells::Vertex)
1254 return Point<dim>();
1255 break;
1256 }
1257 case 1:
1258 {
1259 static const Point<dim> vertices[2] = {
1260 Point<dim>(), // the origin
1261 Point<dim>::unit_vector(0) // unit point along x-axis
1262 };
1263 if (*this == ReferenceCells::Line)
1264 return vertices[v];
1265 break;
1266 }
1267 case 2:
1268 {
1269 switch (this->kind)
1270 {
1272 {
1273 static const Point<dim> vertices[3] = {
1274 Point<dim>(), // the origin
1275 Point<dim>::unit_vector(0), // unit point along x-axis
1276 Point<dim>::unit_vector(1) // unit point along y-axis
1277 };
1278 return vertices[v];
1279 }
1281 {
1282 static const Point<dim> vertices[4] = {
1283 // First the two points on the x-axis
1284 Point<dim>(),
1286 // Then these two points shifted in the y-direction
1289 return vertices[v];
1290 }
1291 }
1292 break;
1293 }
1294 case 3:
1295 {
1296 switch (this->kind)
1297 {
1299 {
1300 static const Point<dim> vertices[4] = {
1301 Point<dim>(), // the origin
1302 Point<dim>::unit_vector(0), // unit point along x-axis
1303 Point<dim>::unit_vector(1), // unit point along y-axis
1304 Point<dim>::unit_vector(2) // unit point along z-axis
1305 };
1306 return vertices[v];
1307 }
1309 {
1310 static const Point<dim> vertices[5] = {
1311 Point<dim>{-1.0, -1.0, 0.0},
1312 Point<dim>{+1.0, -1.0, 0.0},
1313 Point<dim>{-1.0, +1.0, 0.0},
1314 Point<dim>{+1.0, +1.0, 0.0},
1315 Point<dim>{+0.0, +0.0, 1.0}};
1316 return vertices[v];
1317 }
1319 {
1320 static const Point<dim> vertices[6] = {
1321 // First the three points on the triangular base of the
1322 // wedge:
1323 Point<dim>(),
1326 // And now everything shifted in the z-direction again
1330 return vertices[v];
1331 }
1333 {
1334 static const Point<dim> vertices[8] = {
1335 // First the two points on the x-axis
1336 Point<dim>(),
1338 // Then these two points shifted in the y-direction
1341 // And now all four points shifted in the z-direction
1348 return vertices[v];
1349 }
1350 }
1351 break;
1352 }
1353 default:
1355 }
1356
1358 return Point<dim>();
1359}
1360
1361
1362inline unsigned int
1364{
1365 switch (this->kind)
1366 {
1368 return 0;
1370 return 2;
1372 return 3;
1374 return 4;
1376 return 4;
1378 return 5;
1380 return 5;
1382 return 6;
1383 default:
1385 }
1386
1388}
1389
1390
1391
1394{
1395 return {0U, n_faces()};
1396}
1397
1398
1399
1400inline unsigned int
1402{
1403 switch (this->kind)
1404 {
1406 return 0;
1408 return 2;
1410 return 4;
1412 return 4;
1414 return 8;
1416 // We haven't yet decided how to refine pyramids. Update this when we
1417 // have
1421 return 8;
1423 return 8;
1424 default:
1426 }
1427
1429}
1430
1431
1432
1435{
1436 return {0U, n_isotropic_children()};
1437}
1438
1439
1440
1443{
1444 return {0U, n_vertices()};
1445}
1446
1447
1448
1451{
1452 return {0U, n_lines()};
1453}
1454
1455
1456
1457inline ReferenceCell
1458ReferenceCell::face_reference_cell(const unsigned int face_no) const
1459{
1460 AssertIndexRange(face_no, n_faces());
1461
1462 switch (this->kind)
1463 {
1470 return ReferenceCells::Line;
1474 if (face_no == 0)
1476 else
1479 if (face_no > 1)
1481 else
1485 default:
1487 }
1488
1490}
1491
1492
1493
1494inline constexpr unsigned char
1496{
1497 // Our convention is that 'orientation' has a default value of true and
1498 // occupies the least-significant bit while rotate and flip have default
1499 // values of 'false' and occupy the second and third bits.
1500 return 0b001;
1501}
1502
1503
1504
1505inline constexpr unsigned char
1507{
1508 // For a reversed line 'orientation' is false and neither flip nor rotate are
1509 // defined.
1510 return 0b000;
1511}
1512
1513
1514
1515inline unsigned int
1517 const unsigned int face,
1518 const unsigned int subface,
1519 const unsigned char combined_face_orientation) const
1520{
1521 AssertIndexRange(face, n_faces());
1523
1524 switch (this->kind)
1525 {
1528 {
1530 break;
1531 }
1533 {
1534 static constexpr ndarray<unsigned int, 3, 2> subcells = {
1535 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
1536
1537 return subcells[face][subface];
1538 }
1540 {
1541 const auto [face_orientation, face_rotation, face_flip] =
1542 internal::split_face_orientation(combined_face_orientation);
1543
1546 face,
1547 subface,
1548 face_orientation,
1549 face_flip,
1550 face_rotation);
1551 }
1555 {
1557 break;
1558 }
1560 {
1561 const auto [face_orientation, face_rotation, face_flip] =
1562 internal::split_face_orientation(combined_face_orientation);
1563
1566 face,
1567 subface,
1568 face_orientation,
1569 face_flip,
1570 face_rotation);
1571 }
1572 default:
1574 }
1575
1577}
1578
1579
1580
1581inline std::array<unsigned int, 2>
1583 const unsigned int vertex) const
1584{
1586 // Work around a GCC warning at higher optimization levels by making all of
1587 // these tables the same size
1588 constexpr unsigned int X = numbers::invalid_unsigned_int;
1589
1590 switch (this->kind)
1591 {
1595 break;
1597 {
1598 static constexpr ndarray<unsigned int, 6, 2> table = {
1599 {{{0, 0}}, {{0, 1}}, {{1, 1}}, {{X, X}}, {{X, X}}, {{X, X}}}};
1600
1601 return table[vertex];
1602 }
1604 {
1606 vertex);
1607 }
1609 {
1610 static constexpr ndarray<unsigned int, 6, 2> table = {
1611 {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 2}}, {{X, X}}, {{X, X}}}};
1612
1613 return table[vertex];
1614 }
1616 {
1617 static constexpr ndarray<unsigned int, 6, 2> table = {
1618 {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{0, 3}}, {{1, 2}}, {{X, X}}}};
1619
1620 return table[vertex];
1621 }
1623 {
1624 static constexpr ndarray<unsigned int, 6, 2> table = {
1625 {{{0, 1}}, {{0, 0}}, {{0, 2}}, {{1, 0}}, {{1, 1}}, {{1, 2}}}};
1626
1627 return table[vertex];
1628 }
1630 {
1632 vertex);
1633 }
1634 default:
1636 }
1637
1638 return {};
1639}
1640
1641
1642
1643inline std::array<unsigned int, 2>
1645 const unsigned int line) const
1646{
1647 AssertIndexRange(line, n_lines());
1648
1649 switch (this->kind)
1650 {
1655 {
1657 break;
1658 }
1660 {
1661 static const std::array<unsigned int, 2> table[6] = {
1662 {{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}};
1663
1664 return table[line];
1665 }
1667 {
1668 static const std::array<unsigned int, 2> table[8] = {{{0, 0}},
1669 {{0, 1}},
1670 {{0, 2}},
1671 {{0, 3}},
1672 {{1, 2}},
1673 {{2, 1}},
1674 {{1, 1}},
1675 {{2, 2}}};
1676
1677 return table[line];
1678 }
1680 {
1681 static const std::array<unsigned int, 2> table[9] = {{{0, 0}},
1682 {{0, 2}},
1683 {{0, 1}},
1684 {{1, 0}},
1685 {{1, 1}},
1686 {{1, 2}},
1687 {{2, 0}},
1688 {{2, 1}},
1689 {{3, 1}}};
1690
1691 return table[line];
1692 }
1694 {
1696 }
1697 default:
1699 }
1700
1701 return {};
1702}
1703
1704
1705inline unsigned int
1707 const unsigned int vertex) const
1708{
1710 AssertIndexRange(line, n_lines());
1711
1712 switch (this->kind)
1713 {
1716 return vertex;
1718 {
1719 static constexpr ndarray<unsigned int, 3, 2> table = {
1720 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
1721 return table[line][vertex];
1722 }
1724 {
1725 static constexpr ndarray<unsigned int, 4, 2> table = {
1726 {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
1727 return table[line][vertex];
1728 }
1730 {
1731 static constexpr ndarray<unsigned int, 6, 2> table = {
1732 {{{0, 1}}, {{1, 2}}, {{2, 0}}, {{0, 3}}, {{1, 3}}, {{2, 3}}}};
1733 return table[line][vertex];
1734 }
1736 {
1737 static constexpr ndarray<unsigned int, 8, 2> table = {{{{0, 2}},
1738 {{1, 3}},
1739 {{0, 1}},
1740 {{2, 3}},
1741 {{4, 0}},
1742 {{1, 4}},
1743 {{2, 4}},
1744 {{4, 3}}}};
1745 return table[line][vertex];
1746 }
1748 {
1749 static constexpr ndarray<unsigned int, 9, 2> table = {{{{1, 0}},
1750 {{2, 1}},
1751 {{0, 2}},
1752 {{3, 4}},
1753 {{4, 5}},
1754 {{5, 3}},
1755 {{0, 3}},
1756 {{1, 4}},
1757 {{2, 5}}}};
1758 return table[line][vertex];
1759 }
1761 {
1762 // first four lines comprise the bottom face, next four are the top,
1763 // and the last four are 'bottom to top'
1764 static constexpr ndarray<unsigned int, 12, 2> table = {{{{0, 2}},
1765 {{1, 3}},
1766 {{0, 1}},
1767 {{2, 3}},
1768 {{4, 6}},
1769 {{5, 7}},
1770 {{4, 5}},
1771 {{6, 7}},
1772 {{0, 4}},
1773 {{1, 5}},
1774 {{2, 6}},
1775 {{3, 7}}}};
1776 return table[line][vertex];
1777 }
1778
1779 default:
1781 }
1782
1784}
1785
1786
1787inline unsigned int
1789 const unsigned int face,
1790 const unsigned int line,
1791 const unsigned char combined_face_orientation) const
1792{
1793 AssertIndexRange(face, n_faces());
1795
1796 static constexpr unsigned int X = numbers::invalid_unsigned_int;
1797
1798 switch (this->kind)
1799 {
1801 {
1803 break;
1804 }
1806 {
1807 const auto [face_orientation, face_rotation, face_flip] =
1808 internal::split_face_orientation(combined_face_orientation);
1809
1811 face, line, face_orientation, face_flip, face_rotation);
1812 }
1814 {
1815 return face;
1816 }
1818 {
1819 const auto [face_orientation, face_rotation, face_flip] =
1820 internal::split_face_orientation(combined_face_orientation);
1821
1823 face, line, face_orientation, face_flip, face_rotation);
1824 }
1826 {
1827 static constexpr ndarray<unsigned int, 4, 3> table = {
1828 {{{0, 1, 2}}, {{0, 3, 4}}, {{2, 5, 3}}, {{1, 4, 5}}}};
1829
1830 return table[face][standard_to_real_face_line(
1831 line, face, combined_face_orientation)];
1832 }
1834 {
1835 static constexpr ndarray<unsigned int, 5, 4> table = {
1836 {{{0, 1, 2, 3}},
1837 {{0, 6, 4, X}},
1838 {{1, 5, 7, X}},
1839 {{2, 4, 5, X}},
1840 {{3, 7, 6, X}}}};
1841
1842 return table[face][standard_to_real_face_line(
1843 line, face, combined_face_orientation)];
1844 }
1846 {
1847 static constexpr ndarray<unsigned int, 5, 4> table = {
1848 {{{0, 2, 1, X}},
1849 {{3, 4, 5, X}},
1850 {{6, 7, 0, 3}},
1851 {{7, 8, 1, 4}},
1852 {{8, 6, 5, 2}}}};
1853
1854 return table[face][standard_to_real_face_line(
1855 line, face, combined_face_orientation)];
1856 }
1858 {
1859 const auto [face_orientation, face_rotation, face_flip] =
1860 internal::split_face_orientation(combined_face_orientation);
1861
1863 face, line, face_orientation, face_flip, face_rotation);
1864 }
1865 default:
1867 }
1868
1870}
1871
1872
1873
1874inline unsigned int
1876 const unsigned int face,
1877 const unsigned int vertex,
1878 const unsigned char combined_face_orientation) const
1879{
1880 AssertIndexRange(face, n_faces());
1882 // TODO: once the default orientation is switched to 0 then we can remove this
1883 // special case for 1D.
1884 if (get_dimension() == 1)
1885 Assert(combined_face_orientation ==
1887 ExcMessage("In 1D, all faces must have the default orientation."));
1888 else
1889 AssertIndexRange(combined_face_orientation, n_face_orientations(face));
1890
1891 switch (this->kind)
1892 {
1894 {
1896 break;
1897 }
1899 {
1900 const auto [face_orientation, face_rotation, face_flip] =
1901 internal::split_face_orientation(combined_face_orientation);
1902
1904 face, vertex, face_orientation, face_flip, face_rotation);
1905 }
1907 {
1908 static constexpr ndarray<unsigned int, 3, 2> table = {
1909 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
1910
1911 return table[face][combined_face_orientation ==
1913 vertex :
1914 (1 - vertex)];
1915 }
1917 {
1918 const auto [face_orientation, face_rotation, face_flip] =
1919 internal::split_face_orientation(combined_face_orientation);
1920
1922 face, vertex, face_orientation, face_flip, face_rotation);
1923 }
1925 {
1926 static constexpr ndarray<unsigned int, 4, 3> table = {
1927 {{{0, 1, 2}}, {{1, 0, 3}}, {{0, 2, 3}}, {{2, 1, 3}}}};
1928
1929 return table[face][standard_to_real_face_vertex(
1930 vertex, face, combined_face_orientation)];
1931 }
1933 {
1934 constexpr auto X = numbers::invalid_unsigned_int;
1935 static constexpr ndarray<unsigned int, 5, 4> table = {
1936 {{{0, 1, 2, 3}},
1937 {{0, 2, 4, X}},
1938 {{3, 1, 4, X}},
1939 {{1, 0, 4, X}},
1940 {{2, 3, 4, X}}}};
1941
1942 return table[face][standard_to_real_face_vertex(
1943 vertex, face, combined_face_orientation)];
1944 }
1946 {
1947 constexpr auto X = numbers::invalid_unsigned_int;
1948 static constexpr ndarray<unsigned int, 6, 4> table = {
1949 {{{1, 0, 2, X}},
1950 {{3, 4, 5, X}},
1951 {{0, 1, 3, 4}},
1952 {{1, 2, 4, 5}},
1953 {{2, 0, 5, 3}}}};
1954
1955 return table[face][standard_to_real_face_vertex(
1956 vertex, face, combined_face_orientation)];
1957 }
1959 {
1960 const auto [face_orientation, face_rotation, face_flip] =
1961 internal::split_face_orientation(combined_face_orientation);
1962
1964 face, vertex, face_orientation, face_flip, face_rotation);
1965 }
1966 default:
1968 }
1969
1971}
1972
1973
1974
1975template <int dim>
1978 const unsigned int vertex) const
1979{
1980 return this->template vertex<dim>(
1982}
1983
1984
1985
1986inline unsigned int
1988 const unsigned int vertex,
1989 const unsigned int face,
1990 const unsigned char face_orientation) const
1991{
1992 AssertIndexRange(face, n_faces());
1994
1995 switch (this->kind)
1996 {
1999 break;
2001 Assert(face_orientation == default_combined_face_orientation(),
2002 ExcMessage(
2003 "In 1D, all faces must have the default orientation."));
2004 return vertex;
2007 return line_vertex_permutations[face_orientation][vertex];
2009 return triangle_vertex_permutations[face_orientation][vertex];
2011 // face 0 is a quadrilateral
2012 if (face == 0)
2013 return quadrilateral_vertex_permutations[face_orientation][vertex];
2014 else
2015 return triangle_vertex_permutations[face_orientation][vertex];
2017 // faces 0 and 1 are triangles
2018 if (face > 1)
2019 return quadrilateral_vertex_permutations[face_orientation][vertex];
2020 else
2021 return triangle_vertex_permutations[face_orientation][vertex];
2023 return quadrilateral_vertex_permutations[face_orientation][vertex];
2024 default:
2026 }
2027
2030}
2031
2032
2033
2034inline unsigned int
2036 const unsigned int line,
2037 const unsigned int face,
2038 const unsigned char combined_face_orientation) const
2039{
2040 AssertIndexRange(face, n_faces());
2042
2043 static constexpr ndarray<unsigned int, 6, 3> triangle_table = {{{{2, 1, 0}},
2044 {{0, 1, 2}},
2045 {{1, 0, 2}},
2046 {{2, 0, 1}},
2047 {{0, 2, 1}},
2048 {{1, 2, 0}}}};
2049
2050
2051 switch (this->kind)
2052 {
2058 break;
2060 return triangle_table[combined_face_orientation][line];
2062 if (face == 0) // The quadrilateral face
2063 {
2064 const auto [face_orientation, face_rotation, face_flip] =
2065 internal::split_face_orientation(combined_face_orientation);
2066
2068 face_orientation,
2069 face_flip,
2070 face_rotation);
2071 }
2072 else // One of the triangular faces
2073 {
2074 return triangle_table[combined_face_orientation][line];
2075 }
2077 if (face > 1) // One of the quadrilateral faces
2078 {
2079 const auto [face_orientation, face_rotation, face_flip] =
2080 internal::split_face_orientation(combined_face_orientation);
2081
2083 face_orientation,
2084 face_flip,
2085 face_rotation);
2086 }
2087 else // One of the triangular faces
2088 return triangle_table[combined_face_orientation][line];
2090 {
2091 static constexpr ndarray<unsigned int, 8, 4> table = {
2092 {{{2, 3, 0, 1}},
2093 {{0, 1, 2, 3}},
2094 {{0, 1, 3, 2}},
2095 {{3, 2, 0, 1}},
2096 {{3, 2, 1, 0}},
2097 {{1, 0, 3, 2}},
2098 {{1, 0, 2, 3}},
2099 {{2, 3, 1, 0}}}};
2100 return table[combined_face_orientation][line];
2101 }
2102 default:
2104 }
2105
2107}
2108
2109
2110
2111namespace ReferenceCells
2112{
2113 template <int dim>
2114 inline constexpr const ReferenceCell &
2116 {
2117 switch (dim)
2118 {
2119 case 0:
2121 case 1:
2122 return ReferenceCells::Line;
2123 case 2:
2125 case 3:
2127 default:
2129 }
2131 }
2132
2133
2134
2135 template <int dim>
2136 inline constexpr const ReferenceCell &
2138 {
2139 switch (dim)
2140 {
2141 case 0:
2143 case 1:
2144 return ReferenceCells::Line;
2145 case 2:
2147 case 3:
2149 default:
2151 }
2153 }
2154} // namespace ReferenceCells
2155
2156
2157inline ReferenceCell
2158ReferenceCell::n_vertices_to_type(const int dim, const unsigned int n_vertices)
2159{
2160 AssertIndexRange(dim, 4);
2162
2163 const auto X = ReferenceCells::Invalid;
2164 static const ndarray<ReferenceCell, 4, 9> table = {
2165 {// dim 0
2166 {{X, ReferenceCells::Vertex, X, X, X, X, X, X, X}},
2167 // dim 1
2168 {{X, X, ReferenceCells::Line, X, X, X, X, X, X}},
2169 // dim 2
2170 {{X,
2171 X,
2172 X,
2175 X,
2176 X,
2177 X,
2178 X}},
2179 // dim 3
2180 {{X,
2181 X,
2182 X,
2183 X,
2187 X,
2190 ExcMessage("The combination of dim = " + std::to_string(dim) +
2191 " and n_vertices = " + std::to_string(n_vertices) +
2192 " does not correspond to a known reference cell type."));
2193 return table[dim][n_vertices];
2194}
2195
2196
2197
2198template <int dim>
2199inline double
2201 const unsigned int i) const
2202{
2205 switch (this->kind)
2206 {
2212 // see also BarycentricPolynomials<2>::compute_value
2214 {
2215 switch (i)
2216 {
2217 case 0:
2218 return 1.0 - xi[std::min(0, dim - 1)] -
2219 xi[std::min(1, dim - 1)];
2220 case 1:
2221 return xi[std::min(0, dim - 1)];
2222 case 2:
2223 return xi[std::min(1, dim - 1)];
2224 default:
2226 }
2227 }
2228 // see also BarycentricPolynomials<3>::compute_value
2230 {
2231 switch (i)
2232 {
2233 case 0:
2234 return 1.0 - xi[std::min(0, dim - 1)] -
2235 xi[std::min(1, dim - 1)] - xi[std::min(2, dim - 1)];
2236 case 1:
2237 return xi[std::min(0, dim - 1)];
2238 case 2:
2239 return xi[std::min(1, dim - 1)];
2240 case 3:
2241 return xi[std::min(2, dim - 1)];
2242 default:
2244 }
2245 }
2246 // see also ScalarLagrangePolynomialPyramid::compute_value()
2248 {
2249 const double Q14 = 0.25;
2250
2251 const double r = xi[std::min(0, dim - 1)];
2252 const double s = xi[std::min(1, dim - 1)];
2253 const double t = xi[std::min(2, dim - 1)];
2254
2255 const double ratio =
2256 (std::fabs(t - 1.0) > 1.0e-14 ? (r * s * t) / (1.0 - t) : 0.0);
2257
2258 if (i == 0)
2259 return Q14 * ((1.0 - r) * (1.0 - s) - t + ratio);
2260 if (i == 1)
2261 return Q14 * ((1.0 + r) * (1.0 - s) - t - ratio);
2262 if (i == 2)
2263 return Q14 * ((1.0 - r) * (1.0 + s) - t - ratio);
2264 if (i == 3)
2265 return Q14 * ((1.0 + r) * (1.0 + s) - t + ratio);
2266 else
2267 return t;
2268 }
2269 // see also ScalarLagrangePolynomialWedge::compute_value()
2272 .d_linear_shape_function<2>(Point<2>(xi[std::min(0, dim - 1)],
2273 xi[std::min(1, dim - 1)]),
2274 i % 3) *
2276 .d_linear_shape_function<1>(Point<1>(xi[std::min(2, dim - 1)]),
2277 i / 3);
2278 default:
2280 }
2281
2282 return 0.0;
2283}
2284
2285
2286
2287template <int dim>
2288inline Tensor<1, dim>
2290 const unsigned int i) const
2291{
2293 switch (this->kind)
2294 {
2300 // see also BarycentricPolynomials<2>::compute_grad()
2302 switch (i)
2303 {
2304 case 0:
2305 return Point<dim>(-1.0, -1.0);
2306 case 1:
2307 return Point<dim>(+1.0, +0.0);
2308 case 2:
2309 return Point<dim>(+0.0, +1.0);
2310 default:
2312 }
2313 default:
2315 }
2316
2317 return Point<dim>(+0.0, +0.0, +0.0);
2318}
2319
2320
2321
2322inline double
2324{
2325 switch (this->kind)
2326 {
2328 return 0;
2330 return 1;
2332 return 1. / 2.;
2334 return 1;
2336 return 1. / 6.;
2338 return 4. / 3.;
2340 return 1. / 2.;
2342 return 1;
2343 default:
2345 }
2346
2347 return 0.0;
2348}
2349
2350
2351
2352template <int dim>
2353inline Point<dim>
2355{
2357
2358 switch (this->kind)
2359 {
2361 return Point<dim>();
2363 return Point<dim>(1. / 2.);
2365 return Point<dim>(1. / 3., 1. / 3.);
2367 return Point<dim>(1. / 2., 1. / 2.);
2369 return Point<dim>(1. / 4., 1. / 4., 1. / 4.);
2371 return Point<dim>(0, 0, 1. / 4.);
2373 return Point<dim>(1. / 3, 1. / 3, 1. / 2.);
2375 return Point<dim>(1. / 2., 1. / 2., 1. / 2.);
2376 default:
2378 }
2379
2380 return Point<dim>();
2381}
2382
2383
2384
2385template <int dim>
2386inline bool
2387ReferenceCell::contains_point(const Point<dim> &p, const double tolerance) const
2388{
2390
2391 // Introduce abbreviations to silence compiler warnings about invalid
2392 // array accesses (that can't happen, of course, but the compiler
2393 // doesn't know that).
2394 constexpr unsigned int x_coordinate = 0;
2395 constexpr unsigned int y_coordinate = (dim >= 2 ? 1 : 0);
2396 constexpr unsigned int z_coordinate = (dim >= 3 ? 2 : 0);
2397
2398 switch (this->kind)
2399 {
2401 {
2402 // Vertices are special cases in that they do not actually
2403 // have coordinates. Error out if this function is called
2404 // with a vertex:
2405 Assert(false,
2406 ExcMessage("Vertices are zero-dimensional objects and "
2407 "as a consequence have no coordinates. You "
2408 "cannot meaningfully ask whether a point is "
2409 "inside a vertex (within a certain tolerance) "
2410 "without coordinate values."));
2411 return false;
2412 }
2416 {
2417 for (unsigned int d = 0; d < dim; ++d)
2418 if ((p[d] < -tolerance) || (p[d] > 1 + tolerance))
2419 return false;
2420 return true;
2421 }
2424 {
2425 // First make sure that we are in the first quadrant or octant
2426 for (unsigned int d = 0; d < dim; ++d)
2427 if (p[d] < -tolerance)
2428 return false;
2429
2430 // Now we also need to make sure that we are below the diagonal line
2431 // or plane that delineates the simplex. This diagonal is given by
2432 // sum(p[d])<=1, and a diagonal a distance eps away is given by
2433 // sum(p[d])<=1+eps*sqrt(d). (For example, the point at (1,1) is a
2434 // distance of 1/sqrt(2) away from the diagonal. That is, its
2435 // sum satisfies
2436 // sum(p[d]) = 2 <= 1 + (1/sqrt(2)) * sqrt(2)
2437 // in other words, it satisfies the predicate with eps=1/sqrt(2).)
2438 double sum = 0;
2439 for (unsigned int d = 0; d < dim; ++d)
2440 sum += p[d];
2441 return (sum <= 1 + tolerance * std::sqrt(1. * dim));
2442 }
2444 {
2445 // A pyramid only lives in the upper half-space:
2446 if (p[z_coordinate] < -tolerance)
2447 return false;
2448
2449 // It also only lives in the space below z=1:
2450 if (p[z_coordinate] > 1 + tolerance)
2451 return false;
2452
2453 // Within what's left of the space, a pyramid is a cone that tapers
2454 // towards the top. First compute the distance of the point to the
2455 // axis in the max norm (this is the right norm because the vertices
2456 // of the pyramid are at points +/-1, +/-1):
2457 const double distance_from_axis =
2458 std::max(std::fabs(p[x_coordinate]), std::fabs(p[y_coordinate]));
2459
2460 // We are inside the pyramid if the distance from the axis is less
2461 // than (1-z)
2462 return (distance_from_axis <= 1 + tolerance - p[z_coordinate]);
2463 }
2465 {
2466 // The wedge we use is a triangle extruded into the third
2467 // dimension by one unit. So we can use the same logic as for
2468 // triangles above (i.e., for the simplex above, using dim==2)
2469 // and then check the third dimension separately.
2470
2471 if ((p[x_coordinate] < -tolerance) || (p[y_coordinate] < -tolerance))
2472 return false;
2473
2474 const double sum = p[x_coordinate] + p[y_coordinate];
2475 if (sum > 1 + tolerance * std::sqrt(2.0))
2476 return false;
2477
2478 if (p[z_coordinate] < -tolerance)
2479 return false;
2480 if (p[z_coordinate] > 1 + tolerance)
2481 return false;
2482
2483 return true;
2484 }
2485 default:
2487 }
2488
2489 return false;
2490}
2491
2492
2493
2494template <int dim>
2495inline Tensor<1, dim>
2497 const unsigned int i) const
2498{
2500 AssertIndexRange(i, dim - 1);
2501
2502 switch (this->kind)
2503 {
2511 {
2512 AssertIndexRange(face_no, 3);
2513 static const std::array<Tensor<1, dim>, 3> table = {
2514 {Point<dim>(1, 0),
2515 Point<dim>(-std::sqrt(0.5), +std::sqrt(0.5)),
2516 Point<dim>(0, -1)}};
2517
2518 return table[face_no];
2519 }
2521 {
2522 AssertIndexRange(face_no, 4);
2523 static const ndarray<Tensor<1, dim>, 4, 2> table = {
2524 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
2525 {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
2526 {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
2527 {{Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
2528 +std::pow(1.0 / 3.0, 1.0 / 4.0),
2529 0),
2530 Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
2531 0,
2532 +std::pow(1.0 / 3.0, 1.0 / 4.0))}}}};
2533
2534 return table[face_no][i];
2535 }
2537 {
2538 AssertIndexRange(face_no, 5);
2539 static const ndarray<Tensor<1, dim>, 5, 2> table = {
2540 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
2541 {{Point<dim>(+1.0 / sqrt(2.0), 0, +1.0 / sqrt(2.0)),
2542 Point<dim>(0, 1, 0)}},
2543 {{Point<dim>(+1.0 / sqrt(2.0), 0, -1.0 / sqrt(2.0)),
2544 Point<dim>(0, 1, 0)}},
2545 {{Point<dim>(1, 0, 0),
2546 Point<dim>(0, +1.0 / sqrt(2.0), +1.0 / sqrt(2.0))}},
2547 {{Point<dim>(1, 0, 0),
2548 Point<dim>(0, +1.0 / sqrt(2.0), -1.0 / sqrt(2.0))}}}};
2549
2550 return table[face_no][i];
2551 }
2553 {
2554 AssertIndexRange(face_no, 5);
2555 static const ndarray<Tensor<1, dim>, 5, 2> table = {
2556 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
2557 {{Point<dim>(1, 0, 0), Point<dim>(0, 1, 0)}},
2558 {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
2559 {{Point<dim>(-1 / std::sqrt(2.0), +1 / std::sqrt(2.0), 0),
2560 Point<dim>(0, 0, 1)}},
2561 {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}}}};
2562
2563 return table[face_no][i];
2564 }
2565 default:
2567 }
2568
2569
2570 return {};
2571}
2572
2573
2574
2575template <int dim>
2576inline Tensor<1, dim>
2577ReferenceCell::unit_normal_vectors(const unsigned int face_no) const
2578{
2579 AssertDimension(dim, this->get_dimension());
2580
2581 if (is_hyper_cube())
2582 {
2585 }
2586 else if (dim == 2)
2587 {
2589
2590 // Return the rotated vector
2591 return cross_product_2d(unit_tangential_vectors<dim>(face_no, 0));
2592 }
2593 else if (dim == 3)
2594 {
2595 return cross_product_3d(unit_tangential_vectors<dim>(face_no, 0),
2596 unit_tangential_vectors<dim>(face_no, 1));
2597 }
2598
2600
2601 return {};
2602}
2603
2604
2605
2606inline unsigned int
2607ReferenceCell::n_face_orientations(const unsigned int face_no) const
2608{
2609 AssertIndexRange(face_no, n_faces());
2610 if (get_dimension() == 1)
2611 return 1;
2612 if (get_dimension() == 2)
2613 return 2;
2615 return 8;
2616 else if (face_reference_cell(face_no) == ReferenceCells::Triangle)
2617 return 6;
2618
2621}
2622
2623
2624
2625inline bool
2627 const unsigned int line,
2628 const unsigned int face,
2629 const unsigned char combined_face_orientation,
2630 const bool line_orientation) const
2631{
2632 if (*this == ReferenceCells::Hexahedron)
2633 {
2634 static constexpr ::ndarray<bool, 2, 8> bool_table{
2635 {{{true, true, false, true, false, false, true, false}},
2636 {{true, true, true, false, false, false, false, true}}}};
2637
2638 return (line_orientation ==
2639 bool_table[line / 2][combined_face_orientation]);
2640 }
2641 else if (*this == ReferenceCells::Tetrahedron)
2642 {
2643 static constexpr unsigned int X = numbers::invalid_unsigned_int;
2644 static constexpr ::ndarray<unsigned int, 4, 3> combined_lines{
2645 {{{0, 0, 0}}, {{X, 0, 1}}, {{X, 0, X}}, {{X, X, X}}}};
2646
2647 const auto combined_line = combined_lines[face][line];
2648
2649 Assert(combined_line != X,
2650 ExcMessage(
2651 "This function can only be called for following face-line "
2652 "combinations: (0,0), (0,1), (0,2), (1,1), (1,2), (2,1),"));
2653
2654 static constexpr ::ndarray<bool, 2, 6> bool_table{
2655 {{{false, true, false, true, false, true}},
2656 {{true, false, true, false, true, false}}}};
2657
2658 return (line_orientation ==
2659 bool_table[combined_line][combined_face_orientation]);
2660 }
2661 else
2662 // TODO: This might actually be wrong for some of the other
2663 // kinds of objects. We should check this
2664 return true;
2665}
2666
2667
2668
2669namespace internal
2670{
2671 template <typename T>
2673 {
2674 public:
2688
2692 virtual ~NoPermutation() noexcept override = default;
2693
2697 virtual void
2698 print_info(std::ostream &out) const override
2699 {
2700 out << '[';
2701
2702 const unsigned int n_vertices = entity_type.n_vertices();
2703
2704 for (unsigned int i = 0; i < n_vertices; ++i)
2705 {
2706 out << vertices_0[i];
2707 if (i + 1 != n_vertices)
2708 out << ',';
2709 }
2710
2711 out << "] is not a valid permutation of [";
2712
2713 for (unsigned int i = 0; i < n_vertices; ++i)
2714 {
2715 out << vertices_1[i];
2716 if (i + 1 != n_vertices)
2717 out << ',';
2718 }
2719
2720 out << "]." << std::endl;
2721 }
2722
2727
2732
2737 };
2738
2750 << "The reference-cell type used on this cell (" << arg1.to_string()
2751 << ") does not match the reference-cell type of the finite element "
2752 << "associated with this cell (" << arg2.to_string() << "). "
2753 << "Did you accidentally use simplex elements on hypercube meshes "
2754 << "(or the other way around), or are you using a mixed mesh and "
2755 << "assigned a simplex element to a hypercube cell (or the other "
2756 << "way around) via the active_fe_index?");
2757} // namespace internal
2758
2759
2760
2761template <typename T, std::size_t N>
2762inline unsigned char
2763ReferenceCell::compute_orientation(const std::array<T, N> &vertices_0,
2764 const std::array<T, N> &vertices_1) const
2765{
2766 Assert(N >= n_vertices(),
2767 ExcMessage("The number of array elements must be equal to or "
2768 "greater than the number of vertices of the cell "
2769 "referenced by this object."));
2770
2771 // Call the non-deprecated function, taking care of calling it only with
2772 // those array elements that we actually care about (see the note
2773 // in the documentation about the arguments potentially being
2774 // larger arrays than necessary).
2776 make_array_view(vertices_0.begin(), vertices_0.begin() + n_vertices()),
2777 make_array_view(vertices_1.begin(), vertices_1.begin() + n_vertices()));
2778}
2779
2780
2781
2782template <typename T>
2783unsigned char
2785 const ArrayView<const T> &vertices_0,
2786 const ArrayView<const T> &vertices_1) const
2787{
2788 Assert(vertices_0.size() == n_vertices(),
2789 ExcMessage("The number of array elements must be equal to "
2790 "the number of vertices of the cell "
2791 "referenced by this object."));
2792 Assert(vertices_1.size() == n_vertices(),
2793 ExcMessage("The number of array elements must be equal to "
2794 "the number of vertices of the cell "
2795 "referenced by this object."));
2796
2797 auto compute_orientation = [&](const auto &table) {
2798 for (unsigned char o = 0; o < table.size(); ++o)
2799 {
2800 bool match = true;
2801 for (unsigned int j = 0; j < table[o].size(); ++j)
2802 match = (match && vertices_0[j] == vertices_1[table[o][j]]);
2803
2804 if (match)
2805 return o;
2806 }
2807
2808 Assert(false, (internal::NoPermutation<T>(*this, vertices_0, vertices_1)));
2809 return std::numeric_limits<unsigned char>::max();
2810 };
2811
2812 switch (this->kind)
2813 {
2815 // TODO: we can get rid of this special-case and use
2816 // vertex_vertex_permutations once we make 0 the default orientation.
2818 if (vertices_0[0] == vertices_1[0])
2820 break;
2827 default:
2829 }
2830
2831 Assert(false, (internal::NoPermutation<T>(*this, vertices_0, vertices_1)));
2832 return std::numeric_limits<unsigned char>::max();
2833}
2834
2835
2836
2837template <typename T, std::size_t N>
2838inline std::array<T, N>
2840 const std::array<T, N> &vertices,
2841 const unsigned int orientation) const
2842{
2843 Assert(N >= n_vertices(),
2844 ExcMessage("The number of array elements must be equal to or "
2845 "greater than the number of vertices of the cell "
2846 "referenced by this object."));
2847
2848 // Call the non-deprecated function, taking care of calling it only with
2849 // those array elements that we actually care about (see the note
2850 // in the documentation about the arguments potentially being
2851 // larger arrays than necessary).
2852 const auto permutation = permute_by_combined_orientation(
2853 make_array_view(vertices.begin(), vertices.begin() + n_vertices()),
2854 orientation);
2855
2856 std::array<T, N> temp;
2857 std::copy(permutation.begin(), permutation.end(), temp.begin());
2858
2859 return temp;
2860}
2861
2862
2863
2864template <typename T>
2865boost::container::small_vector<T, 8>
2868 const unsigned char orientation) const
2869{
2871 boost::container::small_vector<T, 8> result(n_vertices());
2872
2873 auto permute = [&](const auto &table) {
2874 AssertIndexRange(orientation, table.size());
2875 for (unsigned int j = 0; j < table[orientation].size(); ++j)
2876 result[j] = vertices[table[orientation][j]];
2877 };
2878
2879 switch (this->kind)
2880 {
2882 // TODO: we can get rid of this special-case and use
2883 // vertex_vertex_permutations once we make 0 the default orientation.
2884 std::copy(vertices.begin(), vertices.end(), result.begin());
2885 break;
2887 permute(line_vertex_permutations);
2888 break;
2891 break;
2894 break;
2895 default:
2897 }
2898
2899 return result;
2900}
2901
2902
2903
2904inline unsigned char
2906 const unsigned char orientation) const
2907{
2908 switch (this->kind)
2909 {
2911 // Things are always default-oriented in 1D
2912 return orientation;
2913
2915 // the 1d orientations are the identity and a flip: i.e., the identity
2916 // and an involutory mapping
2917 return orientation;
2918
2920 {
2921 AssertIndexRange(orientation, 6);
2922 constexpr std::array<unsigned char, 6> inverses{{0, 1, 2, 5, 4, 3}};
2923 return inverses[orientation];
2924 }
2926 {
2927 AssertIndexRange(orientation, 8);
2928 constexpr std::array<unsigned char, 8> inverses{
2929 {0, 1, 2, 7, 4, 5, 6, 3}};
2930 return inverses[orientation];
2931 }
2932 default:
2934 }
2935
2936 return std::numeric_limits<unsigned char>::max();
2937}
2938
2939
2940
2941template <>
2942unsigned int
2943ReferenceCell::vtk_lexicographic_to_node_index<0>(
2944 const std::array<unsigned, 0> &node_indices,
2945 const std::array<unsigned, 0> &nodes_per_direction,
2946 const bool legacy_format) const;
2947
2948template <>
2949unsigned int
2950ReferenceCell::vtk_lexicographic_to_node_index<1>(
2951 const std::array<unsigned, 1> &node_indices,
2952 const std::array<unsigned, 1> &nodes_per_direction,
2953 const bool legacy_format) const;
2954
2955template <>
2956unsigned int
2957ReferenceCell::vtk_lexicographic_to_node_index<2>(
2958 const std::array<unsigned, 2> &node_indices,
2959 const std::array<unsigned, 2> &nodes_per_direction,
2960 const bool legacy_format) const;
2961
2962template <>
2963unsigned int
2964ReferenceCell::vtk_lexicographic_to_node_index<3>(
2965 const std::array<unsigned, 3> &node_indices,
2966 const std::array<unsigned, 3> &nodes_per_direction,
2967 const bool legacy_format) const;
2968
2970
2971#endif
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
std::size_t size() const
Definition array_view.h:684
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)
ArrayView< const unsigned int > faces_for_given_vertex(const unsigned int vertex_index) const
unsigned int child_cell_on_face(const unsigned int face, const unsigned int subface, const unsigned char face_orientation=default_combined_face_orientation()) const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices() const
unsigned int n_vertices() const
boost::container::small_vector< T, 8 > permute_by_combined_orientation(const ArrayView< const T > &vertices, const unsigned char orientation) const
Quadrature< dim > get_gauss_type_quadrature(const unsigned n_points_1d) const
bool is_hyper_cube() const
unsigned int vtk_linear_type() const
void serialize(Archive &archive, const unsigned int)
unsigned char get_inverse_combined_orientation(const unsigned char orientation) const
unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const unsigned char face_orientation) const
static constexpr unsigned char default_combined_face_orientation()
unsigned char compute_orientation(const std::array< T, N > &vertices_0, const std::array< T, N > &vertices_1) const
bool standard_vs_true_line_orientation(const unsigned int line, const unsigned int face, const unsigned char face_orientation, const bool line_orientation) const
std::array< unsigned int, 2 > standard_vertex_to_face_and_vertex_index(const unsigned int vertex) 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 unsigned char face_orientation) const
const Quadrature< dim > & get_nodal_type_quadrature() const
Point< dim > face_vertex_location(const unsigned int face, const unsigned int vertex) const
static constexpr ndarray< unsigned int, 1, 1 > vertex_vertex_permutations
std_cxx20::ranges::iota_view< unsigned int, unsigned int > isotropic_child_indices() const
std::array< unsigned int, 2 > standard_line_to_face_and_line_index(const unsigned int line) const
double d_linear_shape_function(const Point< dim > &xi, const unsigned int i) const
Tensor< 1, dim > unit_tangential_vectors(const unsigned int face_no, const unsigned int i) const
static constexpr ndarray< unsigned int, 2, 2 > line_vertex_permutations
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
std::array< T, N > permute_according_orientation(const std::array< T, N > &vertices, const unsigned int orientation) const
bool contains_point(const Point< dim > &p, const double tolerance=0) const
unsigned int n_isotropic_children() const
unsigned int gmsh_element_type() const
unsigned int n_face_orientations(const unsigned int face_no) const
static constexpr ndarray< unsigned int, 8, 4 > quadrilateral_vertex_permutations
double volume() const
unsigned int n_faces() const
unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex) const
std::uint8_t kind
constexpr ReferenceCell()
unsigned int exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const
unsigned int unv_vertex_to_deal_vertex(const unsigned int vertex_n) const
unsigned int vtk_vertex_to_deal_vertex(const unsigned int vertex_index) const
std::unique_ptr< Mapping< dim, spacedim > > get_default_mapping(const unsigned int degree) const
Quadrature< dim > get_midpoint_quadrature() const
unsigned int vtk_quadratic_type() const
static constexpr unsigned char reversed_combined_line_orientation()
unsigned int n_lines() const
unsigned int vtk_lagrange_type() const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > line_indices() const
unsigned int get_dimension() const
ReferenceCell face_reference_cell(const unsigned int face_no) const
unsigned int exodusii_face_to_deal_face(const unsigned int face_n) const
static ReferenceCell n_vertices_to_type(const int dim, const unsigned int n_vertices)
unsigned int standard_to_real_face_line(const unsigned int line, const unsigned int face, const unsigned char face_orientation) const
static constexpr std::size_t memory_consumption()
friend std::istream & operator>>(std::istream &in, ReferenceCell &reference_cell)
const Mapping< dim, spacedim > & get_default_linear_mapping() const
std::string to_string() const
bool is_simplex() const
unsigned int face_to_cell_lines(const unsigned int face, const unsigned int line, const unsigned char face_orientation) const
constexpr bool operator!=(const ReferenceCell &type) const
unsigned char get_combined_orientation(const ArrayView< const T > &vertices_0, const ArrayView< const T > &vertices_1) const
constexpr bool operator==(const ReferenceCell &type) const
Point< dim > barycenter() const
Tensor< 1, dim > unit_normal_vectors(const unsigned int face_no) const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices() const
friend std::ostream & operator<<(std::ostream &out, const ReferenceCell &reference_cell)
Point< dim > closest_point(const Point< dim > &p) const
static constexpr ndarray< unsigned int, 6, 3 > triangle_vertex_permutations
Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i) const
virtual void print_info(std::ostream &out) const override
const ReferenceCell entity_type
virtual ~NoPermutation() noexcept override=default
const ArrayView< const T > vertices_0
const ArrayView< const T > vertices_1
NoPermutation(const ReferenceCell &entity_type, const ArrayView< const T > &vertices_0, const ArrayView< const T > &vertices_1)
#define DEAL_II_DEPRECATED
Definition config.h:207
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 3 > vertices[4]
static ::ExceptionBase & ExcNonMatchingReferenceCellTypes(ReferenceCell arg1, ReferenceCell arg2)
#define Assert(cond, exc)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:539
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Quadrilateral
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Pyramid
constexpr const ReferenceCell Invalid
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Hexahedron
constexpr const ReferenceCell & get_hypercube()
constexpr const ReferenceCell Vertex
constexpr const ReferenceCell Line
constexpr const ReferenceCell & get_simplex()
constexpr ReferenceCell make_reference_cell_from_int(const std::uint8_t kind)
std::tuple< bool, bool, bool > split_face_orientation(const unsigned char combined_face_orientation)
static const unsigned int invalid_unsigned_int
Definition types.h:220
boost::integer_range< IncrementableType > iota_view
Definition iota_view.h:45
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 > pow(const ::VectorizedArray< Number, width > &, const Number p)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:107
std::istream & operator>>(std::istream &in, ReferenceCell &reference_cell)
std::ostream & operator<<(std::ostream &out, const ReferenceCell &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 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 std::array< unsigned int, 2 > standard_hex_line_to_quad_line_index(const unsigned int line)
static unsigned int standard_to_real_face_line(const unsigned int line, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
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 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 Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)