Reference documentation for deal.II version GIT b206511199 2023-01-31 13:40: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 - 2022 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 
28 #include <boost/container/small_vector.hpp>
29 
30 #include <iosfwd>
31 #include <string>
32 
34 
35 // Forward declarations
36 #ifndef DOXYGEN
37 template <int dim, int spacedim>
38 class Mapping;
39 
40 template <int dim>
41 class Quadrature;
42 
43 class ReferenceCell;
44 #endif
45 
46 
47 namespace internal
48 {
61  constexpr ReferenceCell
62  make_reference_cell_from_int(const std::uint8_t kind);
63 } // namespace internal
64 
65 
66 
98 {
99 public:
107  static ReferenceCell
108  n_vertices_to_type(const int dim, const unsigned int n_vertices);
109 
118  constexpr ReferenceCell();
119 
130  bool
131  is_hyper_cube() const;
132 
136  bool
137  is_simplex() const;
138 
143  unsigned int
144  get_dimension() const;
145 
159  template <int dim>
160  double
161  d_linear_shape_function(const Point<dim> &xi, const unsigned int i) const;
162 
167  template <int dim>
170  const unsigned int i) const;
171 
181  template <int dim, int spacedim = dim>
182  std::unique_ptr<Mapping<dim, spacedim>>
183  get_default_mapping(const unsigned int degree) const;
184 
196  template <int dim, int spacedim = dim>
197  const Mapping<dim, spacedim> &
199 
208  template <int dim>
210  get_gauss_type_quadrature(const unsigned n_points_1D) const;
211 
221  template <int dim>
223  get_midpoint_quadrature() const;
224 
238  template <int dim>
239  const Quadrature<dim> &
241 
256  unsigned int
257  n_vertices() const;
258 
264  vertex_indices() const;
265 
273  template <int dim>
274  Point<dim>
275  vertex(const unsigned int v) const;
276 
282  unsigned int
283  n_lines() const;
284 
290  line_indices() const;
291 
297  unsigned int
298  n_faces() const;
299 
305  face_indices() const;
306 
318  unsigned int
319  n_isotropic_children() const;
320 
326  isotropic_child_indices() const;
327 
341  face_reference_cell(const unsigned int face_no) const;
342 
356  static constexpr unsigned char
358 
369  static constexpr unsigned char
371 
413  unsigned int
414  child_cell_on_face(const unsigned int face,
415  const unsigned int subface,
416  const unsigned char face_orientation =
418 
428  std::array<unsigned int, 2>
429  standard_vertex_to_face_and_vertex_index(const unsigned int vertex) const;
430 
440  std::array<unsigned int, 2>
441  standard_line_to_face_and_line_index(const unsigned int line) const;
442 
451  unsigned int
452  line_to_cell_vertices(const unsigned int line,
453  const unsigned int vertex) const;
454 
458  unsigned int
459  face_to_cell_lines(const unsigned int face,
460  const unsigned int line,
461  const unsigned char face_orientation) const;
462 
466  unsigned int
467  face_to_cell_vertices(const unsigned int face,
468  const unsigned int vertex,
469  const unsigned char face_orientation) const;
470 
474  unsigned int
475  standard_to_real_face_vertex(const unsigned int vertex,
476  const unsigned int face,
477  const unsigned char face_orientation) const;
478 
482  unsigned int
483  standard_to_real_face_line(const unsigned int line,
484  const unsigned int face,
485  const unsigned char face_orientation) const;
486 
497  bool
498  standard_vs_true_line_orientation(const unsigned int line,
499  const unsigned char face_orientation,
500  const bool line_orientation) const;
501 
525  double
526  volume() const;
527 
540  template <int dim>
541  Point<dim>
542  barycenter() const;
543 
563  template <int dim>
564  bool
565  contains_point(const Point<dim> &p, const double tolerance = 0) const;
566 
567  /*
568  * Return @f$i@f$-th unit tangential vector of a face of the reference cell.
569  * The vectors are arranged such that the
570  * cross product between the two vectors returns the unit normal vector.
571  *
572  * @pre @f$i@f$ must be between zero and `dim-1`.
573  */
574  template <int dim>
576  unit_tangential_vectors(const unsigned int face_no,
577  const unsigned int i) const;
578 
582  template <int dim>
584  unit_normal_vectors(const unsigned int face_no) const;
585 
591  unsigned int
592  n_face_orientations(const unsigned int face_no) const;
593 
610  template <typename T, std::size_t N>
611  DEAL_II_DEPRECATED_EARLY unsigned char
612  compute_orientation(const std::array<T, N> &vertices_0,
613  const std::array<T, N> &vertices_1) const;
614 
640  template <typename T>
641  unsigned char
643  const ArrayView<const T> &vertices_1) const;
644 
645 
659  template <typename T, std::size_t N>
660  DEAL_II_DEPRECATED_EARLY std::array<T, N>
661  permute_according_orientation(const std::array<T, N> &vertices,
662  const unsigned int orientation) const;
663 
675  template <typename T>
676  boost::container::small_vector<T, 8>
678  const unsigned char orientation) const;
679 
684  faces_for_given_vertex(const unsigned int vertex_index) const;
685 
698  unsigned int
699  exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const;
700 
704  unsigned int
705  exodusii_face_to_deal_face(const unsigned int face_n) const;
706 
710  unsigned int
711  unv_vertex_to_deal_vertex(const unsigned int vertex_n) const;
712 
716  unsigned int
717  vtk_linear_type() const;
718 
723  unsigned int
724  vtk_quadratic_type() const;
725 
730  unsigned int
731  vtk_lagrange_type() const;
732 
745  template <int dim>
746  unsigned int
748  const std::array<unsigned, dim> &node_indices,
749  const std::array<unsigned, dim> &nodes_per_direction,
750  const bool legacy_format) const;
751 
755  unsigned int
756  vtk_vertex_to_deal_vertex(const unsigned int vertex_index) const;
757 
761  unsigned int
762  gmsh_element_type() const;
763 
777  std::string
778  to_string() const;
779 
783  constexpr operator std::uint8_t() const;
784 
788  constexpr bool
789  operator==(const ReferenceCell &type) const;
790 
794  constexpr bool
795  operator!=(const ReferenceCell &type) const;
796 
802  template <class Archive>
803  void
804  serialize(Archive &archive, const unsigned int /*version*/);
805 
809  static constexpr std::size_t
811 
816 private:
820  std::uint8_t kind;
821 
828  constexpr ReferenceCell(const std::uint8_t kind);
829 
834  {{{1, 0}}, {{0, 1}}}};
835 
836 
841  {{{0, 2, 1}},
842  {{0, 1, 2}},
843  {{2, 1, 0}},
844  {{1, 2, 0}},
845  {{1, 0, 2}},
846  {{2, 0, 1}}}};
847 
851  static constexpr ndarray<unsigned int, 8, 4>
853  {{0, 1, 2, 3}},
854  {{2, 3, 0, 1}},
855  {{2, 0, 3, 1}},
856  {{3, 1, 2, 0}},
857  {{3, 2, 1, 0}},
858  {{1, 0, 3, 2}},
859  {{1, 3, 0, 2}}}};
860 
865  friend constexpr ReferenceCell
867 
868  friend std::ostream &
869  operator<<(std::ostream &out, const ReferenceCell &reference_cell);
870 
871  friend std::istream &
872  operator>>(std::istream &in, ReferenceCell &reference_cell);
873 };
874 
875 
883 std::ostream &
884 operator<<(std::ostream &out, const ReferenceCell &reference_cell);
885 
893 std::istream &
894 operator>>(std::istream &in, ReferenceCell &reference_cell);
895 
896 
897 inline constexpr ReferenceCell::ReferenceCell(const std::uint8_t kind)
898  : kind(kind)
899 {}
900 
901 
902 
903 inline constexpr ReferenceCell::operator std::uint8_t() const
904 {
905  return kind;
906 }
907 
908 
909 
910 inline constexpr bool
912 {
913  return kind == type.kind;
914 }
915 
916 
917 
918 inline constexpr bool
920 {
921  return kind != type.kind;
922 }
923 
924 
925 
926 namespace internal
927 {
928  inline constexpr ReferenceCell
929  make_reference_cell_from_int(const std::uint8_t kind)
930  {
931 #ifndef DEAL_II_CXX14_CONSTEXPR_BUG
932  // Make sure these are the only indices from which objects can be
933  // created.
934  Assert((kind == std::numeric_limits<std::uint8_t>::max()) || (kind < 8),
935  ExcInternalError());
936 #endif
937 
938  // Call the private constructor, which we can from here because this
939  // function is a 'friend'.
940  return {kind};
941  }
942 } // namespace internal
943 
944 
945 
953 namespace ReferenceCells
954 {
955  constexpr const ReferenceCell Vertex =
957  constexpr const ReferenceCell Line =
959  constexpr const ReferenceCell Triangle =
961  constexpr const ReferenceCell Quadrilateral =
963  constexpr const ReferenceCell Tetrahedron =
965  constexpr const ReferenceCell Pyramid =
967  constexpr const ReferenceCell Wedge =
969  constexpr const ReferenceCell Hexahedron =
971  constexpr const ReferenceCell Invalid =
974 
980  template <int dim>
981  constexpr const ReferenceCell &
982  get_simplex();
983 
989  template <int dim>
990  constexpr const ReferenceCell &
991  get_hypercube();
992 } // namespace ReferenceCells
993 
994 
995 
998 {}
999 
1000 
1001 
1002 template <class Archive>
1003 inline void
1004 ReferenceCell::serialize(Archive &archive, const unsigned int /*version*/)
1005 {
1006  archive &kind;
1007 }
1008 
1009 
1010 
1011 inline constexpr std::size_t
1013 {
1014  return sizeof(ReferenceCells::Invalid);
1015 }
1016 
1017 
1018 
1020 ReferenceCell::faces_for_given_vertex(const unsigned int vertex) const
1021 {
1023  switch (this->kind)
1024  {
1025  case ReferenceCells::Line:
1026  return {&GeometryInfo<2>::vertex_to_face[vertex][0], 1};
1028  return {&GeometryInfo<2>::vertex_to_face[vertex][0], 2};
1030  {
1031  static constexpr ndarray<unsigned int, 3, 2> table = {
1032  {{{0, 2}}, {{0, 1}}, {{1, 2}}}};
1033  return table[vertex];
1034  }
1036  {
1037  static constexpr ndarray<unsigned int, 4, 3> table = {
1038  {{{0, 1, 2}}, {{0, 1, 3}}, {{0, 2, 3}}, {{1, 2, 3}}}};
1039 
1040  return table[vertex];
1041  }
1043  {
1044  static constexpr unsigned int X = numbers::invalid_unsigned_int;
1045  static constexpr ndarray<unsigned int, 5, 4> table = {
1046  {{{0, 1, 3, X}},
1047  {{0, 2, 3, X}},
1048  {{0, 1, 4, X}},
1049  {{0, 2, 4, X}},
1050  {{1, 2, 3, 4}}}};
1051 
1052  return {&table[vertex][0], vertex == 4 ? 4u : 3u};
1053  }
1054  case ReferenceCells::Wedge:
1055  {
1057  static constexpr ndarray<unsigned int, 6, 3> table = {{{{0, 2, 4}},
1058  {{0, 2, 3}},
1059  {{0, 3, 4}},
1060  {{1, 2, 4}},
1061  {{1, 2, 3}},
1062  {{1, 3, 4}}}};
1063 
1064  return table[vertex];
1065  }
1067  return {&GeometryInfo<3>::vertex_to_face[vertex][0], 3};
1068  default:
1069  Assert(false, ExcNotImplemented());
1070  }
1071 
1072  return {};
1073 }
1074 
1075 
1076 
1077 inline bool
1079 {
1080  return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
1081  *this == ReferenceCells::Quadrilateral ||
1082  *this == ReferenceCells::Hexahedron);
1083 }
1084 
1085 
1086 
1087 inline bool
1089 {
1090  return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
1091  *this == ReferenceCells::Triangle ||
1092  *this == ReferenceCells::Tetrahedron);
1093 }
1094 
1095 
1096 
1097 inline unsigned int
1099 {
1100  switch (this->kind)
1101  {
1103  return 0;
1104  case ReferenceCells::Line:
1105  return 1;
1108  return 2;
1111  case ReferenceCells::Wedge:
1113  return 3;
1114  default:
1115  Assert(false, ExcNotImplemented());
1116  }
1117 
1119 }
1120 
1121 
1122 
1123 template <int dim>
1126 {
1127  return Quadrature<dim>(std::vector<Point<dim>>({barycenter<dim>()}),
1128  std::vector<double>({volume()}));
1129 }
1130 
1131 
1132 
1133 inline unsigned int
1135 {
1136  switch (this->kind)
1137  {
1139  return 1;
1140  case ReferenceCells::Line:
1141  return 2;
1143  return 3;
1145  return 4;
1147  return 4;
1149  return 5;
1150  case ReferenceCells::Wedge:
1151  return 6;
1153  return 8;
1154  default:
1155  Assert(false, ExcNotImplemented());
1156  }
1157 
1159 }
1160 
1161 
1162 
1163 inline unsigned int
1165 {
1166  switch (this->kind)
1167  {
1169  return 0;
1170  case ReferenceCells::Line:
1171  return 1;
1173  return 3;
1175  return 4;
1177  return 6;
1179  return 7;
1180  case ReferenceCells::Wedge:
1181  return 9;
1183  return 12;
1184  default:
1185  Assert(false, ExcNotImplemented());
1186  }
1187 
1189 }
1190 
1191 
1192 
1193 template <int dim>
1194 Point<dim>
1195 ReferenceCell::vertex(const unsigned int v) const
1196 {
1199 
1200  switch (dim)
1201  {
1202  case 0:
1203  {
1204  if (*this == ReferenceCells::Vertex)
1205  return Point<dim>(0);
1206  break;
1207  }
1208  case 1:
1209  {
1210  static const Point<dim> vertices[2] = {
1211  Point<dim>(), // the origin
1212  Point<dim>::unit_vector(0) // unit point along x-axis
1213  };
1214  if (*this == ReferenceCells::Line)
1215  return vertices[v];
1216  break;
1217  }
1218  case 2:
1219  {
1220  switch (this->kind)
1221  {
1223  {
1224  static const Point<dim> vertices[3] = {
1225  Point<dim>(), // the origin
1226  Point<dim>::unit_vector(0), // unit point along x-axis
1227  Point<dim>::unit_vector(1) // unit point along y-axis
1228  };
1229  return vertices[v];
1230  }
1232  {
1233  static const Point<dim> vertices[4] = {
1234  // First the two points on the x-axis
1235  Point<dim>(),
1237  // Then these two points shifted in the y-direction
1240  return vertices[v];
1241  }
1242  }
1243  break;
1244  }
1245  case 3:
1246  {
1247  switch (this->kind)
1248  {
1250  {
1251  static const Point<dim> vertices[4] = {
1252  Point<dim>(), // the origin
1253  Point<dim>::unit_vector(0), // unit point along x-axis
1254  Point<dim>::unit_vector(1), // unit point along y-axis
1255  Point<dim>::unit_vector(2) // unit point along z-axis
1256  };
1257  return vertices[v];
1258  }
1260  {
1261  static const Point<dim> vertices[5] = {
1262  Point<dim>{-1.0, -1.0, 0.0},
1263  Point<dim>{+1.0, -1.0, 0.0},
1264  Point<dim>{-1.0, +1.0, 0.0},
1265  Point<dim>{+1.0, +1.0, 0.0},
1266  Point<dim>{+0.0, +0.0, 1.0}};
1267  return vertices[v];
1268  }
1269  case ReferenceCells::Wedge:
1270  {
1271  static const Point<dim> vertices[6] = {
1272  // First the three points on the triangular base of the
1273  // wedge:
1274  Point<dim>(),
1277  // And now everything shifted in the z-direction again
1281  return vertices[v];
1282  }
1284  {
1285  static const Point<dim> vertices[8] = {
1286  // First the two points on the x-axis
1287  Point<dim>(),
1289  // Then these two points shifted in the y-direction
1292  // And now all four points shifted in the z-direction
1299  return vertices[v];
1300  }
1301  }
1302  break;
1303  }
1304  default:
1305  Assert(false, ExcNotImplemented());
1306  }
1307 
1308  Assert(false, ExcNotImplemented());
1309  return Point<dim>();
1310 }
1311 
1312 
1313 inline unsigned int
1315 {
1316  switch (this->kind)
1317  {
1319  return 0;
1320  case ReferenceCells::Line:
1321  return 2;
1323  return 3;
1325  return 4;
1327  return 4;
1329  return 5;
1330  case ReferenceCells::Wedge:
1331  return 5;
1333  return 6;
1334  default:
1335  Assert(false, ExcNotImplemented());
1336  }
1337 
1339 }
1340 
1341 
1342 
1345 {
1346  return {0U, n_faces()};
1347 }
1348 
1349 
1350 
1351 inline unsigned int
1353 {
1354  switch (this->kind)
1355  {
1357  return 0;
1358  case ReferenceCells::Line:
1359  return 2;
1361  return 4;
1363  return 4;
1365  return 8;
1367  // We haven't yet decided how to refine pyramids. Update this when we
1368  // have
1369  Assert(false, ExcNotImplemented());
1371  case ReferenceCells::Wedge:
1372  return 8;
1374  return 8;
1375  default:
1376  Assert(false, ExcNotImplemented());
1377  }
1378 
1380 }
1381 
1382 
1383 
1386 {
1387  return {0U, n_isotropic_children()};
1388 }
1389 
1390 
1391 
1394 {
1395  return {0U, n_vertices()};
1396 }
1397 
1398 
1399 
1402 {
1403  return {0U, n_lines()};
1404 }
1405 
1406 
1407 
1408 inline ReferenceCell
1409 ReferenceCell::face_reference_cell(const unsigned int face_no) const
1410 {
1411  AssertIndexRange(face_no, n_faces());
1412 
1413  switch (this->kind)
1414  {
1416  return ReferenceCells::Invalid;
1417  case ReferenceCells::Line:
1418  return ReferenceCells::Vertex;
1421  return ReferenceCells::Line;
1423  return ReferenceCells::Triangle;
1425  if (face_no == 0)
1427  else
1428  return ReferenceCells::Triangle;
1429  case ReferenceCells::Wedge:
1430  if (face_no > 1)
1432  else
1433  return ReferenceCells::Triangle;
1436  default:
1437  Assert(false, ExcNotImplemented());
1438  }
1439 
1440  return ReferenceCells::Invalid;
1441 }
1442 
1443 
1444 
1445 inline constexpr unsigned char
1447 {
1448  // Our convention is that 'orientation' has a default value of true and
1449  // occupies the least-significant bit while rotate and flip have default
1450  // values of 'false' and occupy the second and third bits.
1451  return 0b001;
1452 }
1453 
1454 
1455 
1456 inline constexpr unsigned char
1458 {
1459  // For a reversed line 'orientation' is false and neither flip nor rotate are
1460  // defined.
1461  return 0b000;
1462 }
1463 
1464 
1465 
1466 inline unsigned int
1468  const unsigned int face,
1469  const unsigned int subface,
1470  const unsigned char combined_face_orientation) const
1471 {
1472  AssertIndexRange(face, n_faces());
1474 
1475  switch (this->kind)
1476  {
1478  case ReferenceCells::Line:
1479  {
1480  Assert(false, ExcNotImplemented());
1481  break;
1482  }
1484  {
1485  static constexpr ndarray<unsigned int, 3, 2> subcells = {
1486  {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
1487 
1488  return subcells[face][subface];
1489  }
1491  {
1492  const bool face_orientation =
1493  Utilities::get_bit(combined_face_orientation, 0);
1494  const bool face_flip =
1495  Utilities::get_bit(combined_face_orientation, 2);
1496  const bool face_rotation =
1497  Utilities::get_bit(combined_face_orientation, 1);
1498 
1501  face,
1502  subface,
1503  face_orientation,
1504  face_flip,
1505  face_rotation);
1506  }
1509  case ReferenceCells::Wedge:
1510  {
1511  Assert(false, ExcNotImplemented());
1512  break;
1513  }
1515  {
1516  const bool face_orientation =
1517  Utilities::get_bit(combined_face_orientation, 0);
1518  const bool face_flip =
1519  Utilities::get_bit(combined_face_orientation, 2);
1520  const bool face_rotation =
1521  Utilities::get_bit(combined_face_orientation, 1);
1522 
1525  face,
1526  subface,
1527  face_orientation,
1528  face_flip,
1529  face_rotation);
1530  }
1531  default:
1532  Assert(false, ExcNotImplemented());
1533  }
1534 
1536 }
1537 
1538 
1539 
1540 inline std::array<unsigned int, 2>
1542  const unsigned int vertex) const
1543 {
1545  // Work around a GCC warning at higher optimization levels by making all of
1546  // these tables the same size
1547  constexpr unsigned int X = numbers::invalid_unsigned_int;
1548 
1549  switch (this->kind)
1550  {
1552  case ReferenceCells::Line:
1553  Assert(false, ExcNotImplemented());
1554  break;
1556  {
1557  static constexpr ndarray<unsigned int, 6, 2> table = {
1558  {{{0, 0}}, {{0, 1}}, {{1, 1}}, {{X, X}}, {{X, X}}, {{X, X}}}};
1559 
1560  return table[vertex];
1561  }
1563  {
1565  vertex);
1566  }
1568  {
1569  static constexpr ndarray<unsigned int, 6, 2> table = {
1570  {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 2}}, {{X, X}}, {{X, X}}}};
1571 
1572  return table[vertex];
1573  }
1575  {
1576  static constexpr ndarray<unsigned int, 6, 2> table = {
1577  {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{0, 3}}, {{1, 2}}, {{X, X}}}};
1578 
1579  return table[vertex];
1580  }
1581  case ReferenceCells::Wedge:
1582  {
1583  static constexpr ndarray<unsigned int, 6, 2> table = {
1584  {{{0, 1}}, {{0, 0}}, {{0, 2}}, {{1, 0}}, {{1, 1}}, {{1, 2}}}};
1585 
1586  return table[vertex];
1587  }
1589  {
1591  vertex);
1592  }
1593  default:
1594  Assert(false, ExcNotImplemented());
1595  }
1596 
1597  return {};
1598 }
1599 
1600 
1601 
1602 inline std::array<unsigned int, 2>
1604  const unsigned int line) const
1605 {
1606  AssertIndexRange(line, n_lines());
1607 
1608  switch (this->kind)
1609  {
1611  case ReferenceCells::Line:
1614  {
1615  Assert(false, ExcNotImplemented());
1616  break;
1617  }
1619  {
1620  static const std::array<unsigned int, 2> table[6] = {
1621  {{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}};
1622 
1623  return table[line];
1624  }
1626  {
1627  static const std::array<unsigned int, 2> table[8] = {{{0, 0}},
1628  {{0, 1}},
1629  {{0, 2}},
1630  {{0, 3}},
1631  {{1, 2}},
1632  {{2, 1}},
1633  {{1, 1}},
1634  {{2, 2}}};
1635 
1636  return table[line];
1637  }
1638  case ReferenceCells::Wedge:
1639  {
1640  static const std::array<unsigned int, 2> table[9] = {{{0, 0}},
1641  {{0, 2}},
1642  {{0, 1}},
1643  {{1, 0}},
1644  {{1, 1}},
1645  {{1, 2}},
1646  {{2, 0}},
1647  {{2, 1}},
1648  {{3, 1}}};
1649 
1650  return table[line];
1651  }
1653  {
1655  }
1656  default:
1657  Assert(false, ExcNotImplemented());
1658  }
1659 
1660  return {};
1661 }
1662 
1663 
1664 inline unsigned int
1666  const unsigned int vertex) const
1667 {
1669  AssertIndexRange(line, n_lines());
1670 
1671  switch (this->kind)
1672  {
1674  case ReferenceCells::Line:
1675  return vertex;
1677  {
1678  static constexpr ndarray<unsigned int, 3, 2> table = {
1679  {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
1680  return table[line][vertex];
1681  }
1683  {
1684  static constexpr ndarray<unsigned int, 4, 2> table = {
1685  {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
1686  return table[line][vertex];
1687  }
1689  {
1690  static constexpr ndarray<unsigned int, 6, 2> table = {
1691  {{{0, 1}}, {{1, 2}}, {{2, 0}}, {{0, 3}}, {{1, 3}}, {{2, 3}}}};
1692  return table[line][vertex];
1693  }
1695  {
1696  static constexpr ndarray<unsigned int, 8, 2> table = {{{{0, 2}},
1697  {{1, 3}},
1698  {{0, 1}},
1699  {{2, 3}},
1700  {{4, 0}},
1701  {{1, 4}},
1702  {{2, 4}},
1703  {{4, 3}}}};
1704  return table[line][vertex];
1705  }
1706  case ReferenceCells::Wedge:
1707  {
1708  static constexpr ndarray<unsigned int, 9, 2> table = {{{{1, 0}},
1709  {{2, 1}},
1710  {{0, 2}},
1711  {{3, 4}},
1712  {{4, 5}},
1713  {{5, 3}},
1714  {{0, 3}},
1715  {{1, 4}},
1716  {{2, 5}}}};
1717  return table[line][vertex];
1718  }
1720  {
1721  // first four lines comprise the bottom face, next four are the top,
1722  // and the last four are 'bottom to top'
1723  static constexpr ndarray<unsigned int, 12, 2> table = {{{{0, 2}},
1724  {{1, 3}},
1725  {{0, 1}},
1726  {{2, 3}},
1727  {{4, 6}},
1728  {{5, 7}},
1729  {{4, 5}},
1730  {{6, 7}},
1731  {{0, 4}},
1732  {{1, 5}},
1733  {{2, 6}},
1734  {{3, 7}}}};
1735  return table[line][vertex];
1736  }
1737 
1738  default:
1739  Assert(false, ExcNotImplemented());
1740  }
1741 
1743 }
1744 
1745 
1746 inline unsigned int
1747 ReferenceCell::face_to_cell_lines(const unsigned int face,
1748  const unsigned int line,
1749  const unsigned char face_orientation) const
1750 {
1751  AssertIndexRange(face, n_faces());
1753 
1754  static constexpr unsigned int X = numbers::invalid_unsigned_int;
1755 
1756  switch (this->kind)
1757  {
1759  {
1760  Assert(false, ExcNotImplemented());
1761  break;
1762  }
1763  case ReferenceCells::Line:
1764  {
1766  face,
1767  line,
1768  Utilities::get_bit(face_orientation, 0),
1769  Utilities::get_bit(face_orientation, 2),
1770  Utilities::get_bit(face_orientation, 1));
1771  }
1773  {
1774  return face;
1775  }
1777  {
1779  face,
1780  line,
1781  Utilities::get_bit(face_orientation, 0),
1782  Utilities::get_bit(face_orientation, 2),
1783  Utilities::get_bit(face_orientation, 1));
1784  }
1786  {
1787  static constexpr ndarray<unsigned int, 4, 3> table = {
1788  {{{0, 1, 2}}, {{0, 3, 4}}, {{2, 5, 3}}, {{1, 4, 5}}}};
1789 
1790  return table[face][standard_to_real_face_line(
1791  line, face, face_orientation)];
1792  }
1794  {
1795  static constexpr ndarray<unsigned int, 5, 4> table = {
1796  {{{0, 1, 2, 3}},
1797  {{0, 6, 4, X}},
1798  {{1, 5, 7, X}},
1799  {{2, 4, 5, X}},
1800  {{3, 7, 6, 2}}}};
1801 
1802  return table[face][standard_to_real_face_line(
1803  line, face, face_orientation)];
1804  }
1805  case ReferenceCells::Wedge:
1806  {
1807  static constexpr ndarray<unsigned int, 5, 4> table = {
1808  {{{0, 2, 1, X}},
1809  {{3, 4, 5, X}},
1810  {{6, 7, 0, 3}},
1811  {{7, 8, 1, 4}},
1812  {{8, 6, 5, 2}}}};
1813 
1814  return table[face][standard_to_real_face_line(
1815  line, face, face_orientation)];
1816  }
1818  {
1820  face,
1821  line,
1822  Utilities::get_bit(face_orientation, 0),
1823  Utilities::get_bit(face_orientation, 2),
1824  Utilities::get_bit(face_orientation, 1));
1825  }
1826  default:
1827  Assert(false, ExcNotImplemented());
1828  }
1829 
1831 }
1832 
1833 
1834 
1835 inline unsigned int
1837  const unsigned int vertex,
1838  const unsigned char face_orientation) const
1839 {
1840  AssertIndexRange(face, n_faces());
1842 
1843  switch (this->kind)
1844  {
1846  {
1847  Assert(false, ExcNotImplemented());
1848  break;
1849  }
1850  case ReferenceCells::Line:
1851  {
1853  face,
1854  vertex,
1855  Utilities::get_bit(face_orientation, 0),
1856  Utilities::get_bit(face_orientation, 2),
1857  Utilities::get_bit(face_orientation, 1));
1858  }
1860  {
1861  static constexpr ndarray<unsigned int, 3, 2> table = {
1862  {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
1863 
1864  return table[face][face_orientation != 0u ? vertex : (1 - vertex)];
1865  }
1867  {
1869  face,
1870  vertex,
1871  Utilities::get_bit(face_orientation, 0),
1872  Utilities::get_bit(face_orientation, 2),
1873  Utilities::get_bit(face_orientation, 1));
1874  }
1876  {
1877  static constexpr ndarray<unsigned int, 4, 3> table = {
1878  {{{0, 1, 2}}, {{1, 0, 3}}, {{0, 2, 3}}, {{2, 1, 3}}}};
1879 
1880  return table[face][standard_to_real_face_vertex(
1881  vertex, face, face_orientation)];
1882  }
1884  {
1885  constexpr auto X = numbers::invalid_unsigned_int;
1886  static constexpr ndarray<unsigned int, 5, 4> table = {
1887  {{{0, 1, 2, 3}},
1888  {{0, 2, 4, X}},
1889  {{3, 1, 4, X}},
1890  {{1, 0, 4, X}},
1891  {{2, 3, 4, X}}}};
1892 
1893  return table[face][standard_to_real_face_vertex(
1894  vertex, face, face_orientation)];
1895  }
1896  case ReferenceCells::Wedge:
1897  {
1898  constexpr auto X = numbers::invalid_unsigned_int;
1899  static constexpr ndarray<unsigned int, 6, 4> table = {
1900  {{{1, 0, 2, X}},
1901  {{3, 4, 5, X}},
1902  {{0, 1, 3, 4}},
1903  {{1, 2, 4, 5}},
1904  {{2, 0, 5, 3}}}};
1905 
1906  return table[face][standard_to_real_face_vertex(
1907  vertex, face, face_orientation)];
1908  }
1910  {
1912  face,
1913  vertex,
1914  Utilities::get_bit(face_orientation, 0),
1915  Utilities::get_bit(face_orientation, 2),
1916  Utilities::get_bit(face_orientation, 1));
1917  }
1918  default:
1919  Assert(false, ExcNotImplemented());
1920  }
1921 
1923 }
1924 
1925 
1926 
1927 inline unsigned int
1929  const unsigned int vertex,
1930  const unsigned int face,
1931  const unsigned char face_orientation) const
1932 {
1933  AssertIndexRange(face, n_faces());
1935 
1936  switch (this->kind)
1937  {
1939  case ReferenceCells::Line:
1940  Assert(false, ExcNotImplemented());
1941  break;
1944  return line_vertex_permutations[face_orientation][vertex];
1946  return triangle_vertex_permutations[face_orientation][vertex];
1948  // face 0 is a quadrilateral
1949  if (face == 0)
1950  return quadrilateral_vertex_permutations[face_orientation][vertex];
1951  else
1952  return triangle_vertex_permutations[face_orientation][vertex];
1953  case ReferenceCells::Wedge:
1954  // faces 0 and 1 are triangles
1955  if (face > 1)
1956  return quadrilateral_vertex_permutations[face_orientation][vertex];
1957  else
1958  return triangle_vertex_permutations[face_orientation][vertex];
1960  return quadrilateral_vertex_permutations[face_orientation][vertex];
1961  default:
1962  Assert(false, ExcNotImplemented());
1963  }
1964 
1965  Assert(false, ExcNotImplemented());
1967 }
1968 
1969 
1970 
1971 inline unsigned int
1973  const unsigned int line,
1974  const unsigned int face,
1975  const unsigned char face_orientation) const
1976 {
1977  AssertIndexRange(face, n_faces());
1979 
1980  static constexpr ndarray<unsigned int, 6, 3> triangle_table = {{{{2, 1, 0}},
1981  {{0, 1, 2}},
1982  {{1, 0, 2}},
1983  {{1, 2, 0}},
1984  {{0, 2, 1}},
1985  {{2, 0, 1}}}};
1986 
1987 
1988  switch (this->kind)
1989  {
1991  case ReferenceCells::Line:
1994  Assert(false, ExcNotImplemented());
1995  break;
1997  return triangle_table[face_orientation][line];
1999  if (face == 0) // The quadrilateral face
2000  {
2002  line,
2003  Utilities::get_bit(face_orientation, 0),
2004  Utilities::get_bit(face_orientation, 2),
2005  Utilities::get_bit(face_orientation, 1));
2006  }
2007  else // One of the triangular faces
2008  {
2009  return triangle_table[face_orientation][line];
2010  }
2011  case ReferenceCells::Wedge:
2012  if (face > 1) // One of the quadrilateral faces
2013  {
2015  line,
2016  Utilities::get_bit(face_orientation, 0),
2017  Utilities::get_bit(face_orientation, 2),
2018  Utilities::get_bit(face_orientation, 1));
2019  }
2020  else // One of the triangular faces
2021  return triangle_table[face_orientation][line];
2023  {
2024  static constexpr ndarray<unsigned int, 8, 4> table = {
2025  {{{2, 3, 0, 1}},
2026  {{0, 1, 2, 3}},
2027  {{0, 1, 3, 2}},
2028  {{3, 2, 0, 1}},
2029  {{3, 2, 1, 0}},
2030  {{1, 0, 3, 2}},
2031  {{1, 0, 2, 3}},
2032  {{2, 3, 1, 0}}}};
2033  return table[face_orientation][line];
2034  }
2035  default:
2036  Assert(false, ExcNotImplemented());
2037  }
2038 
2040 }
2041 
2042 
2043 
2044 namespace ReferenceCells
2045 {
2046  template <int dim>
2047  inline constexpr const ReferenceCell &
2049  {
2050  switch (dim)
2051  {
2052  case 0:
2053  return ReferenceCells::Vertex;
2054  case 1:
2055  return ReferenceCells::Line;
2056  case 2:
2057  return ReferenceCells::Triangle;
2058  case 3:
2060  default:
2061  Assert(false, ExcNotImplemented());
2062  }
2063  return ReferenceCells::Invalid;
2064  }
2065 
2066 
2067 
2068  template <int dim>
2069  inline constexpr const ReferenceCell &
2071  {
2072  switch (dim)
2073  {
2074  case 0:
2075  return ReferenceCells::Vertex;
2076  case 1:
2077  return ReferenceCells::Line;
2078  case 2:
2080  case 3:
2082  default:
2083  Assert(false, ExcNotImplemented());
2084  }
2085  return ReferenceCells::Invalid;
2086  }
2087 } // namespace ReferenceCells
2088 
2089 
2090 inline ReferenceCell
2091 ReferenceCell::n_vertices_to_type(const int dim, const unsigned int n_vertices)
2092 {
2093  AssertIndexRange(dim, 4);
2095 
2096  const auto X = ReferenceCells::Invalid;
2097  static const ndarray<ReferenceCell, 4, 9> table = {
2098  {// dim 0
2099  {{X, ReferenceCells::Vertex, X, X, X, X, X, X, X}},
2100  // dim 1
2101  {{X, X, ReferenceCells::Line, X, X, X, X, X, X}},
2102  // dim 2
2103  {{X,
2104  X,
2105  X,
2108  X,
2109  X,
2110  X,
2111  X}},
2112  // dim 3
2113  {{X,
2114  X,
2115  X,
2116  X,
2120  X,
2122  Assert(table[dim][n_vertices] != ReferenceCells::Invalid,
2123  ExcMessage("The combination of dim = " + std::to_string(dim) +
2124  " and n_vertices = " + std::to_string(n_vertices) +
2125  " does not correspond to a known reference cell type."));
2126  return table[dim][n_vertices];
2127 }
2128 
2129 
2130 
2131 template <int dim>
2132 inline double
2134  const unsigned int i) const
2135 {
2138  switch (this->kind)
2139  {
2141  case ReferenceCells::Line:
2145  // see also BarycentricPolynomials<2>::compute_value
2147  {
2148  switch (i)
2149  {
2150  case 0:
2151  return 1.0 - xi[std::min(0, dim - 1)] -
2152  xi[std::min(1, dim - 1)];
2153  case 1:
2154  return xi[std::min(0, dim - 1)];
2155  case 2:
2156  return xi[std::min(1, dim - 1)];
2157  default:
2158  Assert(false, ExcInternalError());
2159  }
2160  }
2161  // see also BarycentricPolynomials<3>::compute_value
2163  {
2164  switch (i)
2165  {
2166  case 0:
2167  return 1.0 - xi[std::min(0, dim - 1)] -
2168  xi[std::min(1, dim - 1)] - xi[std::min(2, dim - 1)];
2169  case 1:
2170  return xi[std::min(0, dim - 1)];
2171  case 2:
2172  return xi[std::min(1, dim - 1)];
2173  case 3:
2174  return xi[std::min(2, dim - 1)];
2175  default:
2176  Assert(false, ExcInternalError());
2177  }
2178  }
2179  // see also ScalarLagrangePolynomialPyramid::compute_value()
2181  {
2182  const double Q14 = 0.25;
2183 
2184  const double r = xi[std::min(0, dim - 1)];
2185  const double s = xi[std::min(1, dim - 1)];
2186  const double t = xi[std::min(2, dim - 1)];
2187 
2188  const double ratio =
2189  (std::fabs(t - 1.0) > 1.0e-14 ? (r * s * t) / (1.0 - t) : 0.0);
2190 
2191  if (i == 0)
2192  return Q14 * ((1.0 - r) * (1.0 - s) - t + ratio);
2193  if (i == 1)
2194  return Q14 * ((1.0 + r) * (1.0 - s) - t - ratio);
2195  if (i == 2)
2196  return Q14 * ((1.0 - r) * (1.0 + s) - t - ratio);
2197  if (i == 3)
2198  return Q14 * ((1.0 + r) * (1.0 + s) - t + ratio);
2199  else
2200  return t;
2201  }
2202  // see also ScalarLagrangePolynomialWedge::compute_value()
2203  case ReferenceCells::Wedge:
2205  .d_linear_shape_function<2>(Point<2>(xi[std::min(0, dim - 1)],
2206  xi[std::min(1, dim - 1)]),
2207  i % 3) *
2209  .d_linear_shape_function<1>(Point<1>(xi[std::min(2, dim - 1)]),
2210  i / 3);
2211  default:
2212  Assert(false, ExcNotImplemented());
2213  }
2214 
2215  return 0.0;
2216 }
2217 
2218 
2219 
2220 template <int dim>
2221 inline Tensor<1, dim>
2223  const unsigned int i) const
2224 {
2226  switch (this->kind)
2227  {
2229  case ReferenceCells::Line:
2233  // see also BarycentricPolynomials<2>::compute_grad()
2235  switch (i)
2236  {
2237  case 0:
2238  return Point<dim>(-1.0, -1.0);
2239  case 1:
2240  return Point<dim>(+1.0, +0.0);
2241  case 2:
2242  return Point<dim>(+0.0, +1.0);
2243  default:
2244  Assert(false, ExcInternalError());
2245  }
2246  default:
2247  Assert(false, ExcNotImplemented());
2248  }
2249 
2250  return Point<dim>(+0.0, +0.0, +0.0);
2251 }
2252 
2253 
2254 
2255 inline double
2257 {
2258  switch (this->kind)
2259  {
2261  return 0;
2262  case ReferenceCells::Line:
2263  return 1;
2265  return 1. / 2.;
2267  return 1;
2269  return 1. / 6.;
2271  return 4. / 3.;
2272  case ReferenceCells::Wedge:
2273  return 1. / 2.;
2275  return 1;
2276  default:
2277  Assert(false, ExcNotImplemented());
2278  }
2279 
2280  return 0.0;
2281 }
2282 
2283 
2284 
2285 template <int dim>
2286 inline Point<dim>
2288 {
2290 
2291  switch (this->kind)
2292  {
2294  return Point<dim>();
2295  case ReferenceCells::Line:
2296  return Point<dim>(1. / 2.);
2298  return Point<dim>(1. / 3., 1. / 3.);
2300  return Point<dim>(1. / 2., 1. / 2.);
2302  return Point<dim>(1. / 4., 1. / 4., 1. / 4.);
2304  return Point<dim>(0, 0, 1. / 4.);
2305  case ReferenceCells::Wedge:
2306  return Point<dim>(1. / 3, 1. / 3, 1. / 2.);
2308  return Point<dim>(1. / 2., 1. / 2., 1. / 2.);
2309  default:
2310  Assert(false, ExcNotImplemented());
2311  }
2312 
2313  return Point<dim>();
2314 }
2315 
2316 
2317 
2318 template <int dim>
2319 inline bool
2320 ReferenceCell::contains_point(const Point<dim> &p, const double tolerance) const
2321 {
2323 
2324  // Introduce abbreviations to silence compiler warnings about invalid
2325  // array accesses (that can't happen, of course, but the compiler
2326  // doesn't know that).
2327  constexpr unsigned int x_coordinate = 0;
2328  constexpr unsigned int y_coordinate = (dim >= 2 ? 1 : 0);
2329  constexpr unsigned int z_coordinate = (dim >= 3 ? 2 : 0);
2330 
2331  switch (this->kind)
2332  {
2334  {
2335  // Vertices are special cases in that they do not actually
2336  // have coordinates. Error out if this function is called
2337  // with a vertex:
2338  Assert(false,
2339  ExcMessage("Vertices are zero-dimensional objects and "
2340  "as a consequence have no coordinates. You "
2341  "cannot meaningfully ask whether a point is "
2342  "inside a vertex (within a certain tolerance) "
2343  "without coordinate values."));
2344  return false;
2345  }
2346  case ReferenceCells::Line:
2349  {
2350  for (unsigned int d = 0; d < dim; ++d)
2351  if ((p[d] < -tolerance) || (p[d] > 1 + tolerance))
2352  return false;
2353  return true;
2354  }
2357  {
2358  // First make sure that we are in the first quadrant or octant
2359  for (unsigned int d = 0; d < dim; ++d)
2360  if (p[d] < -tolerance)
2361  return false;
2362 
2363  // Now we also need to make sure that we are below the diagonal line
2364  // or plane that delineates the simplex. This diagonal is given by
2365  // sum(p[d])<=1, and a diagonal a distance eps away is given by
2366  // sum(p[d])<=1+eps*sqrt(d). (For example, the point at (1,1) is a
2367  // distance of 1/sqrt(2) away from the diagonal. That is, its
2368  // sum satisfies
2369  // sum(p[d]) = 2 <= 1 + (1/sqrt(2)) * sqrt(2)
2370  // in other words, it satisfies the predicate with eps=1/sqrt(2).)
2371  double sum = 0;
2372  for (unsigned int d = 0; d < dim; ++d)
2373  sum += p[d];
2374  return (sum <= 1 + tolerance * std::sqrt(1. * dim));
2375  }
2377  {
2378  // A pyramid only lives in the upper half-space:
2379  if (p[z_coordinate] < -tolerance)
2380  return false;
2381 
2382  // It also only lives in the space below z=1:
2383  if (p[z_coordinate] > 1 + tolerance)
2384  return false;
2385 
2386  // Within what's left of the space, a pyramid is a cone that tapers
2387  // towards the top. First compute the distance of the point to the
2388  // axis in the max norm (this is the right norm because the vertices
2389  // of the pyramid are at points +/-1, +/-1):
2390  const double distance_from_axis =
2391  std::max(std::fabs(p[x_coordinate]), std::fabs(p[y_coordinate]));
2392 
2393  // We are inside the pyramid if the distance from the axis is less
2394  // than (1-z)
2395  return (distance_from_axis < 1 + tolerance - p[z_coordinate]);
2396  }
2397  case ReferenceCells::Wedge:
2398  {
2399  // The wedge we use is a triangle extruded into the third
2400  // dimension by one unit. So we can use the same logic as for
2401  // triangles above (i.e., for the simplex above, using dim==2)
2402  // and then check the third dimension separately.
2403 
2404  if ((p[x_coordinate] < -tolerance) || (p[y_coordinate] < -tolerance))
2405  return false;
2406 
2407  const double sum = p[x_coordinate] + p[y_coordinate];
2408  if (sum > 1 + tolerance * std::sqrt(2.0))
2409  return false;
2410 
2411  if (p[z_coordinate] < -tolerance)
2412  return false;
2413  if (p[z_coordinate] > 1 + tolerance)
2414  return false;
2415 
2416  return true;
2417  }
2418  default:
2419  Assert(false, ExcNotImplemented());
2420  }
2421 
2422  return false;
2423 }
2424 
2425 
2426 
2427 template <int dim>
2428 inline Tensor<1, dim>
2429 ReferenceCell::unit_tangential_vectors(const unsigned int face_no,
2430  const unsigned int i) const
2431 {
2433  AssertIndexRange(i, dim - 1);
2434 
2435  switch (this->kind)
2436  {
2438  case ReferenceCells::Line:
2442  return GeometryInfo<dim>::unit_tangential_vectors[face_no][i];
2444  {
2445  AssertIndexRange(face_no, 3);
2446  static const std::array<Tensor<1, dim>, 3> table = {
2447  {Point<dim>(1, 0),
2448  Point<dim>(-std::sqrt(0.5), +std::sqrt(0.5)),
2449  Point<dim>(0, -1)}};
2450 
2451  return table[face_no];
2452  }
2454  {
2455  AssertIndexRange(face_no, 4);
2456  static const ndarray<Tensor<1, dim>, 4, 2> table = {
2457  {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
2458  {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
2459  {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
2460  {{Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
2461  +std::pow(1.0 / 3.0, 1.0 / 4.0),
2462  0),
2463  Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
2464  0,
2465  +std::pow(1.0 / 3.0, 1.0 / 4.0))}}}};
2466 
2467  return table[face_no][i];
2468  }
2470  {
2471  AssertIndexRange(face_no, 5);
2472  static const ndarray<Tensor<1, dim>, 5, 2> table = {
2473  {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
2474  {{Point<dim>(+1.0 / sqrt(2.0), 0, +1.0 / sqrt(2.0)),
2475  Point<dim>(0, 1, 0)}},
2476  {{Point<dim>(+1.0 / sqrt(2.0), 0, -1.0 / sqrt(2.0)),
2477  Point<dim>(0, 1, 0)}},
2478  {{Point<dim>(1, 0, 0),
2479  Point<dim>(0, +1.0 / sqrt(2.0), +1.0 / sqrt(2.0))}},
2480  {{Point<dim>(1, 0, 0),
2481  Point<dim>(0, +1.0 / sqrt(2.0), -1.0 / sqrt(2.0))}}}};
2482 
2483  return table[face_no][i];
2484  }
2485  case ReferenceCells::Wedge:
2486  {
2487  AssertIndexRange(face_no, 5);
2488  static const ndarray<Tensor<1, dim>, 5, 2> table = {
2489  {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
2490  {{Point<dim>(1, 0, 0), Point<dim>(0, 1, 0)}},
2491  {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
2492  {{Point<dim>(-1 / std::sqrt(2.0), +1 / std::sqrt(2.0), 0),
2493  Point<dim>(0, 0, 1)}},
2494  {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}}}};
2495 
2496  return table[face_no][i];
2497  }
2498  default:
2499  Assert(false, ExcNotImplemented());
2500  }
2501 
2502 
2503  return {};
2504 }
2505 
2506 
2507 
2508 template <int dim>
2509 inline Tensor<1, dim>
2510 ReferenceCell::unit_normal_vectors(const unsigned int face_no) const
2511 {
2512  AssertDimension(dim, this->get_dimension());
2513 
2514  if (is_hyper_cube())
2515  {
2517  return GeometryInfo<dim>::unit_normal_vector[face_no];
2518  }
2519  else if (dim == 2)
2520  {
2522 
2523  // Return the rotated vector
2524  return cross_product_2d(unit_tangential_vectors<dim>(face_no, 0));
2525  }
2526  else if (dim == 3)
2527  {
2528  return cross_product_3d(unit_tangential_vectors<dim>(face_no, 0),
2529  unit_tangential_vectors<dim>(face_no, 1));
2530  }
2531 
2532  Assert(false, ExcNotImplemented());
2533 
2534  return {};
2535 }
2536 
2537 
2538 
2539 inline unsigned int
2540 ReferenceCell::n_face_orientations(const unsigned int face_no) const
2541 {
2542  AssertIndexRange(face_no, n_faces());
2543  if (get_dimension() == 1)
2544  return 1;
2545  if (get_dimension() == 2)
2546  return 2;
2548  return 8;
2549  else if (face_reference_cell(face_no) == ReferenceCells::Triangle)
2550  return 6;
2551 
2552  Assert(false, ExcInternalError());
2554 }
2555 
2556 
2557 
2558 inline bool
2560  const unsigned int line,
2561  const unsigned char combined_face_orientation,
2562  const bool line_orientation) const
2563 {
2564  if (*this == ReferenceCells::Hexahedron)
2565  {
2566  static constexpr ::ndarray<bool, 2, 8> bool_table{
2567  {{{true, true, false, true, false, false, true, false}},
2568  {{true, true, true, false, false, false, false, true}}}};
2569 
2570  return (line_orientation ==
2571  bool_table[line / 2][combined_face_orientation]);
2572  }
2573  else
2574  // TODO: This might actually be wrong for some of the other
2575  // kinds of objects. We should check this
2576  return true;
2577 }
2578 
2579 
2580 
2581 namespace internal
2582 {
2583  template <typename T>
2585  {
2586  public:
2596  {
2599  }
2600 
2604  virtual ~NoPermutation() noexcept override = default;
2605 
2609  virtual void
2610  print_info(std::ostream &out) const override
2611  {
2612  out << '[';
2613 
2614  const unsigned int n_vertices = entity_type.n_vertices();
2615 
2616  for (unsigned int i = 0; i < n_vertices; ++i)
2617  {
2618  out << vertices_0[i];
2619  if (i + 1 != n_vertices)
2620  out << ',';
2621  }
2622 
2623  out << "] is not a valid permutation of [";
2624 
2625  for (unsigned int i = 0; i < n_vertices; ++i)
2626  {
2627  out << vertices_1[i];
2628  if (i + 1 != n_vertices)
2629  out << ',';
2630  }
2631 
2632  out << "]." << std::endl;
2633  }
2634 
2639 
2644 
2649  };
2650 } // namespace internal
2651 
2652 
2653 
2654 template <typename T, std::size_t N>
2655 inline unsigned char
2656 ReferenceCell::compute_orientation(const std::array<T, N> &vertices_0,
2657  const std::array<T, N> &vertices_1) const
2658 {
2659  Assert(N >= n_vertices(),
2660  ExcMessage("The number of array elements must be equal to or "
2661  "greater than the number of vertices of the cell "
2662  "referenced by this object."));
2663 
2664  // Call the non-deprecated function, taking care of calling it only with
2665  // those array elements that we actually care about (see the note
2666  // in the documentation about the arguments potentially being
2667  // larger arrays than necessary).
2668  return get_combined_orientation(
2669  make_array_view(vertices_0.begin(), vertices_0.begin() + n_vertices()),
2670  make_array_view(vertices_1.begin(), vertices_1.begin() + n_vertices()));
2671 }
2672 
2673 
2674 
2675 template <typename T>
2676 unsigned char
2678  const ArrayView<const T> &vertices_0,
2679  const ArrayView<const T> &vertices_1) const
2680 {
2681  Assert(vertices_0.size() == n_vertices(),
2682  ExcMessage("The number of array elements must be equal to "
2683  "the number of vertices of the cell "
2684  "referenced by this object."));
2685  Assert(vertices_1.size() == n_vertices(),
2686  ExcMessage("The number of array elements must be equal to "
2687  "the number of vertices of the cell "
2688  "referenced by this object."));
2689 
2690  const auto v0_equals = [&](const std::initializer_list<const T> &list) {
2691  Assert(list.size() == n_vertices(), ExcInternalError());
2692  return std::equal(vertices_0.begin(), vertices_0.end(), std::begin(list));
2693  };
2694 
2695  switch (this->kind)
2696  {
2697  case ReferenceCells::Line:
2698  // line_orientation=true
2699  if (v0_equals({vertices_1[0], vertices_1[1]}))
2700  return 1;
2701 
2702  // line_orientation=false
2703  if (v0_equals({vertices_1[1], vertices_1[0]}))
2704  return 0;
2705  break;
2707  // face_orientation=true, face_rotation=false, face_flip=false
2708  if (v0_equals({vertices_1[0], vertices_1[1], vertices_1[2]}))
2709  return 1;
2710 
2711  // face_orientation=true, face_rotation=true, face_flip=false
2712  if (v0_equals({vertices_1[1], vertices_1[2], vertices_1[0]}))
2713  return 3;
2714 
2715  // face_orientation=true, face_rotation=false, face_flip=true
2716  if (v0_equals({vertices_1[2], vertices_1[0], vertices_1[1]}))
2717  return 5;
2718 
2719  // face_orientation=false, face_rotation=false, face_flip=false
2720  if (v0_equals({vertices_1[0], vertices_1[2], vertices_1[1]}))
2721  return 0;
2722 
2723  // face_orientation=false, face_rotation=true, face_flip=false
2724  if (v0_equals({vertices_1[2], vertices_1[1], vertices_1[0]}))
2725  return 2;
2726 
2727  // face_orientation=false, face_rotation=false, face_flip=true
2728  if (v0_equals({vertices_1[1], vertices_1[0], vertices_1[2]}))
2729  return 4;
2730  break;
2732  // face_orientation=true, face_rotation=false, face_flip=false
2733  if (v0_equals(
2734  {vertices_1[0], vertices_1[1], vertices_1[2], vertices_1[3]}))
2735  return 1;
2736 
2737  // face_orientation=true, face_rotation=true, face_flip=false
2738  if (v0_equals(
2739  {vertices_1[2], vertices_1[0], vertices_1[3], vertices_1[1]}))
2740  return 3;
2741 
2742  // face_orientation=true, face_rotation=false, face_flip=true
2743  if (v0_equals(
2744  {vertices_1[3], vertices_1[2], vertices_1[1], vertices_1[0]}))
2745  return 5;
2746 
2747  // face_orientation=true, face_rotation=true, face_flip=true
2748  if (v0_equals(
2749  {vertices_1[1], vertices_1[3], vertices_1[0], vertices_1[2]}))
2750  return 7;
2751 
2752  // face_orientation=false, face_rotation=false, face_flip=false
2753  if (v0_equals(
2754  {vertices_1[0], vertices_1[2], vertices_1[1], vertices_1[3]}))
2755  return 0;
2756 
2757  // face_orientation=false, face_rotation=true, face_flip=false
2758  if (v0_equals(
2759  {vertices_1[2], vertices_1[3], vertices_1[0], vertices_1[1]}))
2760  return 2;
2761 
2762  // face_orientation=false, face_rotation=false, face_flip=true
2763  if (v0_equals(
2764  {vertices_1[3], vertices_1[1], vertices_1[2], vertices_1[0]}))
2765  return 4;
2766 
2767  // face_orientation=false, face_rotation=true, face_flip=true
2768  if (v0_equals(
2769  {vertices_1[1], vertices_1[0], vertices_1[3], vertices_1[2]}))
2770  return 6;
2771  break;
2772  default:
2773  Assert(false, ExcNotImplemented());
2774  }
2775 
2776  Assert(false, (internal::NoPermutation<T>(*this, vertices_0, vertices_1)));
2778 }
2779 
2780 
2781 
2782 template <typename T, std::size_t N>
2783 inline std::array<T, N>
2785  const std::array<T, N> &vertices,
2786  const unsigned int orientation) const
2787 {
2788  Assert(N >= n_vertices(),
2789  ExcMessage("The number of array elements must be equal to or "
2790  "greater than the number of vertices of the cell "
2791  "referenced by this object."));
2792 
2793  // Call the non-deprecated function, taking care of calling it only with
2794  // those array elements that we actually care about (see the note
2795  // in the documentation about the arguments potentially being
2796  // larger arrays than necessary).
2797  const auto permutation = permute_by_combined_orientation(
2798  make_array_view(vertices.begin(), vertices.begin() + n_vertices()),
2799  orientation);
2800 
2801  std::array<T, N> temp;
2802  std::copy(permutation.begin(), permutation.end(), temp.begin());
2803 
2804  return temp;
2805 }
2806 
2807 
2808 
2809 template <typename T>
2810 boost::container::small_vector<T, 8>
2813  const unsigned char orientation) const
2814 {
2815  Assert(vertices.size() == n_vertices(),
2816  ExcMessage("The number of array elements must be equal to "
2817  "the number of vertices of the cell "
2818  "referenced by this object."));
2819 
2820  switch (this->kind)
2821  {
2822  case ReferenceCells::Line:
2823  switch (orientation)
2824  {
2825  case 1:
2826  return {vertices[0], vertices[1]};
2827  case 0:
2828  return {vertices[1], vertices[0]};
2829  default:
2830  Assert(false, ExcNotImplemented());
2831  }
2832  break;
2834  switch (orientation)
2835  {
2836  case 1:
2837  return {vertices[0], vertices[1], vertices[2]};
2838  case 3:
2839  return {vertices[1], vertices[2], vertices[0]};
2840  case 5:
2841  return {vertices[2], vertices[0], vertices[1]};
2842  case 0:
2843  return {vertices[0], vertices[2], vertices[1]};
2844  case 2:
2845  return {vertices[2], vertices[1], vertices[0]};
2846  case 4:
2847  return {vertices[1], vertices[0], vertices[2]};
2848  default:
2849  Assert(false, ExcNotImplemented());
2850  }
2851  break;
2853  switch (orientation)
2854  {
2855  case 1:
2856  return {vertices[0], vertices[1], vertices[2], vertices[3]};
2857  case 3:
2858  return {vertices[2], vertices[0], vertices[3], vertices[1]};
2859  case 5:
2860  return {vertices[3], vertices[2], vertices[1], vertices[0]};
2861  case 7:
2862  return {vertices[1], vertices[3], vertices[0], vertices[2]};
2863  case 0:
2864  return {vertices[0], vertices[2], vertices[1], vertices[3]};
2865  case 2:
2866  return {vertices[2], vertices[3], vertices[0], vertices[1]};
2867  case 4:
2868  return {vertices[3], vertices[1], vertices[2], vertices[0]};
2869  case 6:
2870  return {vertices[1], vertices[0], vertices[3], vertices[2]};
2871  default:
2872  Assert(false, ExcNotImplemented());
2873  }
2874  break;
2875  default:
2876  AssertThrow(false, ExcNotImplemented());
2877  }
2878 
2879  return {};
2880 }
2881 
2882 
2883 
2884 template <>
2885 unsigned int
2886 ReferenceCell::vtk_lexicographic_to_node_index<0>(
2887  const std::array<unsigned, 0> &node_indices,
2888  const std::array<unsigned, 0> &nodes_per_direction,
2889  const bool legacy_format) const;
2890 
2891 template <>
2892 unsigned int
2893 ReferenceCell::vtk_lexicographic_to_node_index<1>(
2894  const std::array<unsigned, 1> &node_indices,
2895  const std::array<unsigned, 1> &nodes_per_direction,
2896  const bool legacy_format) const;
2897 
2898 template <>
2899 unsigned int
2900 ReferenceCell::vtk_lexicographic_to_node_index<2>(
2901  const std::array<unsigned, 2> &node_indices,
2902  const std::array<unsigned, 2> &nodes_per_direction,
2903  const bool legacy_format) const;
2904 
2905 template <>
2906 unsigned int
2907 ReferenceCell::vtk_lexicographic_to_node_index<3>(
2908  const std::array<unsigned, 3> &node_indices,
2909  const std::array<unsigned, 3> &nodes_per_direction,
2910  const bool legacy_format) const;
2911 
2913 
2914 #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:685
OutputOperator< VectorType > & operator<<(OutputOperator< VectorType > &out, unsigned int step)
Definition: operator.h:165
iterator begin() const
Definition: array_view.h:575
iterator end() const
Definition: array_view.h:584
std::size_t size() const
Definition: array_view.h:566
Abstract base class for mapping classes.
Definition: mapping.h:312
Definition: point.h:111
std::istream & operator>>(std::istream &in, Point< dim, Number > &p)
Definition: point.h:685
static 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
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
std::array< unsigned int, 2 > standard_vertex_to_face_and_vertex_index(const unsigned int vertex) const
Quadrature< dim > get_gauss_type_quadrature(const unsigned n_points_1D) const
bool standard_vs_true_line_orientation(const unsigned int line, const unsigned char face_orientation, const bool line_orientation) 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
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
const Quadrature< dim > & get_nodal_type_quadrature() 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
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_NAMESPACE_OPEN
Definition: config.h:461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
Point< 3 > vertices[4]
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1509
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1703
#define AssertIndexRange(index, range)
Definition: exceptions.h:1768
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1619
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:2394
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)
bool get_bit(const unsigned char number, const unsigned int n)
Definition: utilities.h:1549
void copy(const T *begin, const T *end, U *dest)
constexpr ReferenceCell make_reference_cell_from_int(const std::uint8_t kind)
static const unsigned int invalid_unsigned_int
Definition: types.h:206
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)