Reference documentation for deal.II version GIT 9042b9283b 2023-12-02 14: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\}}\)
reference_cell.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2020 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_tria_reference_cell_h
17 #define dealii_tria_reference_cell_h
18 
19 #include <deal.II/base/config.h>
20 
23 #include <deal.II/base/ndarray.h>
24 #include <deal.II/base/point.h>
25 #include <deal.II/base/tensor.h>
26 #include <deal.II/base/utilities.h>
27 
29 
30 #include <boost/container/small_vector.hpp>
31 
32 #include <iosfwd>
33 #include <string>
34 
36 
37 // Forward declarations
38 #ifndef DOXYGEN
39 template <int dim, int spacedim>
40 class Mapping;
41 
42 template <int dim>
43 class Quadrature;
44 
45 class ReferenceCell;
46 #endif
47 
48 
49 namespace internal
50 {
63  constexpr ReferenceCell
64  make_reference_cell_from_int(const std::uint8_t kind);
65 } // namespace internal
66 
67 
68 
100 {
101 public:
109  static ReferenceCell
110  n_vertices_to_type(const int dim, const unsigned int n_vertices);
111 
120  constexpr ReferenceCell();
121 
132  bool
133  is_hyper_cube() const;
134 
138  bool
139  is_simplex() const;
140 
145  unsigned int
146  get_dimension() const;
147 
161  template <int dim>
162  double
163  d_linear_shape_function(const Point<dim> &xi, const unsigned int i) const;
164 
169  template <int dim>
172  const unsigned int i) const;
173 
183  template <int dim, int spacedim = dim>
184  std::unique_ptr<Mapping<dim, spacedim>>
185  get_default_mapping(const unsigned int degree) const;
186 
198  template <int dim, int spacedim = dim>
199  const Mapping<dim, spacedim> &
201 
210  template <int dim>
212  get_gauss_type_quadrature(const unsigned n_points_1d) const;
213 
223  template <int dim>
225  get_midpoint_quadrature() const;
226 
240  template <int dim>
241  const Quadrature<dim> &
243 
258  unsigned int
259  n_vertices() const;
260 
266  vertex_indices() const;
267 
275  template <int dim>
276  Point<dim>
277  vertex(const unsigned int v) const;
278 
284  unsigned int
285  n_lines() const;
286 
292  line_indices() const;
293 
299  unsigned int
300  n_faces() const;
301 
307  face_indices() const;
308 
320  unsigned int
321  n_isotropic_children() const;
322 
328  isotropic_child_indices() const;
329 
343  face_reference_cell(const unsigned int face_no) const;
344 
358  static constexpr unsigned char
360 
371  static constexpr unsigned char
373 
415  unsigned int
416  child_cell_on_face(const unsigned int face,
417  const unsigned int subface,
418  const unsigned char face_orientation =
420 
430  std::array<unsigned int, 2>
431  standard_vertex_to_face_and_vertex_index(const unsigned int vertex) const;
432 
442  std::array<unsigned int, 2>
443  standard_line_to_face_and_line_index(const unsigned int line) const;
444 
453  unsigned int
454  line_to_cell_vertices(const unsigned int line,
455  const unsigned int vertex) const;
456 
460  unsigned int
461  face_to_cell_lines(const unsigned int face,
462  const unsigned int line,
463  const unsigned char face_orientation) const;
464 
468  unsigned int
469  face_to_cell_vertices(const unsigned int face,
470  const unsigned int vertex,
471  const unsigned char face_orientation) const;
472 
492  template <int dim>
493  Point<dim>
494  face_vertex_location(const unsigned int face,
495  const unsigned int vertex) const;
496 
500  unsigned int
501  standard_to_real_face_vertex(const unsigned int vertex,
502  const unsigned int face,
503  const unsigned char face_orientation) const;
504 
508  unsigned int
509  standard_to_real_face_line(const unsigned int line,
510  const unsigned int face,
511  const unsigned char face_orientation) const;
512 
523  bool
524  standard_vs_true_line_orientation(const unsigned int line,
525  const unsigned int face,
526  const unsigned char face_orientation,
527  const bool line_orientation) const;
528 
552  double
553  volume() const;
554 
567  template <int dim>
568  Point<dim>
569  barycenter() const;
570 
590  template <int dim>
591  bool
592  contains_point(const Point<dim> &p, const double tolerance = 0) const;
593 
598  template <int dim>
599  Point<dim>
600  closest_point(const Point<dim> &p) const;
601 
609  template <int dim>
611  unit_tangential_vectors(const unsigned int face_no,
612  const unsigned int i) const;
613 
617  template <int dim>
619  unit_normal_vectors(const unsigned int face_no) const;
620 
626  unsigned int
627  n_face_orientations(const unsigned int face_no) const;
628 
645  template <typename T, std::size_t N>
646  DEAL_II_DEPRECATED unsigned char
647  compute_orientation(const std::array<T, N> &vertices_0,
648  const std::array<T, N> &vertices_1) const;
649 
675  template <typename T>
676  unsigned char
678  const ArrayView<const T> &vertices_1) const;
679 
680 
694  template <typename T, std::size_t N>
695  DEAL_II_DEPRECATED std::array<T, N>
696  permute_according_orientation(const std::array<T, N> &vertices,
697  const unsigned int orientation) const;
698 
710  template <typename T>
711  boost::container::small_vector<T, 8>
713  const unsigned char orientation) const;
714 
719  faces_for_given_vertex(const unsigned int vertex_index) const;
720 
733  unsigned int
734  exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const;
735 
739  unsigned int
740  exodusii_face_to_deal_face(const unsigned int face_n) const;
741 
745  unsigned int
746  unv_vertex_to_deal_vertex(const unsigned int vertex_n) const;
747 
751  unsigned int
752  vtk_linear_type() const;
753 
758  unsigned int
759  vtk_quadratic_type() const;
760 
765  unsigned int
766  vtk_lagrange_type() const;
767 
780  template <int dim>
781  unsigned int
783  const std::array<unsigned, dim> &node_indices,
784  const std::array<unsigned, dim> &nodes_per_direction,
785  const bool legacy_format) const;
786 
790  unsigned int
791  vtk_vertex_to_deal_vertex(const unsigned int vertex_index) const;
792 
796  unsigned int
797  gmsh_element_type() const;
798 
812  std::string
813  to_string() const;
814 
818  constexpr operator std::uint8_t() const;
819 
823  constexpr bool
824  operator==(const ReferenceCell &type) const;
825 
829  constexpr bool
830  operator!=(const ReferenceCell &type) const;
831 
837  template <class Archive>
838  void
839  serialize(Archive &archive, const unsigned int /*version*/);
840 
844  static constexpr std::size_t
846 
851 private:
855  std::uint8_t kind;
856 
863  constexpr ReferenceCell(const std::uint8_t kind);
864 
869  {{{1, 0}}, {{0, 1}}}};
870 
871 
876  {{{0, 2, 1}},
877  {{0, 1, 2}},
878  {{2, 1, 0}},
879  {{2, 0, 1}},
880  {{1, 0, 2}},
881  {{1, 2, 0}}}};
882 
886  static constexpr ndarray<unsigned int, 8, 4>
888  {{0, 1, 2, 3}},
889  {{2, 3, 0, 1}},
890  {{2, 0, 3, 1}},
891  {{3, 1, 2, 0}},
892  {{3, 2, 1, 0}},
893  {{1, 0, 3, 2}},
894  {{1, 3, 0, 2}}}};
895 
900  friend constexpr ReferenceCell
902 
903  friend std::ostream &
904  operator<<(std::ostream &out, const ReferenceCell &reference_cell);
905 
906  friend std::istream &
907  operator>>(std::istream &in, ReferenceCell &reference_cell);
908 };
909 
910 
918 std::ostream &
919 operator<<(std::ostream &out, const ReferenceCell &reference_cell);
920 
928 std::istream &
929 operator>>(std::istream &in, ReferenceCell &reference_cell);
930 
931 
932 inline constexpr ReferenceCell::ReferenceCell(const std::uint8_t kind)
933  : kind(kind)
934 {}
935 
936 
937 
938 inline constexpr ReferenceCell::operator std::uint8_t() const
939 {
940  return kind;
941 }
942 
943 
944 
945 inline constexpr bool
947 {
948  return kind == type.kind;
949 }
950 
951 
952 
953 inline constexpr bool
955 {
956  return kind != type.kind;
957 }
958 
959 
960 
961 namespace internal
962 {
963  inline constexpr ReferenceCell
964  make_reference_cell_from_int(const std::uint8_t kind)
965  {
966 #ifndef DEAL_II_CXX14_CONSTEXPR_BUG
967  // Make sure these are the only indices from which objects can be
968  // created.
969  Assert((kind == std::numeric_limits<std::uint8_t>::max()) || (kind < 8),
970  ExcInternalError());
971 #endif
972 
973  // Call the private constructor, which we can from here because this
974  // function is a 'friend'.
975  return {kind};
976  }
977 } // namespace internal
978 
979 
980 
988 namespace ReferenceCells
989 {
990  constexpr const ReferenceCell Vertex =
992  constexpr const ReferenceCell Line =
994  constexpr const ReferenceCell Triangle =
996  constexpr const ReferenceCell Quadrilateral =
998  constexpr const ReferenceCell Tetrahedron =
1000  constexpr const ReferenceCell Pyramid =
1002  constexpr const ReferenceCell Wedge =
1004  constexpr const ReferenceCell Hexahedron =
1006  constexpr const ReferenceCell Invalid =
1009 
1015  template <int dim>
1016  constexpr const ReferenceCell &
1017  get_simplex();
1018 
1024  template <int dim>
1025  constexpr const ReferenceCell &
1026  get_hypercube();
1027 } // namespace ReferenceCells
1028 
1029 
1030 
1033 {}
1034 
1035 
1036 
1037 template <class Archive>
1038 inline void
1039 ReferenceCell::serialize(Archive &archive, const unsigned int /*version*/)
1040 {
1041  archive &kind;
1042 }
1043 
1044 
1045 
1046 inline constexpr std::size_t
1048 {
1049  return sizeof(ReferenceCells::Invalid);
1050 }
1051 
1052 
1053 
1055 ReferenceCell::faces_for_given_vertex(const unsigned int vertex) const
1056 {
1058  switch (this->kind)
1059  {
1060  case ReferenceCells::Line:
1061  return {&GeometryInfo<2>::vertex_to_face[vertex][0], 1};
1063  return {&GeometryInfo<2>::vertex_to_face[vertex][0], 2};
1065  {
1066  static constexpr ndarray<unsigned int, 3, 2> table = {
1067  {{{0, 2}}, {{0, 1}}, {{1, 2}}}};
1068  return table[vertex];
1069  }
1071  {
1072  static constexpr ndarray<unsigned int, 4, 3> table = {
1073  {{{0, 1, 2}}, {{0, 1, 3}}, {{0, 2, 3}}, {{1, 2, 3}}}};
1074 
1075  return table[vertex];
1076  }
1078  {
1079  static constexpr unsigned int X = numbers::invalid_unsigned_int;
1080  static constexpr ndarray<unsigned int, 5, 4> table = {
1081  {{{0, 1, 3, X}},
1082  {{0, 2, 3, X}},
1083  {{0, 1, 4, X}},
1084  {{0, 2, 4, X}},
1085  {{1, 2, 3, 4}}}};
1086 
1087  return {&table[vertex][0], vertex == 4 ? 4u : 3u};
1088  }
1089  case ReferenceCells::Wedge:
1090  {
1092  static constexpr ndarray<unsigned int, 6, 3> table = {{{{0, 2, 4}},
1093  {{0, 2, 3}},
1094  {{0, 3, 4}},
1095  {{1, 2, 4}},
1096  {{1, 2, 3}},
1097  {{1, 3, 4}}}};
1098 
1099  return table[vertex];
1100  }
1102  return {&GeometryInfo<3>::vertex_to_face[vertex][0], 3};
1103  default:
1104  Assert(false, ExcNotImplemented());
1105  }
1106 
1107  return {};
1108 }
1109 
1110 
1111 
1112 inline bool
1114 {
1115  return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
1116  *this == ReferenceCells::Quadrilateral ||
1117  *this == ReferenceCells::Hexahedron);
1118 }
1119 
1120 
1121 
1122 inline bool
1124 {
1125  return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
1126  *this == ReferenceCells::Triangle ||
1127  *this == ReferenceCells::Tetrahedron);
1128 }
1129 
1130 
1131 
1132 inline unsigned int
1134 {
1135  switch (this->kind)
1136  {
1138  return 0;
1139  case ReferenceCells::Line:
1140  return 1;
1143  return 2;
1146  case ReferenceCells::Wedge:
1148  return 3;
1149  default:
1150  Assert(false, ExcNotImplemented());
1151  }
1152 
1154 }
1155 
1156 
1157 
1158 template <int dim>
1161 {
1162  return Quadrature<dim>(std::vector<Point<dim>>({barycenter<dim>()}),
1163  std::vector<double>({volume()}));
1164 }
1165 
1166 
1167 
1168 inline unsigned int
1170 {
1171  switch (this->kind)
1172  {
1174  return 1;
1175  case ReferenceCells::Line:
1176  return 2;
1178  return 3;
1180  return 4;
1182  return 4;
1184  return 5;
1185  case ReferenceCells::Wedge:
1186  return 6;
1188  return 8;
1189  default:
1190  Assert(false, ExcNotImplemented());
1191  }
1192 
1194 }
1195 
1196 
1197 
1198 inline unsigned int
1200 {
1201  switch (this->kind)
1202  {
1204  return 0;
1205  case ReferenceCells::Line:
1206  return 1;
1208  return 3;
1210  return 4;
1212  return 6;
1214  return 8;
1215  case ReferenceCells::Wedge:
1216  return 9;
1218  return 12;
1219  default:
1220  Assert(false, ExcNotImplemented());
1221  }
1222 
1224 }
1225 
1226 
1227 
1228 template <int dim>
1229 Point<dim>
1230 ReferenceCell::vertex(const unsigned int v) const
1231 {
1234 
1235  switch (dim)
1236  {
1237  case 0:
1238  {
1239  if (*this == ReferenceCells::Vertex)
1240  return Point<dim>();
1241  break;
1242  }
1243  case 1:
1244  {
1245  static const Point<dim> vertices[2] = {
1246  Point<dim>(), // the origin
1247  Point<dim>::unit_vector(0) // unit point along x-axis
1248  };
1249  if (*this == ReferenceCells::Line)
1250  return vertices[v];
1251  break;
1252  }
1253  case 2:
1254  {
1255  switch (this->kind)
1256  {
1258  {
1259  static const Point<dim> vertices[3] = {
1260  Point<dim>(), // the origin
1261  Point<dim>::unit_vector(0), // unit point along x-axis
1262  Point<dim>::unit_vector(1) // unit point along y-axis
1263  };
1264  return vertices[v];
1265  }
1267  {
1268  static const Point<dim> vertices[4] = {
1269  // First the two points on the x-axis
1270  Point<dim>(),
1272  // Then these two points shifted in the y-direction
1275  return vertices[v];
1276  }
1277  }
1278  break;
1279  }
1280  case 3:
1281  {
1282  switch (this->kind)
1283  {
1285  {
1286  static const Point<dim> vertices[4] = {
1287  Point<dim>(), // the origin
1288  Point<dim>::unit_vector(0), // unit point along x-axis
1289  Point<dim>::unit_vector(1), // unit point along y-axis
1290  Point<dim>::unit_vector(2) // unit point along z-axis
1291  };
1292  return vertices[v];
1293  }
1295  {
1296  static const Point<dim> vertices[5] = {
1297  Point<dim>{-1.0, -1.0, 0.0},
1298  Point<dim>{+1.0, -1.0, 0.0},
1299  Point<dim>{-1.0, +1.0, 0.0},
1300  Point<dim>{+1.0, +1.0, 0.0},
1301  Point<dim>{+0.0, +0.0, 1.0}};
1302  return vertices[v];
1303  }
1304  case ReferenceCells::Wedge:
1305  {
1306  static const Point<dim> vertices[6] = {
1307  // First the three points on the triangular base of the
1308  // wedge:
1309  Point<dim>(),
1312  // And now everything shifted in the z-direction again
1316  return vertices[v];
1317  }
1319  {
1320  static const Point<dim> vertices[8] = {
1321  // First the two points on the x-axis
1322  Point<dim>(),
1324  // Then these two points shifted in the y-direction
1327  // And now all four points shifted in the z-direction
1334  return vertices[v];
1335  }
1336  }
1337  break;
1338  }
1339  default:
1340  Assert(false, ExcNotImplemented());
1341  }
1342 
1343  Assert(false, ExcNotImplemented());
1344  return Point<dim>();
1345 }
1346 
1347 
1348 inline unsigned int
1350 {
1351  switch (this->kind)
1352  {
1354  return 0;
1355  case ReferenceCells::Line:
1356  return 2;
1358  return 3;
1360  return 4;
1362  return 4;
1364  return 5;
1365  case ReferenceCells::Wedge:
1366  return 5;
1368  return 6;
1369  default:
1370  Assert(false, ExcNotImplemented());
1371  }
1372 
1374 }
1375 
1376 
1377 
1380 {
1381  return {0U, n_faces()};
1382 }
1383 
1384 
1385 
1386 inline unsigned int
1388 {
1389  switch (this->kind)
1390  {
1392  return 0;
1393  case ReferenceCells::Line:
1394  return 2;
1396  return 4;
1398  return 4;
1400  return 8;
1402  // We haven't yet decided how to refine pyramids. Update this when we
1403  // have
1404  Assert(false, ExcNotImplemented());
1406  case ReferenceCells::Wedge:
1407  return 8;
1409  return 8;
1410  default:
1411  Assert(false, ExcNotImplemented());
1412  }
1413 
1415 }
1416 
1417 
1418 
1421 {
1422  return {0U, n_isotropic_children()};
1423 }
1424 
1425 
1426 
1429 {
1430  return {0U, n_vertices()};
1431 }
1432 
1433 
1434 
1437 {
1438  return {0U, n_lines()};
1439 }
1440 
1441 
1442 
1443 inline ReferenceCell
1444 ReferenceCell::face_reference_cell(const unsigned int face_no) const
1445 {
1446  AssertIndexRange(face_no, n_faces());
1447 
1448  switch (this->kind)
1449  {
1451  return ReferenceCells::Invalid;
1452  case ReferenceCells::Line:
1453  return ReferenceCells::Vertex;
1456  return ReferenceCells::Line;
1458  return ReferenceCells::Triangle;
1460  if (face_no == 0)
1462  else
1463  return ReferenceCells::Triangle;
1464  case ReferenceCells::Wedge:
1465  if (face_no > 1)
1467  else
1468  return ReferenceCells::Triangle;
1471  default:
1472  Assert(false, ExcNotImplemented());
1473  }
1474 
1475  return ReferenceCells::Invalid;
1476 }
1477 
1478 
1479 
1480 inline constexpr unsigned char
1482 {
1483  // Our convention is that 'orientation' has a default value of true and
1484  // occupies the least-significant bit while rotate and flip have default
1485  // values of 'false' and occupy the second and third bits.
1486  return 0b001;
1487 }
1488 
1489 
1490 
1491 inline constexpr unsigned char
1493 {
1494  // For a reversed line 'orientation' is false and neither flip nor rotate are
1495  // defined.
1496  return 0b000;
1497 }
1498 
1499 
1500 
1501 inline unsigned int
1503  const unsigned int face,
1504  const unsigned int subface,
1505  const unsigned char combined_face_orientation) const
1506 {
1507  AssertIndexRange(face, n_faces());
1509 
1510  switch (this->kind)
1511  {
1513  case ReferenceCells::Line:
1514  {
1515  Assert(false, ExcNotImplemented());
1516  break;
1517  }
1519  {
1520  static constexpr ndarray<unsigned int, 3, 2> subcells = {
1521  {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
1522 
1523  return subcells[face][subface];
1524  }
1526  {
1527  const auto [face_orientation, face_rotation, face_flip] =
1529 
1532  face,
1533  subface,
1534  face_orientation,
1535  face_flip,
1536  face_rotation);
1537  }
1540  case ReferenceCells::Wedge:
1541  {
1542  Assert(false, ExcNotImplemented());
1543  break;
1544  }
1546  {
1547  const auto [face_orientation, face_rotation, face_flip] =
1549 
1552  face,
1553  subface,
1554  face_orientation,
1555  face_flip,
1556  face_rotation);
1557  }
1558  default:
1559  Assert(false, ExcNotImplemented());
1560  }
1561 
1563 }
1564 
1565 
1566 
1567 inline std::array<unsigned int, 2>
1569  const unsigned int vertex) const
1570 {
1572  // Work around a GCC warning at higher optimization levels by making all of
1573  // these tables the same size
1574  constexpr unsigned int X = numbers::invalid_unsigned_int;
1575 
1576  switch (this->kind)
1577  {
1579  case ReferenceCells::Line:
1580  Assert(false, ExcNotImplemented());
1581  break;
1583  {
1584  static constexpr ndarray<unsigned int, 6, 2> table = {
1585  {{{0, 0}}, {{0, 1}}, {{1, 1}}, {{X, X}}, {{X, X}}, {{X, X}}}};
1586 
1587  return table[vertex];
1588  }
1590  {
1592  vertex);
1593  }
1595  {
1596  static constexpr ndarray<unsigned int, 6, 2> table = {
1597  {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 2}}, {{X, X}}, {{X, X}}}};
1598 
1599  return table[vertex];
1600  }
1602  {
1603  static constexpr ndarray<unsigned int, 6, 2> table = {
1604  {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{0, 3}}, {{1, 2}}, {{X, X}}}};
1605 
1606  return table[vertex];
1607  }
1608  case ReferenceCells::Wedge:
1609  {
1610  static constexpr ndarray<unsigned int, 6, 2> table = {
1611  {{{0, 1}}, {{0, 0}}, {{0, 2}}, {{1, 0}}, {{1, 1}}, {{1, 2}}}};
1612 
1613  return table[vertex];
1614  }
1616  {
1618  vertex);
1619  }
1620  default:
1621  Assert(false, ExcNotImplemented());
1622  }
1623 
1624  return {};
1625 }
1626 
1627 
1628 
1629 inline std::array<unsigned int, 2>
1631  const unsigned int line) const
1632 {
1633  AssertIndexRange(line, n_lines());
1634 
1635  switch (this->kind)
1636  {
1638  case ReferenceCells::Line:
1641  {
1642  Assert(false, ExcNotImplemented());
1643  break;
1644  }
1646  {
1647  static const std::array<unsigned int, 2> table[6] = {
1648  {{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}};
1649 
1650  return table[line];
1651  }
1653  {
1654  static const std::array<unsigned int, 2> table[8] = {{{0, 0}},
1655  {{0, 1}},
1656  {{0, 2}},
1657  {{0, 3}},
1658  {{1, 2}},
1659  {{2, 1}},
1660  {{1, 1}},
1661  {{2, 2}}};
1662 
1663  return table[line];
1664  }
1665  case ReferenceCells::Wedge:
1666  {
1667  static const std::array<unsigned int, 2> table[9] = {{{0, 0}},
1668  {{0, 2}},
1669  {{0, 1}},
1670  {{1, 0}},
1671  {{1, 1}},
1672  {{1, 2}},
1673  {{2, 0}},
1674  {{2, 1}},
1675  {{3, 1}}};
1676 
1677  return table[line];
1678  }
1680  {
1682  }
1683  default:
1684  Assert(false, ExcNotImplemented());
1685  }
1686 
1687  return {};
1688 }
1689 
1690 
1691 inline unsigned int
1693  const unsigned int vertex) const
1694 {
1696  AssertIndexRange(line, n_lines());
1697 
1698  switch (this->kind)
1699  {
1701  case ReferenceCells::Line:
1702  return vertex;
1704  {
1705  static constexpr ndarray<unsigned int, 3, 2> table = {
1706  {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
1707  return table[line][vertex];
1708  }
1710  {
1711  static constexpr ndarray<unsigned int, 4, 2> table = {
1712  {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
1713  return table[line][vertex];
1714  }
1716  {
1717  static constexpr ndarray<unsigned int, 6, 2> table = {
1718  {{{0, 1}}, {{1, 2}}, {{2, 0}}, {{0, 3}}, {{1, 3}}, {{2, 3}}}};
1719  return table[line][vertex];
1720  }
1722  {
1723  static constexpr ndarray<unsigned int, 8, 2> table = {{{{0, 2}},
1724  {{1, 3}},
1725  {{0, 1}},
1726  {{2, 3}},
1727  {{4, 0}},
1728  {{1, 4}},
1729  {{2, 4}},
1730  {{4, 3}}}};
1731  return table[line][vertex];
1732  }
1733  case ReferenceCells::Wedge:
1734  {
1735  static constexpr ndarray<unsigned int, 9, 2> table = {{{{1, 0}},
1736  {{2, 1}},
1737  {{0, 2}},
1738  {{3, 4}},
1739  {{4, 5}},
1740  {{5, 3}},
1741  {{0, 3}},
1742  {{1, 4}},
1743  {{2, 5}}}};
1744  return table[line][vertex];
1745  }
1747  {
1748  // first four lines comprise the bottom face, next four are the top,
1749  // and the last four are 'bottom to top'
1750  static constexpr ndarray<unsigned int, 12, 2> table = {{{{0, 2}},
1751  {{1, 3}},
1752  {{0, 1}},
1753  {{2, 3}},
1754  {{4, 6}},
1755  {{5, 7}},
1756  {{4, 5}},
1757  {{6, 7}},
1758  {{0, 4}},
1759  {{1, 5}},
1760  {{2, 6}},
1761  {{3, 7}}}};
1762  return table[line][vertex];
1763  }
1764 
1765  default:
1766  Assert(false, ExcNotImplemented());
1767  }
1768 
1770 }
1771 
1772 
1773 inline unsigned int
1775  const unsigned int face,
1776  const unsigned int line,
1777  const unsigned char combined_face_orientation) const
1778 {
1779  AssertIndexRange(face, n_faces());
1781 
1782  static constexpr unsigned int X = numbers::invalid_unsigned_int;
1783 
1784  switch (this->kind)
1785  {
1787  {
1788  Assert(false, ExcNotImplemented());
1789  break;
1790  }
1791  case ReferenceCells::Line:
1792  {
1793  const auto [face_orientation, face_rotation, face_flip] =
1795 
1797  face, line, face_orientation, face_flip, face_rotation);
1798  }
1800  {
1801  return face;
1802  }
1804  {
1805  const auto [face_orientation, face_rotation, face_flip] =
1807 
1809  face, line, face_orientation, face_flip, face_rotation);
1810  }
1812  {
1813  static constexpr ndarray<unsigned int, 4, 3> table = {
1814  {{{0, 1, 2}}, {{0, 3, 4}}, {{2, 5, 3}}, {{1, 4, 5}}}};
1815 
1816  return table[face][standard_to_real_face_line(
1817  line, face, combined_face_orientation)];
1818  }
1820  {
1821  static constexpr ndarray<unsigned int, 5, 4> table = {
1822  {{{0, 1, 2, 3}},
1823  {{0, 6, 4, X}},
1824  {{1, 5, 7, X}},
1825  {{2, 4, 5, X}},
1826  {{3, 7, 6, X}}}};
1827 
1828  return table[face][standard_to_real_face_line(
1829  line, face, combined_face_orientation)];
1830  }
1831  case ReferenceCells::Wedge:
1832  {
1833  static constexpr ndarray<unsigned int, 5, 4> table = {
1834  {{{0, 2, 1, X}},
1835  {{3, 4, 5, X}},
1836  {{6, 7, 0, 3}},
1837  {{7, 8, 1, 4}},
1838  {{8, 6, 5, 2}}}};
1839 
1840  return table[face][standard_to_real_face_line(
1841  line, face, combined_face_orientation)];
1842  }
1844  {
1845  const auto [face_orientation, face_rotation, face_flip] =
1847 
1849  face, line, face_orientation, face_flip, face_rotation);
1850  }
1851  default:
1852  Assert(false, ExcNotImplemented());
1853  }
1854 
1856 }
1857 
1858 
1859 
1860 inline unsigned int
1862  const unsigned int face,
1863  const unsigned int vertex,
1864  const unsigned char combined_face_orientation) const
1865 {
1866  AssertIndexRange(face, n_faces());
1868 
1869  switch (this->kind)
1870  {
1872  {
1873  Assert(false, ExcNotImplemented());
1874  break;
1875  }
1876  case ReferenceCells::Line:
1877  {
1878  const auto [face_orientation, face_rotation, face_flip] =
1880 
1882  face, vertex, face_orientation, face_flip, face_rotation);
1883  }
1885  {
1886  static constexpr ndarray<unsigned int, 3, 2> table = {
1887  {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
1888 
1889  return table[face][combined_face_orientation !=
1891  vertex :
1892  (1 - vertex)];
1893  }
1895  {
1896  const auto [face_orientation, face_rotation, face_flip] =
1898 
1900  face, vertex, face_orientation, face_flip, face_rotation);
1901  }
1903  {
1904  static constexpr ndarray<unsigned int, 4, 3> table = {
1905  {{{0, 1, 2}}, {{1, 0, 3}}, {{0, 2, 3}}, {{2, 1, 3}}}};
1906 
1907  return table[face][standard_to_real_face_vertex(
1909  }
1911  {
1912  constexpr auto X = numbers::invalid_unsigned_int;
1913  static constexpr ndarray<unsigned int, 5, 4> table = {
1914  {{{0, 1, 2, 3}},
1915  {{0, 2, 4, X}},
1916  {{3, 1, 4, X}},
1917  {{1, 0, 4, X}},
1918  {{2, 3, 4, X}}}};
1919 
1920  return table[face][standard_to_real_face_vertex(
1922  }
1923  case ReferenceCells::Wedge:
1924  {
1925  constexpr auto X = numbers::invalid_unsigned_int;
1926  static constexpr ndarray<unsigned int, 6, 4> table = {
1927  {{{1, 0, 2, X}},
1928  {{3, 4, 5, X}},
1929  {{0, 1, 3, 4}},
1930  {{1, 2, 4, 5}},
1931  {{2, 0, 5, 3}}}};
1932 
1933  return table[face][standard_to_real_face_vertex(
1935  }
1937  {
1938  const auto [face_orientation, face_rotation, face_flip] =
1940 
1942  face, vertex, face_orientation, face_flip, face_rotation);
1943  }
1944  default:
1945  Assert(false, ExcNotImplemented());
1946  }
1947 
1949 }
1950 
1951 
1952 
1953 template <int dim>
1954 Point<dim>
1955 ReferenceCell::face_vertex_location(const unsigned int face,
1956  const unsigned int vertex) const
1957 {
1958  return this->template vertex<dim>(
1960 }
1961 
1962 
1963 
1964 inline unsigned int
1966  const unsigned int vertex,
1967  const unsigned int face,
1968  const unsigned char face_orientation) const
1969 {
1970  AssertIndexRange(face, n_faces());
1972 
1973  switch (this->kind)
1974  {
1976  case ReferenceCells::Line:
1977  Assert(false, ExcNotImplemented());
1978  break;
1981  return line_vertex_permutations[face_orientation][vertex];
1983  return triangle_vertex_permutations[face_orientation][vertex];
1985  // face 0 is a quadrilateral
1986  if (face == 0)
1987  return quadrilateral_vertex_permutations[face_orientation][vertex];
1988  else
1989  return triangle_vertex_permutations[face_orientation][vertex];
1990  case ReferenceCells::Wedge:
1991  // faces 0 and 1 are triangles
1992  if (face > 1)
1993  return quadrilateral_vertex_permutations[face_orientation][vertex];
1994  else
1995  return triangle_vertex_permutations[face_orientation][vertex];
1997  return quadrilateral_vertex_permutations[face_orientation][vertex];
1998  default:
1999  Assert(false, ExcNotImplemented());
2000  }
2001 
2002  Assert(false, ExcNotImplemented());
2004 }
2005 
2006 
2007 
2008 inline unsigned int
2010  const unsigned int line,
2011  const unsigned int face,
2012  const unsigned char combined_face_orientation) const
2013 {
2014  AssertIndexRange(face, n_faces());
2016 
2017  static constexpr ndarray<unsigned int, 6, 3> triangle_table = {{{{2, 1, 0}},
2018  {{0, 1, 2}},
2019  {{1, 0, 2}},
2020  {{2, 0, 1}},
2021  {{0, 2, 1}},
2022  {{1, 2, 0}}}};
2023 
2024 
2025  switch (this->kind)
2026  {
2028  case ReferenceCells::Line:
2031  Assert(false, ExcNotImplemented());
2032  break;
2034  return triangle_table[combined_face_orientation][line];
2036  if (face == 0) // The quadrilateral face
2037  {
2038  const auto [face_orientation, face_rotation, face_flip] =
2040 
2042  face_orientation,
2043  face_flip,
2044  face_rotation);
2045  }
2046  else // One of the triangular faces
2047  {
2048  return triangle_table[combined_face_orientation][line];
2049  }
2050  case ReferenceCells::Wedge:
2051  if (face > 1) // One of the quadrilateral faces
2052  {
2053  const auto [face_orientation, face_rotation, face_flip] =
2055 
2057  face_orientation,
2058  face_flip,
2059  face_rotation);
2060  }
2061  else // One of the triangular faces
2062  return triangle_table[combined_face_orientation][line];
2064  {
2065  static constexpr ndarray<unsigned int, 8, 4> table = {
2066  {{{2, 3, 0, 1}},
2067  {{0, 1, 2, 3}},
2068  {{0, 1, 3, 2}},
2069  {{3, 2, 0, 1}},
2070  {{3, 2, 1, 0}},
2071  {{1, 0, 3, 2}},
2072  {{1, 0, 2, 3}},
2073  {{2, 3, 1, 0}}}};
2074  return table[combined_face_orientation][line];
2075  }
2076  default:
2077  Assert(false, ExcNotImplemented());
2078  }
2079 
2081 }
2082 
2083 
2084 
2085 namespace ReferenceCells
2086 {
2087  template <int dim>
2088  inline constexpr const ReferenceCell &
2090  {
2091  switch (dim)
2092  {
2093  case 0:
2094  return ReferenceCells::Vertex;
2095  case 1:
2096  return ReferenceCells::Line;
2097  case 2:
2098  return ReferenceCells::Triangle;
2099  case 3:
2101  default:
2102  Assert(false, ExcNotImplemented());
2103  }
2104  return ReferenceCells::Invalid;
2105  }
2106 
2107 
2108 
2109  template <int dim>
2110  inline constexpr const ReferenceCell &
2112  {
2113  switch (dim)
2114  {
2115  case 0:
2116  return ReferenceCells::Vertex;
2117  case 1:
2118  return ReferenceCells::Line;
2119  case 2:
2121  case 3:
2123  default:
2124  Assert(false, ExcNotImplemented());
2125  }
2126  return ReferenceCells::Invalid;
2127  }
2128 } // namespace ReferenceCells
2129 
2130 
2131 inline ReferenceCell
2132 ReferenceCell::n_vertices_to_type(const int dim, const unsigned int n_vertices)
2133 {
2134  AssertIndexRange(dim, 4);
2136 
2137  const auto X = ReferenceCells::Invalid;
2138  static const ndarray<ReferenceCell, 4, 9> table = {
2139  {// dim 0
2140  {{X, ReferenceCells::Vertex, X, X, X, X, X, X, X}},
2141  // dim 1
2142  {{X, X, ReferenceCells::Line, X, X, X, X, X, X}},
2143  // dim 2
2144  {{X,
2145  X,
2146  X,
2149  X,
2150  X,
2151  X,
2152  X}},
2153  // dim 3
2154  {{X,
2155  X,
2156  X,
2157  X,
2161  X,
2163  Assert(table[dim][n_vertices] != ReferenceCells::Invalid,
2164  ExcMessage("The combination of dim = " + std::to_string(dim) +
2165  " and n_vertices = " + std::to_string(n_vertices) +
2166  " does not correspond to a known reference cell type."));
2167  return table[dim][n_vertices];
2168 }
2169 
2170 
2171 
2172 template <int dim>
2173 inline double
2175  const unsigned int i) const
2176 {
2179  switch (this->kind)
2180  {
2182  case ReferenceCells::Line:
2186  // see also BarycentricPolynomials<2>::compute_value
2188  {
2189  switch (i)
2190  {
2191  case 0:
2192  return 1.0 - xi[std::min(0, dim - 1)] -
2193  xi[std::min(1, dim - 1)];
2194  case 1:
2195  return xi[std::min(0, dim - 1)];
2196  case 2:
2197  return xi[std::min(1, dim - 1)];
2198  default:
2199  Assert(false, ExcInternalError());
2200  }
2201  }
2202  // see also BarycentricPolynomials<3>::compute_value
2204  {
2205  switch (i)
2206  {
2207  case 0:
2208  return 1.0 - xi[std::min(0, dim - 1)] -
2209  xi[std::min(1, dim - 1)] - xi[std::min(2, dim - 1)];
2210  case 1:
2211  return xi[std::min(0, dim - 1)];
2212  case 2:
2213  return xi[std::min(1, dim - 1)];
2214  case 3:
2215  return xi[std::min(2, dim - 1)];
2216  default:
2217  Assert(false, ExcInternalError());
2218  }
2219  }
2220  // see also ScalarLagrangePolynomialPyramid::compute_value()
2222  {
2223  const double Q14 = 0.25;
2224 
2225  const double r = xi[std::min(0, dim - 1)];
2226  const double s = xi[std::min(1, dim - 1)];
2227  const double t = xi[std::min(2, dim - 1)];
2228 
2229  const double ratio =
2230  (std::fabs(t - 1.0) > 1.0e-14 ? (r * s * t) / (1.0 - t) : 0.0);
2231 
2232  if (i == 0)
2233  return Q14 * ((1.0 - r) * (1.0 - s) - t + ratio);
2234  if (i == 1)
2235  return Q14 * ((1.0 + r) * (1.0 - s) - t - ratio);
2236  if (i == 2)
2237  return Q14 * ((1.0 - r) * (1.0 + s) - t - ratio);
2238  if (i == 3)
2239  return Q14 * ((1.0 + r) * (1.0 + s) - t + ratio);
2240  else
2241  return t;
2242  }
2243  // see also ScalarLagrangePolynomialWedge::compute_value()
2244  case ReferenceCells::Wedge:
2246  .d_linear_shape_function<2>(Point<2>(xi[std::min(0, dim - 1)],
2247  xi[std::min(1, dim - 1)]),
2248  i % 3) *
2250  .d_linear_shape_function<1>(Point<1>(xi[std::min(2, dim - 1)]),
2251  i / 3);
2252  default:
2253  Assert(false, ExcNotImplemented());
2254  }
2255 
2256  return 0.0;
2257 }
2258 
2259 
2260 
2261 template <int dim>
2262 inline Tensor<1, dim>
2264  const unsigned int i) const
2265 {
2267  switch (this->kind)
2268  {
2270  case ReferenceCells::Line:
2274  // see also BarycentricPolynomials<2>::compute_grad()
2276  switch (i)
2277  {
2278  case 0:
2279  return Point<dim>(-1.0, -1.0);
2280  case 1:
2281  return Point<dim>(+1.0, +0.0);
2282  case 2:
2283  return Point<dim>(+0.0, +1.0);
2284  default:
2285  Assert(false, ExcInternalError());
2286  }
2287  default:
2288  Assert(false, ExcNotImplemented());
2289  }
2290 
2291  return Point<dim>(+0.0, +0.0, +0.0);
2292 }
2293 
2294 
2295 
2296 inline double
2298 {
2299  switch (this->kind)
2300  {
2302  return 0;
2303  case ReferenceCells::Line:
2304  return 1;
2306  return 1. / 2.;
2308  return 1;
2310  return 1. / 6.;
2312  return 4. / 3.;
2313  case ReferenceCells::Wedge:
2314  return 1. / 2.;
2316  return 1;
2317  default:
2318  Assert(false, ExcNotImplemented());
2319  }
2320 
2321  return 0.0;
2322 }
2323 
2324 
2325 
2326 template <int dim>
2327 inline Point<dim>
2329 {
2331 
2332  switch (this->kind)
2333  {
2335  return Point<dim>();
2336  case ReferenceCells::Line:
2337  return Point<dim>(1. / 2.);
2339  return Point<dim>(1. / 3., 1. / 3.);
2341  return Point<dim>(1. / 2., 1. / 2.);
2343  return Point<dim>(1. / 4., 1. / 4., 1. / 4.);
2345  return Point<dim>(0, 0, 1. / 4.);
2346  case ReferenceCells::Wedge:
2347  return Point<dim>(1. / 3, 1. / 3, 1. / 2.);
2349  return Point<dim>(1. / 2., 1. / 2., 1. / 2.);
2350  default:
2351  Assert(false, ExcNotImplemented());
2352  }
2353 
2354  return Point<dim>();
2355 }
2356 
2357 
2358 
2359 template <int dim>
2360 inline bool
2361 ReferenceCell::contains_point(const Point<dim> &p, const double tolerance) const
2362 {
2364 
2365  // Introduce abbreviations to silence compiler warnings about invalid
2366  // array accesses (that can't happen, of course, but the compiler
2367  // doesn't know that).
2368  constexpr unsigned int x_coordinate = 0;
2369  constexpr unsigned int y_coordinate = (dim >= 2 ? 1 : 0);
2370  constexpr unsigned int z_coordinate = (dim >= 3 ? 2 : 0);
2371 
2372  switch (this->kind)
2373  {
2375  {
2376  // Vertices are special cases in that they do not actually
2377  // have coordinates. Error out if this function is called
2378  // with a vertex:
2379  Assert(false,
2380  ExcMessage("Vertices are zero-dimensional objects and "
2381  "as a consequence have no coordinates. You "
2382  "cannot meaningfully ask whether a point is "
2383  "inside a vertex (within a certain tolerance) "
2384  "without coordinate values."));
2385  return false;
2386  }
2387  case ReferenceCells::Line:
2390  {
2391  for (unsigned int d = 0; d < dim; ++d)
2392  if ((p[d] < -tolerance) || (p[d] > 1 + tolerance))
2393  return false;
2394  return true;
2395  }
2398  {
2399  // First make sure that we are in the first quadrant or octant
2400  for (unsigned int d = 0; d < dim; ++d)
2401  if (p[d] < -tolerance)
2402  return false;
2403 
2404  // Now we also need to make sure that we are below the diagonal line
2405  // or plane that delineates the simplex. This diagonal is given by
2406  // sum(p[d])<=1, and a diagonal a distance eps away is given by
2407  // sum(p[d])<=1+eps*sqrt(d). (For example, the point at (1,1) is a
2408  // distance of 1/sqrt(2) away from the diagonal. That is, its
2409  // sum satisfies
2410  // sum(p[d]) = 2 <= 1 + (1/sqrt(2)) * sqrt(2)
2411  // in other words, it satisfies the predicate with eps=1/sqrt(2).)
2412  double sum = 0;
2413  for (unsigned int d = 0; d < dim; ++d)
2414  sum += p[d];
2415  return (sum <= 1 + tolerance * std::sqrt(1. * dim));
2416  }
2418  {
2419  // A pyramid only lives in the upper half-space:
2420  if (p[z_coordinate] < -tolerance)
2421  return false;
2422 
2423  // It also only lives in the space below z=1:
2424  if (p[z_coordinate] > 1 + tolerance)
2425  return false;
2426 
2427  // Within what's left of the space, a pyramid is a cone that tapers
2428  // towards the top. First compute the distance of the point to the
2429  // axis in the max norm (this is the right norm because the vertices
2430  // of the pyramid are at points +/-1, +/-1):
2431  const double distance_from_axis =
2432  std::max(std::fabs(p[x_coordinate]), std::fabs(p[y_coordinate]));
2433 
2434  // We are inside the pyramid if the distance from the axis is less
2435  // than (1-z)
2436  return (distance_from_axis <= 1 + tolerance - p[z_coordinate]);
2437  }
2438  case ReferenceCells::Wedge:
2439  {
2440  // The wedge we use is a triangle extruded into the third
2441  // dimension by one unit. So we can use the same logic as for
2442  // triangles above (i.e., for the simplex above, using dim==2)
2443  // and then check the third dimension separately.
2444 
2445  if ((p[x_coordinate] < -tolerance) || (p[y_coordinate] < -tolerance))
2446  return false;
2447 
2448  const double sum = p[x_coordinate] + p[y_coordinate];
2449  if (sum > 1 + tolerance * std::sqrt(2.0))
2450  return false;
2451 
2452  if (p[z_coordinate] < -tolerance)
2453  return false;
2454  if (p[z_coordinate] > 1 + tolerance)
2455  return false;
2456 
2457  return true;
2458  }
2459  default:
2460  Assert(false, ExcNotImplemented());
2461  }
2462 
2463  return false;
2464 }
2465 
2466 
2467 
2468 template <int dim>
2469 inline Tensor<1, dim>
2470 ReferenceCell::unit_tangential_vectors(const unsigned int face_no,
2471  const unsigned int i) const
2472 {
2474  AssertIndexRange(i, dim - 1);
2475 
2476  switch (this->kind)
2477  {
2479  case ReferenceCells::Line:
2483  return GeometryInfo<dim>::unit_tangential_vectors[face_no][i];
2485  {
2486  AssertIndexRange(face_no, 3);
2487  static const std::array<Tensor<1, dim>, 3> table = {
2488  {Point<dim>(1, 0),
2489  Point<dim>(-std::sqrt(0.5), +std::sqrt(0.5)),
2490  Point<dim>(0, -1)}};
2491 
2492  return table[face_no];
2493  }
2495  {
2496  AssertIndexRange(face_no, 4);
2497  static const ndarray<Tensor<1, dim>, 4, 2> table = {
2498  {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
2499  {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
2500  {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
2501  {{Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
2502  +std::pow(1.0 / 3.0, 1.0 / 4.0),
2503  0),
2504  Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
2505  0,
2506  +std::pow(1.0 / 3.0, 1.0 / 4.0))}}}};
2507 
2508  return table[face_no][i];
2509  }
2511  {
2512  AssertIndexRange(face_no, 5);
2513  static const ndarray<Tensor<1, dim>, 5, 2> table = {
2514  {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
2515  {{Point<dim>(+1.0 / sqrt(2.0), 0, +1.0 / sqrt(2.0)),
2516  Point<dim>(0, 1, 0)}},
2517  {{Point<dim>(+1.0 / sqrt(2.0), 0, -1.0 / sqrt(2.0)),
2518  Point<dim>(0, 1, 0)}},
2519  {{Point<dim>(1, 0, 0),
2520  Point<dim>(0, +1.0 / sqrt(2.0), +1.0 / sqrt(2.0))}},
2521  {{Point<dim>(1, 0, 0),
2522  Point<dim>(0, +1.0 / sqrt(2.0), -1.0 / sqrt(2.0))}}}};
2523 
2524  return table[face_no][i];
2525  }
2526  case ReferenceCells::Wedge:
2527  {
2528  AssertIndexRange(face_no, 5);
2529  static const ndarray<Tensor<1, dim>, 5, 2> table = {
2530  {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
2531  {{Point<dim>(1, 0, 0), Point<dim>(0, 1, 0)}},
2532  {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
2533  {{Point<dim>(-1 / std::sqrt(2.0), +1 / std::sqrt(2.0), 0),
2534  Point<dim>(0, 0, 1)}},
2535  {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}}}};
2536 
2537  return table[face_no][i];
2538  }
2539  default:
2540  Assert(false, ExcNotImplemented());
2541  }
2542 
2543 
2544  return {};
2545 }
2546 
2547 
2548 
2549 template <int dim>
2550 inline Tensor<1, dim>
2551 ReferenceCell::unit_normal_vectors(const unsigned int face_no) const
2552 {
2553  AssertDimension(dim, this->get_dimension());
2554 
2555  if (is_hyper_cube())
2556  {
2558  return GeometryInfo<dim>::unit_normal_vector[face_no];
2559  }
2560  else if (dim == 2)
2561  {
2563 
2564  // Return the rotated vector
2565  return cross_product_2d(unit_tangential_vectors<dim>(face_no, 0));
2566  }
2567  else if (dim == 3)
2568  {
2569  return cross_product_3d(unit_tangential_vectors<dim>(face_no, 0),
2570  unit_tangential_vectors<dim>(face_no, 1));
2571  }
2572 
2573  Assert(false, ExcNotImplemented());
2574 
2575  return {};
2576 }
2577 
2578 
2579 
2580 inline unsigned int
2581 ReferenceCell::n_face_orientations(const unsigned int face_no) const
2582 {
2583  AssertIndexRange(face_no, n_faces());
2584  if (get_dimension() == 1)
2585  return 1;
2586  if (get_dimension() == 2)
2587  return 2;
2589  return 8;
2590  else if (face_reference_cell(face_no) == ReferenceCells::Triangle)
2591  return 6;
2592 
2593  Assert(false, ExcInternalError());
2595 }
2596 
2597 
2598 
2599 inline bool
2601  const unsigned int line,
2602  const unsigned int face,
2603  const unsigned char combined_face_orientation,
2604  const bool line_orientation) const
2605 {
2606  if (*this == ReferenceCells::Hexahedron)
2607  {
2608  static constexpr ::ndarray<bool, 2, 8> bool_table{
2609  {{{true, true, false, true, false, false, true, false}},
2610  {{true, true, true, false, false, false, false, true}}}};
2611 
2612  return (line_orientation ==
2613  bool_table[line / 2][combined_face_orientation]);
2614  }
2615  else if (*this == ReferenceCells::Tetrahedron)
2616  {
2617  static constexpr unsigned int X = numbers::invalid_unsigned_int;
2618  static constexpr ::ndarray<unsigned int, 4, 3> combined_lines{
2619  {{{0, 0, 0}}, {{X, 0, 1}}, {{X, 0, X}}, {{X, X, X}}}};
2620 
2621  const auto combined_line = combined_lines[face][line];
2622 
2623  Assert(combined_line != X,
2624  ExcMessage(
2625  "This function can only be called for following face-line "
2626  "combinations: (0,0), (0,1), (0,2), (1,1), (1,2), (2,1),"));
2627 
2628  static constexpr ::ndarray<bool, 2, 6> bool_table{
2629  {{{false, true, false, true, false, true}},
2630  {{true, false, true, false, true, false}}}};
2631 
2632  return (line_orientation ==
2633  bool_table[combined_line][combined_face_orientation]);
2634  }
2635  else
2636  // TODO: This might actually be wrong for some of the other
2637  // kinds of objects. We should check this
2638  return true;
2639 }
2640 
2641 
2642 
2643 namespace internal
2644 {
2645  template <typename T>
2647  {
2648  public:
2658  {
2661  }
2662 
2666  virtual ~NoPermutation() noexcept override = default;
2667 
2671  virtual void
2672  print_info(std::ostream &out) const override
2673  {
2674  out << '[';
2675 
2676  const unsigned int n_vertices = entity_type.n_vertices();
2677 
2678  for (unsigned int i = 0; i < n_vertices; ++i)
2679  {
2680  out << vertices_0[i];
2681  if (i + 1 != n_vertices)
2682  out << ',';
2683  }
2684 
2685  out << "] is not a valid permutation of [";
2686 
2687  for (unsigned int i = 0; i < n_vertices; ++i)
2688  {
2689  out << vertices_1[i];
2690  if (i + 1 != n_vertices)
2691  out << ',';
2692  }
2693 
2694  out << "]." << std::endl;
2695  }
2696 
2701 
2706 
2711  };
2712 
2722  ReferenceCell,
2723  ReferenceCell,
2724  << "The reference-cell type used on this cell (" << arg1.to_string()
2725  << ") does not match the reference-cell type of the finite element "
2726  << "associated with this cell (" << arg2.to_string() << "). "
2727  << "Did you accidentally use simplex elements on hypercube meshes "
2728  << "(or the other way around), or are you using a mixed mesh and "
2729  << "assigned a simplex element to a hypercube cell (or the other "
2730  << "way around) via the active_fe_index?");
2731 } // namespace internal
2732 
2733 
2734 
2735 template <typename T, std::size_t N>
2736 inline unsigned char
2737 ReferenceCell::compute_orientation(const std::array<T, N> &vertices_0,
2738  const std::array<T, N> &vertices_1) const
2739 {
2740  Assert(N >= n_vertices(),
2741  ExcMessage("The number of array elements must be equal to or "
2742  "greater than the number of vertices of the cell "
2743  "referenced by this object."));
2744 
2745  // Call the non-deprecated function, taking care of calling it only with
2746  // those array elements that we actually care about (see the note
2747  // in the documentation about the arguments potentially being
2748  // larger arrays than necessary).
2749  return get_combined_orientation(
2750  make_array_view(vertices_0.begin(), vertices_0.begin() + n_vertices()),
2751  make_array_view(vertices_1.begin(), vertices_1.begin() + n_vertices()));
2752 }
2753 
2754 
2755 
2756 template <typename T>
2757 unsigned char
2759  const ArrayView<const T> &vertices_0,
2760  const ArrayView<const T> &vertices_1) const
2761 {
2762  Assert(vertices_0.size() == n_vertices(),
2763  ExcMessage("The number of array elements must be equal to "
2764  "the number of vertices of the cell "
2765  "referenced by this object."));
2766  Assert(vertices_1.size() == n_vertices(),
2767  ExcMessage("The number of array elements must be equal to "
2768  "the number of vertices of the cell "
2769  "referenced by this object."));
2770 
2771  const auto v0_equals = [&](const std::initializer_list<const T> &list) {
2772  Assert(list.size() == n_vertices(), ExcInternalError());
2773  return std::equal(vertices_0.begin(), vertices_0.end(), std::begin(list));
2774  };
2775 
2776  switch (this->kind)
2777  {
2779  // Things are always default-oriented in 1D
2780  if (v0_equals({vertices_1[0]}))
2782  break;
2783 
2784  case ReferenceCells::Line:
2785  // line_orientation=true
2786  if (v0_equals({vertices_1[0], vertices_1[1]}))
2788 
2789  // line_orientation=false
2790  if (v0_equals({vertices_1[1], vertices_1[0]}))
2792  break;
2794  // face_orientation=true, face_rotation=false, face_flip=false
2795  if (v0_equals({vertices_1[0], vertices_1[1], vertices_1[2]}))
2796  return 1;
2797 
2798  // face_orientation=true, face_rotation=true, face_flip=false
2799  if (v0_equals({vertices_1[1], vertices_1[2], vertices_1[0]}))
2800  return 5;
2801 
2802  // face_orientation=true, face_rotation=false, face_flip=true
2803  if (v0_equals({vertices_1[2], vertices_1[0], vertices_1[1]}))
2804  return 3;
2805 
2806  // face_orientation=false, face_rotation=false, face_flip=false
2807  if (v0_equals({vertices_1[0], vertices_1[2], vertices_1[1]}))
2808  return 0;
2809 
2810  // face_orientation=false, face_rotation=true, face_flip=false
2811  if (v0_equals({vertices_1[2], vertices_1[1], vertices_1[0]}))
2812  return 2;
2813 
2814  // face_orientation=false, face_rotation=false, face_flip=true
2815  if (v0_equals({vertices_1[1], vertices_1[0], vertices_1[2]}))
2816  return 4;
2817  break;
2819  // face_orientation=true, face_rotation=false, face_flip=false
2820  if (v0_equals(
2821  {vertices_1[0], vertices_1[1], vertices_1[2], vertices_1[3]}))
2822  return 1;
2823 
2824  // face_orientation=true, face_rotation=true, face_flip=false
2825  if (v0_equals(
2826  {vertices_1[2], vertices_1[0], vertices_1[3], vertices_1[1]}))
2827  return 3;
2828 
2829  // face_orientation=true, face_rotation=false, face_flip=true
2830  if (v0_equals(
2831  {vertices_1[3], vertices_1[2], vertices_1[1], vertices_1[0]}))
2832  return 5;
2833 
2834  // face_orientation=true, face_rotation=true, face_flip=true
2835  if (v0_equals(
2836  {vertices_1[1], vertices_1[3], vertices_1[0], vertices_1[2]}))
2837  return 7;
2838 
2839  // face_orientation=false, face_rotation=false, face_flip=false
2840  if (v0_equals(
2841  {vertices_1[0], vertices_1[2], vertices_1[1], vertices_1[3]}))
2842  return 0;
2843 
2844  // face_orientation=false, face_rotation=true, face_flip=false
2845  if (v0_equals(
2846  {vertices_1[2], vertices_1[3], vertices_1[0], vertices_1[1]}))
2847  return 2;
2848 
2849  // face_orientation=false, face_rotation=false, face_flip=true
2850  if (v0_equals(
2851  {vertices_1[3], vertices_1[1], vertices_1[2], vertices_1[0]}))
2852  return 4;
2853 
2854  // face_orientation=false, face_rotation=true, face_flip=true
2855  if (v0_equals(
2856  {vertices_1[1], vertices_1[0], vertices_1[3], vertices_1[2]}))
2857  return 6;
2858  break;
2859  default:
2860  Assert(false, ExcNotImplemented());
2861  }
2862 
2863  Assert(false, (internal::NoPermutation<T>(*this, vertices_0, vertices_1)));
2865 }
2866 
2867 
2868 
2869 template <typename T, std::size_t N>
2870 inline std::array<T, N>
2872  const std::array<T, N> &vertices,
2873  const unsigned int orientation) const
2874 {
2875  Assert(N >= n_vertices(),
2876  ExcMessage("The number of array elements must be equal to or "
2877  "greater than the number of vertices of the cell "
2878  "referenced by this object."));
2879 
2880  // Call the non-deprecated function, taking care of calling it only with
2881  // those array elements that we actually care about (see the note
2882  // in the documentation about the arguments potentially being
2883  // larger arrays than necessary).
2884  const auto permutation = permute_by_combined_orientation(
2885  make_array_view(vertices.begin(), vertices.begin() + n_vertices()),
2886  orientation);
2887 
2888  std::array<T, N> temp;
2889  std::copy(permutation.begin(), permutation.end(), temp.begin());
2890 
2891  return temp;
2892 }
2893 
2894 
2895 
2896 template <typename T>
2897 boost::container::small_vector<T, 8>
2900  const unsigned char orientation) const
2901 {
2902  Assert(vertices.size() == n_vertices(),
2903  ExcMessage("The number of array elements must be equal to "
2904  "the number of vertices of the cell "
2905  "referenced by this object."));
2906 
2907  switch (this->kind)
2908  {
2909  case ReferenceCells::Line:
2910  switch (orientation)
2911  {
2912  case 1:
2913  return {vertices[0], vertices[1]};
2914  case 0:
2915  return {vertices[1], vertices[0]};
2916  default:
2917  Assert(false, ExcNotImplemented());
2918  }
2919  break;
2921  switch (orientation)
2922  {
2923  case 1:
2924  return {vertices[0], vertices[1], vertices[2]};
2925  case 3:
2926  return {vertices[1], vertices[2], vertices[0]};
2927  case 5:
2928  return {vertices[2], vertices[0], vertices[1]};
2929  case 0:
2930  return {vertices[0], vertices[2], vertices[1]};
2931  case 2:
2932  return {vertices[2], vertices[1], vertices[0]};
2933  case 4:
2934  return {vertices[1], vertices[0], vertices[2]};
2935  default:
2936  Assert(false, ExcNotImplemented());
2937  }
2938  break;
2940  switch (orientation)
2941  {
2942  case 1:
2943  return {vertices[0], vertices[1], vertices[2], vertices[3]};
2944  case 3:
2945  return {vertices[2], vertices[0], vertices[3], vertices[1]};
2946  case 5:
2947  return {vertices[3], vertices[2], vertices[1], vertices[0]};
2948  case 7:
2949  return {vertices[1], vertices[3], vertices[0], vertices[2]};
2950  case 0:
2951  return {vertices[0], vertices[2], vertices[1], vertices[3]};
2952  case 2:
2953  return {vertices[2], vertices[3], vertices[0], vertices[1]};
2954  case 4:
2955  return {vertices[3], vertices[1], vertices[2], vertices[0]};
2956  case 6:
2957  return {vertices[1], vertices[0], vertices[3], vertices[2]};
2958  default:
2959  Assert(false, ExcNotImplemented());
2960  }
2961  break;
2962  default:
2963  AssertThrow(false, ExcNotImplemented());
2964  }
2965 
2966  return {};
2967 }
2968 
2969 
2970 
2971 template <>
2972 unsigned int
2973 ReferenceCell::vtk_lexicographic_to_node_index<0>(
2974  const std::array<unsigned, 0> &node_indices,
2975  const std::array<unsigned, 0> &nodes_per_direction,
2976  const bool legacy_format) const;
2977 
2978 template <>
2979 unsigned int
2980 ReferenceCell::vtk_lexicographic_to_node_index<1>(
2981  const std::array<unsigned, 1> &node_indices,
2982  const std::array<unsigned, 1> &nodes_per_direction,
2983  const bool legacy_format) const;
2984 
2985 template <>
2986 unsigned int
2987 ReferenceCell::vtk_lexicographic_to_node_index<2>(
2988  const std::array<unsigned, 2> &node_indices,
2989  const std::array<unsigned, 2> &nodes_per_direction,
2990  const bool legacy_format) const;
2991 
2992 template <>
2993 unsigned int
2994 ReferenceCell::vtk_lexicographic_to_node_index<3>(
2995  const std::array<unsigned, 3> &node_indices,
2996  const std::array<unsigned, 3> &nodes_per_direction,
2997  const bool legacy_format) const;
2998 
3000 
3001 #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:950
OutputOperator< VectorType > & operator<<(OutputOperator< VectorType > &out, unsigned int step)
Definition: operator.h:165
iterator begin() const
Definition: array_view.h:703
iterator end() const
Definition: array_view.h:712
std::size_t size() const
Definition: array_view.h:685
Abstract base class for mapping classes.
Definition: mapping.h:317
Definition: point.h:112
static constexpr Point< dim, Number > unit_vector(const unsigned int i)
std::istream & operator>>(std::istream &in, Point< dim, Number > &p)
Definition: point.h:703
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
friend std::ostream & operator<<(std::ostream &out, const ReferenceCell &reference_cell)
void serialize(Archive &archive, const unsigned int)
friend std::istream & operator>>(std::istream &in, ReferenceCell &reference_cell)
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()
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
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:178
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
Point< 3 > vertices[4]
static ::ExceptionBase & ExcNonMatchingReferenceCellTypes(ReferenceCell arg1, ReferenceCell arg2)
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:540
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
Expression fabs(const Expression &x)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
static const char U
static const char N
std::string to_string(const T &t)
Definition: patterns.h:2391
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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()
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
void copy(const T *begin, const T *end, U *dest)
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)
unsigned char combined_face_orientation(const bool face_orientation, const bool face_rotation, const bool face_flip)
static const unsigned int invalid_unsigned_int
Definition: types.h:221
boost::integer_range< IncrementableType > iota_view
Definition: iota_view.h:46
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108
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 double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
static Tensor< 1, dim > d_linear_shape_function_gradient(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 std::array< unsigned int, 2 > standard_quad_vertex_to_line_vertex_index(const unsigned int vertex)
static std::array< unsigned int, 2 > standard_hex_vertex_to_quad_vertex_index(const unsigned int vertex)