Reference documentation for deal.II version GIT relicensing-362-gab68047079 2024-04-11 17:50: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
876 {{{1, 0}}, {{0, 1}}}};
877
878
883 {{{0, 2, 1}},
884 {{0, 1, 2}},
885 {{2, 1, 0}},
886 {{2, 0, 1}},
887 {{1, 0, 2}},
888 {{1, 2, 0}}}};
889
893 static constexpr ndarray<unsigned int, 8, 4>
895 {{0, 1, 2, 3}},
896 {{2, 3, 0, 1}},
897 {{2, 0, 3, 1}},
898 {{3, 1, 2, 0}},
899 {{3, 2, 1, 0}},
900 {{1, 0, 3, 2}},
901 {{1, 3, 0, 2}}}};
902
907 friend constexpr ReferenceCell
909
910 friend std::ostream &
911 operator<<(std::ostream &out, const ReferenceCell &reference_cell);
912
913 friend std::istream &
914 operator>>(std::istream &in, ReferenceCell &reference_cell);
915};
916
917
925std::ostream &
926operator<<(std::ostream &out, const ReferenceCell &reference_cell);
927
935std::istream &
936operator>>(std::istream &in, ReferenceCell &reference_cell);
937
938
939inline constexpr ReferenceCell::ReferenceCell(const std::uint8_t kind)
940 : kind(kind)
941{}
942
943
944
945inline constexpr ReferenceCell::operator std::uint8_t() const
946{
947 return kind;
948}
949
950
951
952inline constexpr bool
954{
955 return kind == type.kind;
956}
957
958
959
960inline constexpr bool
962{
963 return kind != type.kind;
964}
965
966
967
968namespace internal
969{
970 inline constexpr ReferenceCell
971 make_reference_cell_from_int(const std::uint8_t kind)
972 {
973#ifndef DEAL_II_CXX14_CONSTEXPR_BUG
974 // Make sure these are the only indices from which objects can be
975 // created.
976 Assert((kind == std::numeric_limits<std::uint8_t>::max()) || (kind < 8),
978#endif
979
980 // Call the private constructor, which we can from here because this
981 // function is a 'friend'.
982 return {kind};
983 }
984} // namespace internal
985
986
987
996{
997 constexpr const ReferenceCell Vertex =
999 constexpr const ReferenceCell Line =
1001 constexpr const ReferenceCell Triangle =
1005 constexpr const ReferenceCell Tetrahedron =
1007 constexpr const ReferenceCell Pyramid =
1009 constexpr const ReferenceCell Wedge =
1011 constexpr const ReferenceCell Hexahedron =
1013 constexpr const ReferenceCell Invalid =
1015 std::numeric_limits<std::uint8_t>::max());
1016
1022 template <int dim>
1023 constexpr const ReferenceCell &
1024 get_simplex();
1025
1031 template <int dim>
1032 constexpr const ReferenceCell &
1033 get_hypercube();
1034} // namespace ReferenceCells
1035
1036
1037
1039 : ReferenceCell(ReferenceCells::Invalid)
1040{}
1041
1042
1043
1044template <class Archive>
1045inline void
1046ReferenceCell::serialize(Archive &archive, const unsigned int /*version*/)
1047{
1048 archive &kind;
1049}
1050
1051
1052
1053inline constexpr std::size_t
1058
1059
1060
1062ReferenceCell::faces_for_given_vertex(const unsigned int vertex) const
1063{
1065 switch (this->kind)
1066 {
1068 return {&GeometryInfo<2>::vertex_to_face[vertex][0], 1};
1070 return {&GeometryInfo<2>::vertex_to_face[vertex][0], 2};
1072 {
1073 static constexpr ndarray<unsigned int, 3, 2> table = {
1074 {{{0, 2}}, {{0, 1}}, {{1, 2}}}};
1075 return table[vertex];
1076 }
1078 {
1079 static constexpr ndarray<unsigned int, 4, 3> table = {
1080 {{{0, 1, 2}}, {{0, 1, 3}}, {{0, 2, 3}}, {{1, 2, 3}}}};
1081
1082 return table[vertex];
1083 }
1085 {
1086 static constexpr unsigned int X = numbers::invalid_unsigned_int;
1087 static constexpr ndarray<unsigned int, 5, 4> table = {
1088 {{{0, 1, 3, X}},
1089 {{0, 2, 3, X}},
1090 {{0, 1, 4, X}},
1091 {{0, 2, 4, X}},
1092 {{1, 2, 3, 4}}}};
1093
1094 return {&table[vertex][0], vertex == 4 ? 4u : 3u};
1095 }
1097 {
1099 static constexpr ndarray<unsigned int, 6, 3> table = {{{{0, 2, 4}},
1100 {{0, 2, 3}},
1101 {{0, 3, 4}},
1102 {{1, 2, 4}},
1103 {{1, 2, 3}},
1104 {{1, 3, 4}}}};
1105
1106 return table[vertex];
1107 }
1109 return {&GeometryInfo<3>::vertex_to_face[vertex][0], 3};
1110 default:
1112 }
1113
1114 return {};
1115}
1116
1117
1118
1119inline bool
1121{
1122 return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
1125}
1126
1127
1128
1129inline bool
1131{
1132 return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
1133 *this == ReferenceCells::Triangle ||
1135}
1136
1137
1138
1139inline unsigned int
1141{
1142 switch (this->kind)
1143 {
1145 return 0;
1147 return 1;
1150 return 2;
1155 return 3;
1156 default:
1158 }
1159
1161}
1162
1163
1164
1165template <int dim>
1168{
1169 return Quadrature<dim>(std::vector<Point<dim>>({barycenter<dim>()}),
1170 std::vector<double>({volume()}));
1171}
1172
1173
1174
1175inline unsigned int
1177{
1178 switch (this->kind)
1179 {
1181 return 1;
1183 return 2;
1185 return 3;
1187 return 4;
1189 return 4;
1191 return 5;
1193 return 6;
1195 return 8;
1196 default:
1198 }
1199
1201}
1202
1203
1204
1205inline unsigned int
1207{
1208 switch (this->kind)
1209 {
1211 return 0;
1213 return 1;
1215 return 3;
1217 return 4;
1219 return 6;
1221 return 8;
1223 return 9;
1225 return 12;
1226 default:
1228 }
1229
1231}
1232
1233
1234
1235template <int dim>
1237ReferenceCell::vertex(const unsigned int v) const
1238{
1241
1242 switch (dim)
1243 {
1244 case 0:
1245 {
1246 if (*this == ReferenceCells::Vertex)
1247 return Point<dim>();
1248 break;
1249 }
1250 case 1:
1251 {
1252 static const Point<dim> vertices[2] = {
1253 Point<dim>(), // the origin
1254 Point<dim>::unit_vector(0) // unit point along x-axis
1255 };
1256 if (*this == ReferenceCells::Line)
1257 return vertices[v];
1258 break;
1259 }
1260 case 2:
1261 {
1262 switch (this->kind)
1263 {
1265 {
1266 static const Point<dim> vertices[3] = {
1267 Point<dim>(), // the origin
1268 Point<dim>::unit_vector(0), // unit point along x-axis
1269 Point<dim>::unit_vector(1) // unit point along y-axis
1270 };
1271 return vertices[v];
1272 }
1274 {
1275 static const Point<dim> vertices[4] = {
1276 // First the two points on the x-axis
1277 Point<dim>(),
1279 // Then these two points shifted in the y-direction
1282 return vertices[v];
1283 }
1284 }
1285 break;
1286 }
1287 case 3:
1288 {
1289 switch (this->kind)
1290 {
1292 {
1293 static const Point<dim> vertices[4] = {
1294 Point<dim>(), // the origin
1295 Point<dim>::unit_vector(0), // unit point along x-axis
1296 Point<dim>::unit_vector(1), // unit point along y-axis
1297 Point<dim>::unit_vector(2) // unit point along z-axis
1298 };
1299 return vertices[v];
1300 }
1302 {
1303 static const Point<dim> vertices[5] = {
1304 Point<dim>{-1.0, -1.0, 0.0},
1305 Point<dim>{+1.0, -1.0, 0.0},
1306 Point<dim>{-1.0, +1.0, 0.0},
1307 Point<dim>{+1.0, +1.0, 0.0},
1308 Point<dim>{+0.0, +0.0, 1.0}};
1309 return vertices[v];
1310 }
1312 {
1313 static const Point<dim> vertices[6] = {
1314 // First the three points on the triangular base of the
1315 // wedge:
1316 Point<dim>(),
1319 // And now everything shifted in the z-direction again
1323 return vertices[v];
1324 }
1326 {
1327 static const Point<dim> vertices[8] = {
1328 // First the two points on the x-axis
1329 Point<dim>(),
1331 // Then these two points shifted in the y-direction
1334 // And now all four points shifted in the z-direction
1341 return vertices[v];
1342 }
1343 }
1344 break;
1345 }
1346 default:
1348 }
1349
1351 return Point<dim>();
1352}
1353
1354
1355inline unsigned int
1357{
1358 switch (this->kind)
1359 {
1361 return 0;
1363 return 2;
1365 return 3;
1367 return 4;
1369 return 4;
1371 return 5;
1373 return 5;
1375 return 6;
1376 default:
1378 }
1379
1381}
1382
1383
1384
1387{
1388 return {0U, n_faces()};
1389}
1390
1391
1392
1393inline unsigned int
1395{
1396 switch (this->kind)
1397 {
1399 return 0;
1401 return 2;
1403 return 4;
1405 return 4;
1407 return 8;
1409 // We haven't yet decided how to refine pyramids. Update this when we
1410 // have
1414 return 8;
1416 return 8;
1417 default:
1419 }
1420
1422}
1423
1424
1425
1428{
1429 return {0U, n_isotropic_children()};
1430}
1431
1432
1433
1436{
1437 return {0U, n_vertices()};
1438}
1439
1440
1441
1444{
1445 return {0U, n_lines()};
1446}
1447
1448
1449
1450inline ReferenceCell
1451ReferenceCell::face_reference_cell(const unsigned int face_no) const
1452{
1453 AssertIndexRange(face_no, n_faces());
1454
1455 switch (this->kind)
1456 {
1463 return ReferenceCells::Line;
1467 if (face_no == 0)
1469 else
1472 if (face_no > 1)
1474 else
1478 default:
1480 }
1481
1483}
1484
1485
1486
1487inline constexpr unsigned char
1489{
1490 // Our convention is that 'orientation' has a default value of true and
1491 // occupies the least-significant bit while rotate and flip have default
1492 // values of 'false' and occupy the second and third bits.
1493 return 0b001;
1494}
1495
1496
1497
1498inline constexpr unsigned char
1500{
1501 // For a reversed line 'orientation' is false and neither flip nor rotate are
1502 // defined.
1503 return 0b000;
1504}
1505
1506
1507
1508inline unsigned int
1510 const unsigned int face,
1511 const unsigned int subface,
1512 const unsigned char combined_face_orientation) const
1513{
1514 AssertIndexRange(face, n_faces());
1516
1517 switch (this->kind)
1518 {
1521 {
1523 break;
1524 }
1526 {
1527 static constexpr ndarray<unsigned int, 3, 2> subcells = {
1528 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
1529
1530 return subcells[face][subface];
1531 }
1533 {
1534 const auto [face_orientation, face_rotation, face_flip] =
1535 internal::split_face_orientation(combined_face_orientation);
1536
1539 face,
1540 subface,
1541 face_orientation,
1542 face_flip,
1543 face_rotation);
1544 }
1548 {
1550 break;
1551 }
1553 {
1554 const auto [face_orientation, face_rotation, face_flip] =
1555 internal::split_face_orientation(combined_face_orientation);
1556
1559 face,
1560 subface,
1561 face_orientation,
1562 face_flip,
1563 face_rotation);
1564 }
1565 default:
1567 }
1568
1570}
1571
1572
1573
1574inline std::array<unsigned int, 2>
1576 const unsigned int vertex) const
1577{
1579 // Work around a GCC warning at higher optimization levels by making all of
1580 // these tables the same size
1581 constexpr unsigned int X = numbers::invalid_unsigned_int;
1582
1583 switch (this->kind)
1584 {
1588 break;
1590 {
1591 static constexpr ndarray<unsigned int, 6, 2> table = {
1592 {{{0, 0}}, {{0, 1}}, {{1, 1}}, {{X, X}}, {{X, X}}, {{X, X}}}};
1593
1594 return table[vertex];
1595 }
1597 {
1599 vertex);
1600 }
1602 {
1603 static constexpr ndarray<unsigned int, 6, 2> table = {
1604 {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 2}}, {{X, X}}, {{X, X}}}};
1605
1606 return table[vertex];
1607 }
1609 {
1610 static constexpr ndarray<unsigned int, 6, 2> table = {
1611 {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{0, 3}}, {{1, 2}}, {{X, X}}}};
1612
1613 return table[vertex];
1614 }
1616 {
1617 static constexpr ndarray<unsigned int, 6, 2> table = {
1618 {{{0, 1}}, {{0, 0}}, {{0, 2}}, {{1, 0}}, {{1, 1}}, {{1, 2}}}};
1619
1620 return table[vertex];
1621 }
1623 {
1625 vertex);
1626 }
1627 default:
1629 }
1630
1631 return {};
1632}
1633
1634
1635
1636inline std::array<unsigned int, 2>
1638 const unsigned int line) const
1639{
1640 AssertIndexRange(line, n_lines());
1641
1642 switch (this->kind)
1643 {
1648 {
1650 break;
1651 }
1653 {
1654 static const std::array<unsigned int, 2> table[6] = {
1655 {{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}};
1656
1657 return table[line];
1658 }
1660 {
1661 static const std::array<unsigned int, 2> table[8] = {{{0, 0}},
1662 {{0, 1}},
1663 {{0, 2}},
1664 {{0, 3}},
1665 {{1, 2}},
1666 {{2, 1}},
1667 {{1, 1}},
1668 {{2, 2}}};
1669
1670 return table[line];
1671 }
1673 {
1674 static const std::array<unsigned int, 2> table[9] = {{{0, 0}},
1675 {{0, 2}},
1676 {{0, 1}},
1677 {{1, 0}},
1678 {{1, 1}},
1679 {{1, 2}},
1680 {{2, 0}},
1681 {{2, 1}},
1682 {{3, 1}}};
1683
1684 return table[line];
1685 }
1687 {
1689 }
1690 default:
1692 }
1693
1694 return {};
1695}
1696
1697
1698inline unsigned int
1700 const unsigned int vertex) const
1701{
1703 AssertIndexRange(line, n_lines());
1704
1705 switch (this->kind)
1706 {
1709 return vertex;
1711 {
1712 static constexpr ndarray<unsigned int, 3, 2> table = {
1713 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
1714 return table[line][vertex];
1715 }
1717 {
1718 static constexpr ndarray<unsigned int, 4, 2> table = {
1719 {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
1720 return table[line][vertex];
1721 }
1723 {
1724 static constexpr ndarray<unsigned int, 6, 2> table = {
1725 {{{0, 1}}, {{1, 2}}, {{2, 0}}, {{0, 3}}, {{1, 3}}, {{2, 3}}}};
1726 return table[line][vertex];
1727 }
1729 {
1730 static constexpr ndarray<unsigned int, 8, 2> table = {{{{0, 2}},
1731 {{1, 3}},
1732 {{0, 1}},
1733 {{2, 3}},
1734 {{4, 0}},
1735 {{1, 4}},
1736 {{2, 4}},
1737 {{4, 3}}}};
1738 return table[line][vertex];
1739 }
1741 {
1742 static constexpr ndarray<unsigned int, 9, 2> table = {{{{1, 0}},
1743 {{2, 1}},
1744 {{0, 2}},
1745 {{3, 4}},
1746 {{4, 5}},
1747 {{5, 3}},
1748 {{0, 3}},
1749 {{1, 4}},
1750 {{2, 5}}}};
1751 return table[line][vertex];
1752 }
1754 {
1755 // first four lines comprise the bottom face, next four are the top,
1756 // and the last four are 'bottom to top'
1757 static constexpr ndarray<unsigned int, 12, 2> table = {{{{0, 2}},
1758 {{1, 3}},
1759 {{0, 1}},
1760 {{2, 3}},
1761 {{4, 6}},
1762 {{5, 7}},
1763 {{4, 5}},
1764 {{6, 7}},
1765 {{0, 4}},
1766 {{1, 5}},
1767 {{2, 6}},
1768 {{3, 7}}}};
1769 return table[line][vertex];
1770 }
1771
1772 default:
1774 }
1775
1777}
1778
1779
1780inline unsigned int
1782 const unsigned int face,
1783 const unsigned int line,
1784 const unsigned char combined_face_orientation) const
1785{
1786 AssertIndexRange(face, n_faces());
1788
1789 static constexpr unsigned int X = numbers::invalid_unsigned_int;
1790
1791 switch (this->kind)
1792 {
1794 {
1796 break;
1797 }
1799 {
1800 const auto [face_orientation, face_rotation, face_flip] =
1801 internal::split_face_orientation(combined_face_orientation);
1802
1804 face, line, face_orientation, face_flip, face_rotation);
1805 }
1807 {
1808 return face;
1809 }
1811 {
1812 const auto [face_orientation, face_rotation, face_flip] =
1813 internal::split_face_orientation(combined_face_orientation);
1814
1816 face, line, face_orientation, face_flip, face_rotation);
1817 }
1819 {
1820 static constexpr ndarray<unsigned int, 4, 3> table = {
1821 {{{0, 1, 2}}, {{0, 3, 4}}, {{2, 5, 3}}, {{1, 4, 5}}}};
1822
1823 return table[face][standard_to_real_face_line(
1824 line, face, combined_face_orientation)];
1825 }
1827 {
1828 static constexpr ndarray<unsigned int, 5, 4> table = {
1829 {{{0, 1, 2, 3}},
1830 {{0, 6, 4, X}},
1831 {{1, 5, 7, X}},
1832 {{2, 4, 5, X}},
1833 {{3, 7, 6, X}}}};
1834
1835 return table[face][standard_to_real_face_line(
1836 line, face, combined_face_orientation)];
1837 }
1839 {
1840 static constexpr ndarray<unsigned int, 5, 4> table = {
1841 {{{0, 2, 1, X}},
1842 {{3, 4, 5, X}},
1843 {{6, 7, 0, 3}},
1844 {{7, 8, 1, 4}},
1845 {{8, 6, 5, 2}}}};
1846
1847 return table[face][standard_to_real_face_line(
1848 line, face, combined_face_orientation)];
1849 }
1851 {
1852 const auto [face_orientation, face_rotation, face_flip] =
1853 internal::split_face_orientation(combined_face_orientation);
1854
1856 face, line, face_orientation, face_flip, face_rotation);
1857 }
1858 default:
1860 }
1861
1863}
1864
1865
1866
1867inline unsigned int
1869 const unsigned int face,
1870 const unsigned int vertex,
1871 const unsigned char combined_face_orientation) const
1872{
1873 AssertIndexRange(face, n_faces());
1875 // TODO: once the default orientation is switched to 0 then we can remove this
1876 // special case for 1D.
1877 if (get_dimension() == 1)
1878 Assert(combined_face_orientation ==
1880 ExcMessage("In 1D, all faces must have the default orientation."));
1881 else
1882 AssertIndexRange(combined_face_orientation, n_face_orientations(face));
1883
1884 switch (this->kind)
1885 {
1887 {
1889 break;
1890 }
1892 {
1893 const auto [face_orientation, face_rotation, face_flip] =
1894 internal::split_face_orientation(combined_face_orientation);
1895
1897 face, vertex, face_orientation, face_flip, face_rotation);
1898 }
1900 {
1901 static constexpr ndarray<unsigned int, 3, 2> table = {
1902 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
1903
1904 return table[face][combined_face_orientation !=
1906 vertex :
1907 (1 - vertex)];
1908 }
1910 {
1911 const auto [face_orientation, face_rotation, face_flip] =
1912 internal::split_face_orientation(combined_face_orientation);
1913
1915 face, vertex, face_orientation, face_flip, face_rotation);
1916 }
1918 {
1919 static constexpr ndarray<unsigned int, 4, 3> table = {
1920 {{{0, 1, 2}}, {{1, 0, 3}}, {{0, 2, 3}}, {{2, 1, 3}}}};
1921
1922 return table[face][standard_to_real_face_vertex(
1923 vertex, face, combined_face_orientation)];
1924 }
1926 {
1927 constexpr auto X = numbers::invalid_unsigned_int;
1928 static constexpr ndarray<unsigned int, 5, 4> table = {
1929 {{{0, 1, 2, 3}},
1930 {{0, 2, 4, X}},
1931 {{3, 1, 4, X}},
1932 {{1, 0, 4, X}},
1933 {{2, 3, 4, X}}}};
1934
1935 return table[face][standard_to_real_face_vertex(
1936 vertex, face, combined_face_orientation)];
1937 }
1939 {
1940 constexpr auto X = numbers::invalid_unsigned_int;
1941 static constexpr ndarray<unsigned int, 6, 4> table = {
1942 {{{1, 0, 2, X}},
1943 {{3, 4, 5, X}},
1944 {{0, 1, 3, 4}},
1945 {{1, 2, 4, 5}},
1946 {{2, 0, 5, 3}}}};
1947
1948 return table[face][standard_to_real_face_vertex(
1949 vertex, face, combined_face_orientation)];
1950 }
1952 {
1953 const auto [face_orientation, face_rotation, face_flip] =
1954 internal::split_face_orientation(combined_face_orientation);
1955
1957 face, vertex, face_orientation, face_flip, face_rotation);
1958 }
1959 default:
1961 }
1962
1964}
1965
1966
1967
1968template <int dim>
1971 const unsigned int vertex) const
1972{
1973 return this->template vertex<dim>(
1975}
1976
1977
1978
1979inline unsigned int
1981 const unsigned int vertex,
1982 const unsigned int face,
1983 const unsigned char face_orientation) const
1984{
1985 AssertIndexRange(face, n_faces());
1987
1988 switch (this->kind)
1989 {
1992 break;
1994 Assert(face_orientation == default_combined_face_orientation(),
1995 ExcMessage(
1996 "In 1D, all faces must have the default orientation."));
1997 return vertex;
2000 return line_vertex_permutations[face_orientation][vertex];
2002 return triangle_vertex_permutations[face_orientation][vertex];
2004 // face 0 is a quadrilateral
2005 if (face == 0)
2006 return quadrilateral_vertex_permutations[face_orientation][vertex];
2007 else
2008 return triangle_vertex_permutations[face_orientation][vertex];
2010 // faces 0 and 1 are triangles
2011 if (face > 1)
2012 return quadrilateral_vertex_permutations[face_orientation][vertex];
2013 else
2014 return triangle_vertex_permutations[face_orientation][vertex];
2016 return quadrilateral_vertex_permutations[face_orientation][vertex];
2017 default:
2019 }
2020
2023}
2024
2025
2026
2027inline unsigned int
2029 const unsigned int line,
2030 const unsigned int face,
2031 const unsigned char combined_face_orientation) const
2032{
2033 AssertIndexRange(face, n_faces());
2035
2036 static constexpr ndarray<unsigned int, 6, 3> triangle_table = {{{{2, 1, 0}},
2037 {{0, 1, 2}},
2038 {{1, 0, 2}},
2039 {{2, 0, 1}},
2040 {{0, 2, 1}},
2041 {{1, 2, 0}}}};
2042
2043
2044 switch (this->kind)
2045 {
2051 break;
2053 return triangle_table[combined_face_orientation][line];
2055 if (face == 0) // The quadrilateral face
2056 {
2057 const auto [face_orientation, face_rotation, face_flip] =
2058 internal::split_face_orientation(combined_face_orientation);
2059
2061 face_orientation,
2062 face_flip,
2063 face_rotation);
2064 }
2065 else // One of the triangular faces
2066 {
2067 return triangle_table[combined_face_orientation][line];
2068 }
2070 if (face > 1) // One of the quadrilateral faces
2071 {
2072 const auto [face_orientation, face_rotation, face_flip] =
2073 internal::split_face_orientation(combined_face_orientation);
2074
2076 face_orientation,
2077 face_flip,
2078 face_rotation);
2079 }
2080 else // One of the triangular faces
2081 return triangle_table[combined_face_orientation][line];
2083 {
2084 static constexpr ndarray<unsigned int, 8, 4> table = {
2085 {{{2, 3, 0, 1}},
2086 {{0, 1, 2, 3}},
2087 {{0, 1, 3, 2}},
2088 {{3, 2, 0, 1}},
2089 {{3, 2, 1, 0}},
2090 {{1, 0, 3, 2}},
2091 {{1, 0, 2, 3}},
2092 {{2, 3, 1, 0}}}};
2093 return table[combined_face_orientation][line];
2094 }
2095 default:
2097 }
2098
2100}
2101
2102
2103
2104namespace ReferenceCells
2105{
2106 template <int dim>
2107 inline constexpr const ReferenceCell &
2109 {
2110 switch (dim)
2111 {
2112 case 0:
2114 case 1:
2115 return ReferenceCells::Line;
2116 case 2:
2118 case 3:
2120 default:
2122 }
2124 }
2125
2126
2127
2128 template <int dim>
2129 inline constexpr const ReferenceCell &
2131 {
2132 switch (dim)
2133 {
2134 case 0:
2136 case 1:
2137 return ReferenceCells::Line;
2138 case 2:
2140 case 3:
2142 default:
2144 }
2146 }
2147} // namespace ReferenceCells
2148
2149
2150inline ReferenceCell
2151ReferenceCell::n_vertices_to_type(const int dim, const unsigned int n_vertices)
2152{
2153 AssertIndexRange(dim, 4);
2155
2156 const auto X = ReferenceCells::Invalid;
2157 static const ndarray<ReferenceCell, 4, 9> table = {
2158 {// dim 0
2159 {{X, ReferenceCells::Vertex, X, X, X, X, X, X, X}},
2160 // dim 1
2161 {{X, X, ReferenceCells::Line, X, X, X, X, X, X}},
2162 // dim 2
2163 {{X,
2164 X,
2165 X,
2168 X,
2169 X,
2170 X,
2171 X}},
2172 // dim 3
2173 {{X,
2174 X,
2175 X,
2176 X,
2180 X,
2183 ExcMessage("The combination of dim = " + std::to_string(dim) +
2184 " and n_vertices = " + std::to_string(n_vertices) +
2185 " does not correspond to a known reference cell type."));
2186 return table[dim][n_vertices];
2187}
2188
2189
2190
2191template <int dim>
2192inline double
2194 const unsigned int i) const
2195{
2198 switch (this->kind)
2199 {
2205 // see also BarycentricPolynomials<2>::compute_value
2207 {
2208 switch (i)
2209 {
2210 case 0:
2211 return 1.0 - xi[std::min(0, dim - 1)] -
2212 xi[std::min(1, dim - 1)];
2213 case 1:
2214 return xi[std::min(0, dim - 1)];
2215 case 2:
2216 return xi[std::min(1, dim - 1)];
2217 default:
2219 }
2220 }
2221 // see also BarycentricPolynomials<3>::compute_value
2223 {
2224 switch (i)
2225 {
2226 case 0:
2227 return 1.0 - xi[std::min(0, dim - 1)] -
2228 xi[std::min(1, dim - 1)] - xi[std::min(2, dim - 1)];
2229 case 1:
2230 return xi[std::min(0, dim - 1)];
2231 case 2:
2232 return xi[std::min(1, dim - 1)];
2233 case 3:
2234 return xi[std::min(2, dim - 1)];
2235 default:
2237 }
2238 }
2239 // see also ScalarLagrangePolynomialPyramid::compute_value()
2241 {
2242 const double Q14 = 0.25;
2243
2244 const double r = xi[std::min(0, dim - 1)];
2245 const double s = xi[std::min(1, dim - 1)];
2246 const double t = xi[std::min(2, dim - 1)];
2247
2248 const double ratio =
2249 (std::fabs(t - 1.0) > 1.0e-14 ? (r * s * t) / (1.0 - t) : 0.0);
2250
2251 if (i == 0)
2252 return Q14 * ((1.0 - r) * (1.0 - s) - t + ratio);
2253 if (i == 1)
2254 return Q14 * ((1.0 + r) * (1.0 - s) - t - ratio);
2255 if (i == 2)
2256 return Q14 * ((1.0 - r) * (1.0 + s) - t - ratio);
2257 if (i == 3)
2258 return Q14 * ((1.0 + r) * (1.0 + s) - t + ratio);
2259 else
2260 return t;
2261 }
2262 // see also ScalarLagrangePolynomialWedge::compute_value()
2265 .d_linear_shape_function<2>(Point<2>(xi[std::min(0, dim - 1)],
2266 xi[std::min(1, dim - 1)]),
2267 i % 3) *
2269 .d_linear_shape_function<1>(Point<1>(xi[std::min(2, dim - 1)]),
2270 i / 3);
2271 default:
2273 }
2274
2275 return 0.0;
2276}
2277
2278
2279
2280template <int dim>
2281inline Tensor<1, dim>
2283 const unsigned int i) const
2284{
2286 switch (this->kind)
2287 {
2293 // see also BarycentricPolynomials<2>::compute_grad()
2295 switch (i)
2296 {
2297 case 0:
2298 return Point<dim>(-1.0, -1.0);
2299 case 1:
2300 return Point<dim>(+1.0, +0.0);
2301 case 2:
2302 return Point<dim>(+0.0, +1.0);
2303 default:
2305 }
2306 default:
2308 }
2309
2310 return Point<dim>(+0.0, +0.0, +0.0);
2311}
2312
2313
2314
2315inline double
2317{
2318 switch (this->kind)
2319 {
2321 return 0;
2323 return 1;
2325 return 1. / 2.;
2327 return 1;
2329 return 1. / 6.;
2331 return 4. / 3.;
2333 return 1. / 2.;
2335 return 1;
2336 default:
2338 }
2339
2340 return 0.0;
2341}
2342
2343
2344
2345template <int dim>
2346inline Point<dim>
2348{
2350
2351 switch (this->kind)
2352 {
2354 return Point<dim>();
2356 return Point<dim>(1. / 2.);
2358 return Point<dim>(1. / 3., 1. / 3.);
2360 return Point<dim>(1. / 2., 1. / 2.);
2362 return Point<dim>(1. / 4., 1. / 4., 1. / 4.);
2364 return Point<dim>(0, 0, 1. / 4.);
2366 return Point<dim>(1. / 3, 1. / 3, 1. / 2.);
2368 return Point<dim>(1. / 2., 1. / 2., 1. / 2.);
2369 default:
2371 }
2372
2373 return Point<dim>();
2374}
2375
2376
2377
2378template <int dim>
2379inline bool
2380ReferenceCell::contains_point(const Point<dim> &p, const double tolerance) const
2381{
2383
2384 // Introduce abbreviations to silence compiler warnings about invalid
2385 // array accesses (that can't happen, of course, but the compiler
2386 // doesn't know that).
2387 constexpr unsigned int x_coordinate = 0;
2388 constexpr unsigned int y_coordinate = (dim >= 2 ? 1 : 0);
2389 constexpr unsigned int z_coordinate = (dim >= 3 ? 2 : 0);
2390
2391 switch (this->kind)
2392 {
2394 {
2395 // Vertices are special cases in that they do not actually
2396 // have coordinates. Error out if this function is called
2397 // with a vertex:
2398 Assert(false,
2399 ExcMessage("Vertices are zero-dimensional objects and "
2400 "as a consequence have no coordinates. You "
2401 "cannot meaningfully ask whether a point is "
2402 "inside a vertex (within a certain tolerance) "
2403 "without coordinate values."));
2404 return false;
2405 }
2409 {
2410 for (unsigned int d = 0; d < dim; ++d)
2411 if ((p[d] < -tolerance) || (p[d] > 1 + tolerance))
2412 return false;
2413 return true;
2414 }
2417 {
2418 // First make sure that we are in the first quadrant or octant
2419 for (unsigned int d = 0; d < dim; ++d)
2420 if (p[d] < -tolerance)
2421 return false;
2422
2423 // Now we also need to make sure that we are below the diagonal line
2424 // or plane that delineates the simplex. This diagonal is given by
2425 // sum(p[d])<=1, and a diagonal a distance eps away is given by
2426 // sum(p[d])<=1+eps*sqrt(d). (For example, the point at (1,1) is a
2427 // distance of 1/sqrt(2) away from the diagonal. That is, its
2428 // sum satisfies
2429 // sum(p[d]) = 2 <= 1 + (1/sqrt(2)) * sqrt(2)
2430 // in other words, it satisfies the predicate with eps=1/sqrt(2).)
2431 double sum = 0;
2432 for (unsigned int d = 0; d < dim; ++d)
2433 sum += p[d];
2434 return (sum <= 1 + tolerance * std::sqrt(1. * dim));
2435 }
2437 {
2438 // A pyramid only lives in the upper half-space:
2439 if (p[z_coordinate] < -tolerance)
2440 return false;
2441
2442 // It also only lives in the space below z=1:
2443 if (p[z_coordinate] > 1 + tolerance)
2444 return false;
2445
2446 // Within what's left of the space, a pyramid is a cone that tapers
2447 // towards the top. First compute the distance of the point to the
2448 // axis in the max norm (this is the right norm because the vertices
2449 // of the pyramid are at points +/-1, +/-1):
2450 const double distance_from_axis =
2451 std::max(std::fabs(p[x_coordinate]), std::fabs(p[y_coordinate]));
2452
2453 // We are inside the pyramid if the distance from the axis is less
2454 // than (1-z)
2455 return (distance_from_axis <= 1 + tolerance - p[z_coordinate]);
2456 }
2458 {
2459 // The wedge we use is a triangle extruded into the third
2460 // dimension by one unit. So we can use the same logic as for
2461 // triangles above (i.e., for the simplex above, using dim==2)
2462 // and then check the third dimension separately.
2463
2464 if ((p[x_coordinate] < -tolerance) || (p[y_coordinate] < -tolerance))
2465 return false;
2466
2467 const double sum = p[x_coordinate] + p[y_coordinate];
2468 if (sum > 1 + tolerance * std::sqrt(2.0))
2469 return false;
2470
2471 if (p[z_coordinate] < -tolerance)
2472 return false;
2473 if (p[z_coordinate] > 1 + tolerance)
2474 return false;
2475
2476 return true;
2477 }
2478 default:
2480 }
2481
2482 return false;
2483}
2484
2485
2486
2487template <int dim>
2488inline Tensor<1, dim>
2490 const unsigned int i) const
2491{
2493 AssertIndexRange(i, dim - 1);
2494
2495 switch (this->kind)
2496 {
2504 {
2505 AssertIndexRange(face_no, 3);
2506 static const std::array<Tensor<1, dim>, 3> table = {
2507 {Point<dim>(1, 0),
2508 Point<dim>(-std::sqrt(0.5), +std::sqrt(0.5)),
2509 Point<dim>(0, -1)}};
2510
2511 return table[face_no];
2512 }
2514 {
2515 AssertIndexRange(face_no, 4);
2516 static const ndarray<Tensor<1, dim>, 4, 2> table = {
2517 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
2518 {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
2519 {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
2520 {{Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
2521 +std::pow(1.0 / 3.0, 1.0 / 4.0),
2522 0),
2523 Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
2524 0,
2525 +std::pow(1.0 / 3.0, 1.0 / 4.0))}}}};
2526
2527 return table[face_no][i];
2528 }
2530 {
2531 AssertIndexRange(face_no, 5);
2532 static const ndarray<Tensor<1, dim>, 5, 2> table = {
2533 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
2534 {{Point<dim>(+1.0 / sqrt(2.0), 0, +1.0 / sqrt(2.0)),
2535 Point<dim>(0, 1, 0)}},
2536 {{Point<dim>(+1.0 / sqrt(2.0), 0, -1.0 / sqrt(2.0)),
2537 Point<dim>(0, 1, 0)}},
2538 {{Point<dim>(1, 0, 0),
2539 Point<dim>(0, +1.0 / sqrt(2.0), +1.0 / sqrt(2.0))}},
2540 {{Point<dim>(1, 0, 0),
2541 Point<dim>(0, +1.0 / sqrt(2.0), -1.0 / sqrt(2.0))}}}};
2542
2543 return table[face_no][i];
2544 }
2546 {
2547 AssertIndexRange(face_no, 5);
2548 static const ndarray<Tensor<1, dim>, 5, 2> table = {
2549 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
2550 {{Point<dim>(1, 0, 0), Point<dim>(0, 1, 0)}},
2551 {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
2552 {{Point<dim>(-1 / std::sqrt(2.0), +1 / std::sqrt(2.0), 0),
2553 Point<dim>(0, 0, 1)}},
2554 {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}}}};
2555
2556 return table[face_no][i];
2557 }
2558 default:
2560 }
2561
2562
2563 return {};
2564}
2565
2566
2567
2568template <int dim>
2569inline Tensor<1, dim>
2570ReferenceCell::unit_normal_vectors(const unsigned int face_no) const
2571{
2572 AssertDimension(dim, this->get_dimension());
2573
2574 if (is_hyper_cube())
2575 {
2578 }
2579 else if (dim == 2)
2580 {
2582
2583 // Return the rotated vector
2584 return cross_product_2d(unit_tangential_vectors<dim>(face_no, 0));
2585 }
2586 else if (dim == 3)
2587 {
2588 return cross_product_3d(unit_tangential_vectors<dim>(face_no, 0),
2589 unit_tangential_vectors<dim>(face_no, 1));
2590 }
2591
2593
2594 return {};
2595}
2596
2597
2598
2599inline unsigned int
2600ReferenceCell::n_face_orientations(const unsigned int face_no) const
2601{
2602 AssertIndexRange(face_no, n_faces());
2603 if (get_dimension() == 1)
2604 return 1;
2605 if (get_dimension() == 2)
2606 return 2;
2608 return 8;
2609 else if (face_reference_cell(face_no) == ReferenceCells::Triangle)
2610 return 6;
2611
2614}
2615
2616
2617
2618inline bool
2620 const unsigned int line,
2621 const unsigned int face,
2622 const unsigned char combined_face_orientation,
2623 const bool line_orientation) const
2624{
2625 if (*this == ReferenceCells::Hexahedron)
2626 {
2627 static constexpr ::ndarray<bool, 2, 8> bool_table{
2628 {{{true, true, false, true, false, false, true, false}},
2629 {{true, true, true, false, false, false, false, true}}}};
2630
2631 return (line_orientation ==
2632 bool_table[line / 2][combined_face_orientation]);
2633 }
2634 else if (*this == ReferenceCells::Tetrahedron)
2635 {
2636 static constexpr unsigned int X = numbers::invalid_unsigned_int;
2637 static constexpr ::ndarray<unsigned int, 4, 3> combined_lines{
2638 {{{0, 0, 0}}, {{X, 0, 1}}, {{X, 0, X}}, {{X, X, X}}}};
2639
2640 const auto combined_line = combined_lines[face][line];
2641
2642 Assert(combined_line != X,
2643 ExcMessage(
2644 "This function can only be called for following face-line "
2645 "combinations: (0,0), (0,1), (0,2), (1,1), (1,2), (2,1),"));
2646
2647 static constexpr ::ndarray<bool, 2, 6> bool_table{
2648 {{{false, true, false, true, false, true}},
2649 {{true, false, true, false, true, false}}}};
2650
2651 return (line_orientation ==
2652 bool_table[combined_line][combined_face_orientation]);
2653 }
2654 else
2655 // TODO: This might actually be wrong for some of the other
2656 // kinds of objects. We should check this
2657 return true;
2658}
2659
2660
2661
2662namespace internal
2663{
2664 template <typename T>
2666 {
2667 public:
2681
2685 virtual ~NoPermutation() noexcept override = default;
2686
2690 virtual void
2691 print_info(std::ostream &out) const override
2692 {
2693 out << '[';
2694
2695 const unsigned int n_vertices = entity_type.n_vertices();
2696
2697 for (unsigned int i = 0; i < n_vertices; ++i)
2698 {
2699 out << vertices_0[i];
2700 if (i + 1 != n_vertices)
2701 out << ',';
2702 }
2703
2704 out << "] is not a valid permutation of [";
2705
2706 for (unsigned int i = 0; i < n_vertices; ++i)
2707 {
2708 out << vertices_1[i];
2709 if (i + 1 != n_vertices)
2710 out << ',';
2711 }
2712
2713 out << "]." << std::endl;
2714 }
2715
2720
2725
2730 };
2731
2743 << "The reference-cell type used on this cell (" << arg1.to_string()
2744 << ") does not match the reference-cell type of the finite element "
2745 << "associated with this cell (" << arg2.to_string() << "). "
2746 << "Did you accidentally use simplex elements on hypercube meshes "
2747 << "(or the other way around), or are you using a mixed mesh and "
2748 << "assigned a simplex element to a hypercube cell (or the other "
2749 << "way around) via the active_fe_index?");
2750} // namespace internal
2751
2752
2753
2754template <typename T, std::size_t N>
2755inline unsigned char
2756ReferenceCell::compute_orientation(const std::array<T, N> &vertices_0,
2757 const std::array<T, N> &vertices_1) const
2758{
2759 Assert(N >= n_vertices(),
2760 ExcMessage("The number of array elements must be equal to or "
2761 "greater than the number of vertices of the cell "
2762 "referenced by this object."));
2763
2764 // Call the non-deprecated function, taking care of calling it only with
2765 // those array elements that we actually care about (see the note
2766 // in the documentation about the arguments potentially being
2767 // larger arrays than necessary).
2769 make_array_view(vertices_0.begin(), vertices_0.begin() + n_vertices()),
2770 make_array_view(vertices_1.begin(), vertices_1.begin() + n_vertices()));
2771}
2772
2773
2774
2775template <typename T>
2776unsigned char
2778 const ArrayView<const T> &vertices_0,
2779 const ArrayView<const T> &vertices_1) const
2780{
2781 Assert(vertices_0.size() == n_vertices(),
2782 ExcMessage("The number of array elements must be equal to "
2783 "the number of vertices of the cell "
2784 "referenced by this object."));
2785 Assert(vertices_1.size() == n_vertices(),
2786 ExcMessage("The number of array elements must be equal to "
2787 "the number of vertices of the cell "
2788 "referenced by this object."));
2789
2790 const auto v0_equals = [&](const std::initializer_list<const T> &list) {
2791 Assert(list.size() == n_vertices(), ExcInternalError());
2792 return std::equal(vertices_0.begin(), vertices_0.end(), std::begin(list));
2793 };
2794
2795 switch (this->kind)
2796 {
2798 // Things are always default-oriented in 1D
2799 if (v0_equals({vertices_1[0]}))
2801 break;
2802
2804 // line_orientation=true
2805 if (v0_equals({vertices_1[0], vertices_1[1]}))
2807
2808 // line_orientation=false
2809 if (v0_equals({vertices_1[1], vertices_1[0]}))
2811 break;
2813 // face_orientation=true, face_rotation=false, face_flip=false
2814 if (v0_equals({vertices_1[0], vertices_1[1], vertices_1[2]}))
2815 return 1;
2816
2817 // face_orientation=true, face_rotation=true, face_flip=false
2818 if (v0_equals({vertices_1[1], vertices_1[2], vertices_1[0]}))
2819 return 5;
2820
2821 // face_orientation=true, face_rotation=false, face_flip=true
2822 if (v0_equals({vertices_1[2], vertices_1[0], vertices_1[1]}))
2823 return 3;
2824
2825 // face_orientation=false, face_rotation=false, face_flip=false
2826 if (v0_equals({vertices_1[0], vertices_1[2], vertices_1[1]}))
2827 return 0;
2828
2829 // face_orientation=false, face_rotation=true, face_flip=false
2830 if (v0_equals({vertices_1[2], vertices_1[1], vertices_1[0]}))
2831 return 2;
2832
2833 // face_orientation=false, face_rotation=false, face_flip=true
2834 if (v0_equals({vertices_1[1], vertices_1[0], vertices_1[2]}))
2835 return 4;
2836 break;
2838 // face_orientation=true, face_rotation=false, face_flip=false
2839 if (v0_equals(
2840 {vertices_1[0], vertices_1[1], vertices_1[2], vertices_1[3]}))
2841 return 1;
2842
2843 // face_orientation=true, face_rotation=true, face_flip=false
2844 if (v0_equals(
2845 {vertices_1[2], vertices_1[0], vertices_1[3], vertices_1[1]}))
2846 return 3;
2847
2848 // face_orientation=true, face_rotation=false, face_flip=true
2849 if (v0_equals(
2850 {vertices_1[3], vertices_1[2], vertices_1[1], vertices_1[0]}))
2851 return 5;
2852
2853 // face_orientation=true, face_rotation=true, face_flip=true
2854 if (v0_equals(
2855 {vertices_1[1], vertices_1[3], vertices_1[0], vertices_1[2]}))
2856 return 7;
2857
2858 // face_orientation=false, face_rotation=false, face_flip=false
2859 if (v0_equals(
2860 {vertices_1[0], vertices_1[2], vertices_1[1], vertices_1[3]}))
2861 return 0;
2862
2863 // face_orientation=false, face_rotation=true, face_flip=false
2864 if (v0_equals(
2865 {vertices_1[2], vertices_1[3], vertices_1[0], vertices_1[1]}))
2866 return 2;
2867
2868 // face_orientation=false, face_rotation=false, face_flip=true
2869 if (v0_equals(
2870 {vertices_1[3], vertices_1[1], vertices_1[2], vertices_1[0]}))
2871 return 4;
2872
2873 // face_orientation=false, face_rotation=true, face_flip=true
2874 if (v0_equals(
2875 {vertices_1[1], vertices_1[0], vertices_1[3], vertices_1[2]}))
2876 return 6;
2877 break;
2878 default:
2880 }
2881
2882 Assert(false, (internal::NoPermutation<T>(*this, vertices_0, vertices_1)));
2883 return std::numeric_limits<unsigned char>::max();
2884}
2885
2886
2887
2888template <typename T, std::size_t N>
2889inline std::array<T, N>
2891 const std::array<T, N> &vertices,
2892 const unsigned int orientation) const
2893{
2894 Assert(N >= n_vertices(),
2895 ExcMessage("The number of array elements must be equal to or "
2896 "greater than the number of vertices of the cell "
2897 "referenced by this object."));
2898
2899 // Call the non-deprecated function, taking care of calling it only with
2900 // those array elements that we actually care about (see the note
2901 // in the documentation about the arguments potentially being
2902 // larger arrays than necessary).
2903 const auto permutation = permute_by_combined_orientation(
2904 make_array_view(vertices.begin(), vertices.begin() + n_vertices()),
2905 orientation);
2906
2907 std::array<T, N> temp;
2908 std::copy(permutation.begin(), permutation.end(), temp.begin());
2909
2910 return temp;
2911}
2912
2913
2914
2915template <typename T>
2916boost::container::small_vector<T, 8>
2919 const unsigned char orientation) const
2920{
2921 Assert(vertices.size() == n_vertices(),
2922 ExcMessage("The number of array elements must be equal to "
2923 "the number of vertices of the cell "
2924 "referenced by this object."));
2925
2926 switch (this->kind)
2927 {
2929 switch (orientation)
2930 {
2931 case 1:
2932 return {vertices[0], vertices[1]};
2933 case 0:
2934 return {vertices[1], vertices[0]};
2935 default:
2937 }
2938 break;
2940 switch (orientation)
2941 {
2942 case 1:
2943 return {vertices[0], vertices[1], vertices[2]};
2944 case 3:
2945 return {vertices[1], vertices[2], vertices[0]};
2946 case 5:
2947 return {vertices[2], vertices[0], vertices[1]};
2948 case 0:
2949 return {vertices[0], vertices[2], vertices[1]};
2950 case 2:
2951 return {vertices[2], vertices[1], vertices[0]};
2952 case 4:
2953 return {vertices[1], vertices[0], vertices[2]};
2954 default:
2956 }
2957 break;
2959 switch (orientation)
2960 {
2961 case 1:
2962 return {vertices[0], vertices[1], vertices[2], vertices[3]};
2963 case 3:
2964 return {vertices[2], vertices[0], vertices[3], vertices[1]};
2965 case 5:
2966 return {vertices[3], vertices[2], vertices[1], vertices[0]};
2967 case 7:
2968 return {vertices[1], vertices[3], vertices[0], vertices[2]};
2969 case 0:
2970 return {vertices[0], vertices[2], vertices[1], vertices[3]};
2971 case 2:
2972 return {vertices[2], vertices[3], vertices[0], vertices[1]};
2973 case 4:
2974 return {vertices[3], vertices[1], vertices[2], vertices[0]};
2975 case 6:
2976 return {vertices[1], vertices[0], vertices[3], vertices[2]};
2977 default:
2979 }
2980 break;
2981 default:
2983 }
2984
2985 return {};
2986}
2987
2988
2989
2990inline unsigned char
2992 const unsigned char orientation) const
2993{
2994 switch (this->kind)
2995 {
2997 // Things are always default-oriented in 1D
2998 return orientation;
2999
3001 // the 1d orientations are the identity and a flip: i.e., the identity
3002 // and an involutory mapping
3003 return orientation;
3004
3006 {
3007 AssertIndexRange(orientation, 6);
3008 constexpr std::array<unsigned char, 6> inverses{{0, 1, 2, 5, 4, 3}};
3009 return inverses[orientation];
3010 }
3012 {
3013 AssertIndexRange(orientation, 8);
3014 constexpr std::array<unsigned char, 8> inverses{
3015 {0, 1, 2, 7, 4, 5, 6, 3}};
3016 return inverses[orientation];
3017 }
3018 default:
3020 }
3021
3022 return std::numeric_limits<unsigned char>::max();
3023}
3024
3025
3026
3027template <>
3028unsigned int
3029ReferenceCell::vtk_lexicographic_to_node_index<0>(
3030 const std::array<unsigned, 0> &node_indices,
3031 const std::array<unsigned, 0> &nodes_per_direction,
3032 const bool legacy_format) const;
3033
3034template <>
3035unsigned int
3036ReferenceCell::vtk_lexicographic_to_node_index<1>(
3037 const std::array<unsigned, 1> &node_indices,
3038 const std::array<unsigned, 1> &nodes_per_direction,
3039 const bool legacy_format) const;
3040
3041template <>
3042unsigned int
3043ReferenceCell::vtk_lexicographic_to_node_index<2>(
3044 const std::array<unsigned, 2> &node_indices,
3045 const std::array<unsigned, 2> &nodes_per_direction,
3046 const bool legacy_format) const;
3047
3048template <>
3049unsigned int
3050ReferenceCell::vtk_lexicographic_to_node_index<3>(
3051 const std::array<unsigned, 3> &node_indices,
3052 const std::array<unsigned, 3> &nodes_per_direction,
3053 const bool legacy_format) const;
3054
3056
3057#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
iterator begin() const
Definition array_view.h:702
iterator end() const
Definition array_view.h:711
std::size_t size() const
Definition array_view.h:684
Abstract base class for mapping classes.
Definition mapping.h:316
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
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)
static ::ExceptionBase & ExcNotImplemented()
#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 AssertThrow(cond, exc)
#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)