Reference documentation for deal.II version GIT dad323def1 2022-06-25 19:00:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 <iosfwd>
29 #include <string>
30 
32 
33 // Forward declarations
34 #ifndef DOXYGEN
35 template <int dim, int spacedim>
36 class Mapping;
37 
38 template <int dim>
39 class Quadrature;
40 
41 class ReferenceCell;
42 #endif
43 
44 
45 namespace internal
46 {
47  namespace ReferenceCell
48  {
62  DEAL_II_CONSTEXPR ::ReferenceCell
63  make_reference_cell_from_int(const std::uint8_t kind);
64 
65  } // namespace ReferenceCell
66 } // namespace internal
67 
68 
69 
101 {
102 public:
110  static ReferenceCell
111  n_vertices_to_type(const int dim, const unsigned int n_vertices);
112 
122  ReferenceCell();
123 
134  bool
135  is_hyper_cube() const;
136 
140  bool
141  is_simplex() const;
142 
147  unsigned int
148  get_dimension() const;
149 
163  template <int dim>
164  double
165  d_linear_shape_function(const Point<dim> &xi, const unsigned int i) const;
166 
171  template <int dim>
174  const unsigned int i) const;
175 
185  template <int dim, int spacedim = dim>
186  std::unique_ptr<Mapping<dim, spacedim>>
187  get_default_mapping(const unsigned int degree) const;
188 
200  template <int dim, int spacedim = dim>
201  const Mapping<dim, spacedim> &
203 
212  template <int dim>
214  get_gauss_type_quadrature(const unsigned n_points_1D) const;
215 
225  template <int dim>
227  get_midpoint_quadrature() const;
228 
242  template <int dim>
243  const Quadrature<dim> &
245 
260  unsigned int
261  n_vertices() const;
262 
268  vertex_indices() const;
269 
277  template <int dim>
278  Point<dim>
279  vertex(const unsigned int v) const;
280 
286  unsigned int
287  n_lines() const;
288 
294  line_indices() const;
295 
301  unsigned int
302  n_faces() const;
303 
309  face_indices() const;
310 
322  unsigned int
323  n_isotropic_children() const;
324 
330  isotropic_child_indices() const;
331 
345  face_reference_cell(const unsigned int face_no) const;
346 
397  unsigned int
398  child_cell_on_face(const unsigned int face,
399  const unsigned int subface,
400  const unsigned char face_orientation = 1) const;
401 
411  std::array<unsigned int, 2>
412  standard_vertex_to_face_and_vertex_index(const unsigned int vertex) const;
413 
423  std::array<unsigned int, 2>
424  standard_line_to_face_and_line_index(const unsigned int line) const;
425 
429  unsigned int
430  face_to_cell_lines(const unsigned int face,
431  const unsigned int line,
432  const unsigned char face_orientation) const;
433 
437  unsigned int
438  face_to_cell_vertices(const unsigned int face,
439  const unsigned int vertex,
440  const unsigned char face_orientation) const;
441 
445  unsigned int
446  standard_to_real_face_vertex(const unsigned int vertex,
447  const unsigned int face,
448  const unsigned char face_orientation) const;
449 
453  unsigned int
454  standard_to_real_face_line(const unsigned int line,
455  const unsigned int face,
456  const unsigned char face_orientation) const;
457 
468  bool
469  standard_vs_true_line_orientation(const unsigned int line,
470  const unsigned char face_orientation,
471  const bool line_orientation) const;
472 
496  double
497  volume() const;
498 
511  template <int dim>
512  Point<dim>
513  barycenter() const;
514 
534  template <int dim>
535  bool
536  contains_point(const Point<dim> &p, const double tolerance = 0) const;
537 
538  /*
539  * Return @f$i@f$-th unit tangential vector of a face of the reference cell.
540  * The vectors are arranged such that the
541  * cross product between the two vectors returns the unit normal vector.
542  *
543  * @pre @f$i@f$ must be between zero and `dim-1`.
544  */
545  template <int dim>
547  unit_tangential_vectors(const unsigned int face_no,
548  const unsigned int i) const;
549 
553  template <int dim>
555  unit_normal_vectors(const unsigned int face_no) const;
556 
561  template <typename T, std::size_t N>
562  unsigned char
563  compute_orientation(const std::array<T, N> &vertices_0,
564  const std::array<T, N> &vertices_1) const;
565 
569  template <typename T, std::size_t N>
570  std::array<T, N>
571  permute_according_orientation(const std::array<T, N> &vertices,
572  const unsigned int orientation) const;
573 
578  faces_for_given_vertex(const unsigned int vertex_index) const;
579 
592  unsigned int
593  exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const;
594 
598  unsigned int
599  exodusii_face_to_deal_face(const unsigned int face_n) const;
600 
604  unsigned int
605  unv_vertex_to_deal_vertex(const unsigned int vertex_n) const;
606 
610  unsigned int
611  vtk_linear_type() const;
612 
617  unsigned int
618  vtk_quadratic_type() const;
619 
624  unsigned int
625  vtk_lagrange_type() const;
626 
630  unsigned int
631  gmsh_element_type() const;
632 
646  std::string
647  to_string() const;
648 
652  constexpr operator std::uint8_t() const;
653 
657  constexpr bool
658  operator==(const ReferenceCell &type) const;
659 
663  constexpr bool
664  operator!=(const ReferenceCell &type) const;
665 
671  template <class Archive>
672  void
673  serialize(Archive &archive, const unsigned int /*version*/);
674 
679 private:
683  std::uint8_t kind;
684 
691  constexpr ReferenceCell(const std::uint8_t kind);
692 
699 
700  friend std::ostream &
701  operator<<(std::ostream &out, const ReferenceCell &reference_cell);
702 
703  friend std::istream &
704  operator>>(std::istream &in, ReferenceCell &reference_cell);
705 };
706 
707 
715 std::ostream &
716 operator<<(std::ostream &out, const ReferenceCell &reference_cell);
717 
725 std::istream &
726 operator>>(std::istream &in, ReferenceCell &reference_cell);
727 
728 
729 inline constexpr ReferenceCell::ReferenceCell(const std::uint8_t kind)
730  : kind(kind)
731 {}
732 
733 
734 
735 inline constexpr ReferenceCell::operator std::uint8_t() const
736 {
737  return kind;
738 }
739 
740 
741 
742 inline constexpr bool
744 {
745  return kind == type.kind;
746 }
747 
748 
749 
750 inline constexpr bool
752 {
753  return kind != type.kind;
754 }
755 
756 
757 
758 namespace internal
759 {
760  namespace ReferenceCell
761  {
762  inline DEAL_II_CONSTEXPR ::ReferenceCell
763  make_reference_cell_from_int(const std::uint8_t kind)
764  {
765  // Make sure these are the only indices from which objects can be
766  // created.
767  Assert((kind == static_cast<std::uint8_t>(-1)) || (kind < 8),
768  ExcInternalError());
769 
770  // Call the private constructor, which we can from here because this
771  // function is a 'friend'.
772  return {kind};
773  }
774  } // namespace ReferenceCell
775 } // namespace internal
776 
777 
778 
786 namespace ReferenceCells
787 {
806  static_cast<std::uint8_t>(-1));
807 
813  template <int dim>
814  constexpr const ReferenceCell &
815  get_simplex();
816 
822  template <int dim>
823  constexpr const ReferenceCell &
824  get_hypercube();
825 } // namespace ReferenceCells
826 
827 
828 
829 inline DEAL_II_CONSTEXPR
832 {}
833 
834 
835 
836 template <class Archive>
837 inline void
838 ReferenceCell::serialize(Archive &archive, const unsigned int /*version*/)
839 {
840  archive &kind;
841 }
842 
843 
844 
846 ReferenceCell::faces_for_given_vertex(const unsigned int vertex) const
847 {
848  if (*this == ReferenceCells::Line)
849  {
851  return {&GeometryInfo<2>::vertex_to_face[vertex][0], 1};
852  }
853  else if (*this == ReferenceCells::Quadrilateral)
854  {
856  return {&GeometryInfo<2>::vertex_to_face[vertex][0], 2};
857  }
858  else if (*this == ReferenceCells::Hexahedron)
859  {
861  return {&GeometryInfo<3>::vertex_to_face[vertex][0], 3};
862  }
863  else if (*this == ReferenceCells::Triangle)
864  {
866  static const ndarray<unsigned int, 3, 2> table = {
867  {{{0, 2}}, {{0, 1}}, {{1, 2}}}};
868 
869  return table[vertex];
870  }
871  else if (*this == ReferenceCells::Tetrahedron)
872  {
874  static const ndarray<unsigned int, 4, 3> table = {
875  {{{0, 1, 2}}, {{0, 1, 3}}, {{0, 2, 3}}, {{1, 2, 3}}}};
876 
877  return table[vertex];
878  }
879  else if (*this == ReferenceCells::Wedge)
880  {
882  static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 4}},
883  {{0, 2, 3}},
884  {{0, 3, 4}},
885  {{1, 2, 4}},
886  {{1, 2, 3}},
887  {{1, 3, 4}}}};
888 
889  return table[vertex];
890  }
891  else if (*this == ReferenceCells::Pyramid)
892  {
894  static const unsigned int X = numbers::invalid_unsigned_int;
895  static const ndarray<unsigned int, 5, 4> table = {{{{0, 1, 3, X}},
896  {{0, 2, 3, X}},
897  {{0, 1, 4, X}},
898  {{0, 2, 4, X}},
899  {{1, 2, 3, 4}}}};
900 
901  return {&table[vertex][0], vertex == 4 ? 4u : 3u};
902  }
903 
904  Assert(false, ExcNotImplemented());
905 
906  return {};
907 }
908 
909 
910 
911 inline bool
913 {
914  return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
916  *this == ReferenceCells::Hexahedron);
917 }
918 
919 
920 
921 inline bool
923 {
924  return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
925  *this == ReferenceCells::Triangle ||
926  *this == ReferenceCells::Tetrahedron);
927 }
928 
929 
930 
931 inline unsigned int
933 {
934  if (*this == ReferenceCells::Vertex)
935  return 0;
936  else if (*this == ReferenceCells::Line)
937  return 1;
938  else if ((*this == ReferenceCells::Triangle) ||
940  return 2;
941  else if ((*this == ReferenceCells::Tetrahedron) ||
942  (*this == ReferenceCells::Pyramid) ||
943  (*this == ReferenceCells::Wedge) ||
944  (*this == ReferenceCells::Hexahedron))
945  return 3;
946 
947  Assert(false, ExcNotImplemented());
949 }
950 
951 
952 
953 template <int dim>
956 {
957  return Quadrature<dim>(std::vector<Point<dim>>({barycenter<dim>()}),
958  std::vector<double>({volume()}));
959 }
960 
961 
962 
963 inline unsigned int
965 {
966  if (*this == ReferenceCells::Vertex)
967  return 1;
968  else if (*this == ReferenceCells::Line)
969  return 2;
970  else if (*this == ReferenceCells::Triangle)
971  return 3;
972  else if (*this == ReferenceCells::Quadrilateral)
973  return 4;
974  else if (*this == ReferenceCells::Tetrahedron)
975  return 4;
976  else if (*this == ReferenceCells::Pyramid)
977  return 5;
978  else if (*this == ReferenceCells::Wedge)
979  return 6;
980  else if (*this == ReferenceCells::Hexahedron)
981  return 8;
982 
983  Assert(false, ExcNotImplemented());
985 }
986 
987 
988 
989 inline unsigned int
991 {
992  if (*this == ReferenceCells::Vertex)
993  return 0;
994  else if (*this == ReferenceCells::Line)
995  return 1;
996  else if (*this == ReferenceCells::Triangle)
997  return 3;
998  else if (*this == ReferenceCells::Quadrilateral)
999  return 4;
1000  else if (*this == ReferenceCells::Tetrahedron)
1001  return 6;
1002  else if (*this == ReferenceCells::Pyramid)
1003  return 7;
1004  else if (*this == ReferenceCells::Wedge)
1005  return 9;
1006  else if (*this == ReferenceCells::Hexahedron)
1007  return 12;
1008 
1009  Assert(false, ExcNotImplemented());
1011 }
1012 
1013 
1014 
1015 template <int dim>
1016 Point<dim>
1017 ReferenceCell::vertex(const unsigned int v) const
1018 {
1021 
1022  if ((dim == 0) && (*this == ReferenceCells::Vertex))
1023  {
1024  return Point<dim>(0);
1025  }
1026  else if ((dim == 1) && (*this == ReferenceCells::Line))
1027  {
1028  static const Point<dim> vertices[2] = {
1029  Point<dim>(), // the origin
1030  Point<dim>::unit_vector(0) // unit point along x-axis
1031  };
1032  return vertices[v];
1033  }
1034  else if ((dim == 2) && (*this == ReferenceCells::Quadrilateral))
1035  {
1036  static const Point<dim> vertices[4] = {
1037  // First the two points on the x-axis
1038  Point<dim>(),
1040  // Then these two points shifted in the y-direction
1043  return vertices[v];
1044  }
1045  else if ((dim == 3) && (*this == ReferenceCells::Hexahedron))
1046  {
1047  static const Point<dim> vertices[8] = {
1048  // First the two points on the x-axis
1049  Point<dim>(),
1051  // Then these two points shifted in the y-direction
1054  // And now all four points shifted in the z-direction
1060  return vertices[v];
1061  }
1062  else if ((dim == 2) && (*this == ReferenceCells::Triangle))
1063  {
1064  static const Point<dim> vertices[3] = {
1065  Point<dim>(), // the origin
1066  Point<dim>::unit_vector(0), // unit point along x-axis
1067  Point<dim>::unit_vector(1) // unit point along y-axis
1068  };
1069  return vertices[v];
1070  }
1071  else if ((dim == 3) && (*this == ReferenceCells::Tetrahedron))
1072  {
1073  static const Point<dim> vertices[4] = {
1074  Point<dim>(), // the origin
1075  Point<dim>::unit_vector(0), // unit point along x-axis
1076  Point<dim>::unit_vector(1), // unit point along y-axis
1077  Point<dim>::unit_vector(2) // unit point along z-axis
1078  };
1079  return vertices[v];
1080  }
1081  else if ((dim == 3) && (*this == ReferenceCells::Pyramid))
1082  {
1083  static const Point<dim> vertices[5] = {Point<dim>{-1.0, -1.0, 0.0},
1084  Point<dim>{+1.0, -1.0, 0.0},
1085  Point<dim>{-1.0, +1.0, 0.0},
1086  Point<dim>{+1.0, +1.0, 0.0},
1087  Point<dim>{+0.0, +0.0, 1.0}};
1088  return vertices[v];
1089  }
1090  else if ((dim == 3) && (*this == ReferenceCells::Wedge))
1091  {
1092  static const Point<dim> vertices[6] = {
1093  // First the three points on the triangular base of the wedge:
1094  Point<dim>(),
1097  // And now everything shifted in the z-direction again
1101  return vertices[v];
1102  }
1103  else
1104  {
1105  Assert(false, ExcNotImplemented());
1106  return Point<dim>();
1107  }
1108 }
1109 
1110 
1111 inline unsigned int
1113 {
1114  if (*this == ReferenceCells::Vertex)
1115  return 0;
1116  else if (*this == ReferenceCells::Line)
1117  return 2;
1118  else if (*this == ReferenceCells::Triangle)
1119  return 3;
1120  else if (*this == ReferenceCells::Quadrilateral)
1121  return 4;
1122  else if (*this == ReferenceCells::Tetrahedron)
1123  return 4;
1124  else if (*this == ReferenceCells::Pyramid)
1125  return 5;
1126  else if (*this == ReferenceCells::Wedge)
1127  return 5;
1128  else if (*this == ReferenceCells::Hexahedron)
1129  return 6;
1130 
1131  Assert(false, ExcNotImplemented());
1133 }
1134 
1135 
1136 
1139 {
1140  return {0U, n_faces()};
1141 }
1142 
1143 
1144 
1145 inline unsigned int
1147 {
1148  if (*this == ReferenceCells::Vertex)
1149  return 0;
1150  else if (*this == ReferenceCells::Line)
1151  return 2;
1152  else if (*this == ReferenceCells::Triangle)
1153  return 4;
1154  else if (*this == ReferenceCells::Quadrilateral)
1155  return 4;
1156  else if (*this == ReferenceCells::Tetrahedron)
1157  return 8;
1158  else if (*this == ReferenceCells::Pyramid)
1159  {
1160  // We haven't yet decided how to refine pyramids. Update
1161  // this when we have
1162  Assert(false, ExcNotImplemented());
1164  }
1165  else if (*this == ReferenceCells::Wedge)
1166  return 8;
1167  else if (*this == ReferenceCells::Hexahedron)
1168  return 8;
1169 
1170  Assert(false, ExcNotImplemented());
1172 }
1173 
1174 
1175 
1178 {
1179  return {0U, n_isotropic_children()};
1180 }
1181 
1182 
1183 
1186 {
1187  return {0U, n_vertices()};
1188 }
1189 
1190 
1191 
1194 {
1195  return {0U, n_lines()};
1196 }
1197 
1198 
1199 
1200 inline ReferenceCell
1201 ReferenceCell::face_reference_cell(const unsigned int face_no) const
1202 {
1203  AssertIndexRange(face_no, n_faces());
1204 
1205  if (*this == ReferenceCells::Vertex)
1206  return ReferenceCells::Invalid;
1207  else if (*this == ReferenceCells::Line)
1208  return ReferenceCells::Vertex;
1209  else if (*this == ReferenceCells::Triangle)
1210  return ReferenceCells::Line;
1211  else if (*this == ReferenceCells::Quadrilateral)
1212  return ReferenceCells::Line;
1213  else if (*this == ReferenceCells::Tetrahedron)
1214  return ReferenceCells::Triangle;
1215  else if (*this == ReferenceCells::Pyramid)
1216  {
1217  if (face_no == 0)
1219  else
1220  return ReferenceCells::Triangle;
1221  }
1222  else if (*this == ReferenceCells::Wedge)
1223  {
1224  if (face_no > 1)
1226  else
1227  return ReferenceCells::Triangle;
1228  }
1229  else if (*this == ReferenceCells::Hexahedron)
1231 
1232  Assert(false, ExcNotImplemented());
1233  return ReferenceCells::Invalid;
1234 }
1235 
1236 
1237 
1238 inline unsigned int
1240  const unsigned int face,
1241  const unsigned int subface,
1242  const unsigned char face_orientation_raw) const
1243 {
1244  AssertIndexRange(face, n_faces());
1246 
1247  if (*this == ReferenceCells::Vertex)
1248  {
1249  Assert(false, ExcNotImplemented());
1250  }
1251  else if (*this == ReferenceCells::Line)
1252  {
1253  Assert(false, ExcNotImplemented());
1254  }
1255  else if (*this == ReferenceCells::Triangle)
1256  {
1257  static const ndarray<unsigned int, 3, 2> subcells = {
1258  {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
1259 
1260  return subcells[face][subface];
1261  }
1262  else if (*this == ReferenceCells::Quadrilateral)
1263  {
1264  const bool face_orientation = Utilities::get_bit(face_orientation_raw, 0);
1265  const bool face_flip = Utilities::get_bit(face_orientation_raw, 2);
1266  const bool face_rotation = Utilities::get_bit(face_orientation_raw, 1);
1267 
1270  face,
1271  subface,
1272  face_orientation,
1273  face_flip,
1274  face_rotation);
1275  }
1276  else if (*this == ReferenceCells::Tetrahedron)
1277  {
1278  Assert(false, ExcNotImplemented());
1279  }
1280  else if (*this == ReferenceCells::Pyramid)
1281  {
1282  Assert(false, ExcNotImplemented());
1283  }
1284  else if (*this == ReferenceCells::Wedge)
1285  {
1286  Assert(false, ExcNotImplemented());
1287  }
1288  else if (*this == ReferenceCells::Hexahedron)
1289  {
1290  const bool face_orientation = Utilities::get_bit(face_orientation_raw, 0);
1291  const bool face_flip = Utilities::get_bit(face_orientation_raw, 2);
1292  const bool face_rotation = Utilities::get_bit(face_orientation_raw, 1);
1293 
1296  face,
1297  subface,
1298  face_orientation,
1299  face_flip,
1300  face_rotation);
1301  }
1302 
1303  Assert(false, ExcNotImplemented());
1304  return {};
1305 }
1306 
1307 
1308 
1309 inline std::array<unsigned int, 2>
1311  const unsigned int vertex) const
1312 {
1314  // Work around a GCC warning at higher optimization levels by making all of
1315  // these tables the same size
1316  constexpr unsigned int X = numbers::invalid_unsigned_int;
1317 
1318  if (*this == ReferenceCells::Vertex)
1319  {
1320  Assert(false, ExcNotImplemented());
1321  }
1322  else if (*this == ReferenceCells::Line)
1323  {
1324  Assert(false, ExcNotImplemented());
1325  }
1326  else if (*this == ReferenceCells::Triangle)
1327  {
1328  static const ndarray<unsigned int, 6, 2> table = {
1329  {{{0, 0}}, {{0, 1}}, {{1, 1}}, {{X, X}}, {{X, X}}, {{X, X}}}};
1330 
1331  return table[vertex];
1332  }
1333  else if (*this == ReferenceCells::Quadrilateral)
1334  {
1336  }
1337  else if (*this == ReferenceCells::Tetrahedron)
1338  {
1339  static const ndarray<unsigned int, 6, 2> table = {
1340  {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 2}}, {{X, X}}, {{X, X}}}};
1341 
1342  return table[vertex];
1343  }
1344  else if (*this == ReferenceCells::Pyramid)
1345  {
1346  static const ndarray<unsigned int, 6, 2> table = {
1347  {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{0, 3}}, {{1, 2}}, {{X, X}}}};
1348 
1349  return table[vertex];
1350  }
1351  else if (*this == ReferenceCells::Wedge)
1352  {
1353  static const ndarray<unsigned int, 6, 2> table = {
1354  {{{0, 1}}, {{0, 0}}, {{0, 2}}, {{1, 0}}, {{1, 1}}, {{1, 2}}}};
1355 
1356  return table[vertex];
1357  }
1358  else if (*this == ReferenceCells::Hexahedron)
1359  {
1361  }
1362 
1363  Assert(false, ExcNotImplemented());
1364  return {};
1365 }
1366 
1367 
1368 
1369 inline std::array<unsigned int, 2>
1371  const unsigned int line) const
1372 {
1373  AssertIndexRange(line, n_lines());
1374 
1375  // start with most common cases
1376  if (*this == ReferenceCells::Hexahedron)
1377  {
1379  }
1380  else if (*this == ReferenceCells::Tetrahedron)
1381  {
1382  static const std::array<unsigned int, 2> table[6] = {
1383  {{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}};
1384 
1385  return table[line];
1386  }
1387  else if (*this == ReferenceCells::Pyramid)
1388  {
1389  static const std::array<unsigned int, 2> table[8] = {{{0, 0}},
1390  {{0, 1}},
1391  {{0, 2}},
1392  {{0, 3}},
1393  {{1, 2}},
1394  {{2, 1}},
1395  {{1, 1}},
1396  {{2, 2}}};
1397 
1398  return table[line];
1399  }
1400  else if (*this == ReferenceCells::Wedge)
1401  {
1402  static const std::array<unsigned int, 2> table[9] = {{{0, 0}},
1403  {{0, 2}},
1404  {{0, 1}},
1405  {{1, 0}},
1406  {{1, 1}},
1407  {{1, 2}},
1408  {{2, 0}},
1409  {{2, 1}},
1410  {{3, 1}}};
1411 
1412  return table[line];
1413  }
1414  else if (*this == ReferenceCells::Vertex)
1415  {
1416  Assert(false, ExcNotImplemented());
1417  }
1418  else if (*this == ReferenceCells::Line)
1419  {
1420  Assert(false, ExcNotImplemented());
1421  }
1422  else if (*this == ReferenceCells::Triangle)
1423  {
1424  Assert(false, ExcNotImplemented());
1425  }
1426  else if (*this == ReferenceCells::Quadrilateral)
1427  {
1428  Assert(false, ExcNotImplemented());
1429  }
1430 
1431  Assert(false, ExcNotImplemented());
1432  return {};
1433 }
1434 
1435 
1436 
1437 inline unsigned int
1438 ReferenceCell::face_to_cell_lines(const unsigned int face,
1439  const unsigned int line,
1440  const unsigned char face_orientation) const
1441 {
1442  AssertIndexRange(face, n_faces());
1444 
1445  static constexpr unsigned int X = numbers::invalid_unsigned_int;
1446 
1447  if (*this == ReferenceCells::Vertex)
1448  {
1449  Assert(false, ExcNotImplemented());
1450  }
1451  else if (*this == ReferenceCells::Line)
1452  {
1454  face,
1455  line,
1456  Utilities::get_bit(face_orientation, 0),
1457  Utilities::get_bit(face_orientation, 2),
1458  Utilities::get_bit(face_orientation, 1));
1459  }
1460  else if (*this == ReferenceCells::Triangle)
1461  {
1462  return face;
1463  }
1464  else if (*this == ReferenceCells::Quadrilateral)
1465  {
1467  face,
1468  line,
1469  Utilities::get_bit(face_orientation, 0),
1470  Utilities::get_bit(face_orientation, 2),
1471  Utilities::get_bit(face_orientation, 1));
1472  }
1473  else if (*this == ReferenceCells::Tetrahedron)
1474  {
1475  const static ndarray<unsigned int, 4, 3> table = {
1476  {{{0, 1, 2}}, {{0, 3, 4}}, {{2, 5, 3}}, {{1, 4, 5}}}};
1477 
1478  return table[face]
1479  [standard_to_real_face_line(line, face, face_orientation)];
1480  }
1481  else if (*this == ReferenceCells::Pyramid)
1482  {
1483  static const ndarray<unsigned int, 5, 4> table = {{{{0, 1, 2, 3}},
1484  {{0, 6, 4, X}},
1485  {{1, 5, 7, X}},
1486  {{2, 4, 5, X}},
1487  {{3, 7, 6, 2}}}};
1488 
1489  return table[face]
1490  [standard_to_real_face_line(line, face, face_orientation)];
1491  }
1492  else if (*this == ReferenceCells::Wedge)
1493  {
1494  static const ndarray<unsigned int, 5, 4> table = {{{{0, 2, 1, X}},
1495  {{3, 4, 5, X}},
1496  {{6, 7, 0, 3}},
1497  {{7, 8, 1, 4}},
1498  {{8, 6, 5, 2}}}};
1499 
1500  return table[face]
1501  [standard_to_real_face_line(line, face, face_orientation)];
1502  }
1503  else if (*this == ReferenceCells::Hexahedron)
1504  {
1506  face,
1507  line,
1508  Utilities::get_bit(face_orientation, 0),
1509  Utilities::get_bit(face_orientation, 2),
1510  Utilities::get_bit(face_orientation, 1));
1511  }
1512 
1513  Assert(false, ExcNotImplemented());
1515 }
1516 
1517 
1518 
1519 inline unsigned int
1521  const unsigned int vertex,
1522  const unsigned char face_orientation) const
1523 {
1524  AssertIndexRange(face, n_faces());
1526 
1527  if (*this == ReferenceCells::Vertex)
1528  {
1529  Assert(false, ExcNotImplemented());
1530  }
1531  else if (*this == ReferenceCells::Line)
1532  {
1534  face,
1535  vertex,
1536  Utilities::get_bit(face_orientation, 0),
1537  Utilities::get_bit(face_orientation, 2),
1538  Utilities::get_bit(face_orientation, 1));
1539  }
1540  else if (*this == ReferenceCells::Triangle)
1541  {
1542  static const ndarray<unsigned int, 3, 2> table = {
1543  {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
1544 
1545  return table[face][face_orientation != 0u ? vertex : (1 - vertex)];
1546  }
1547  else if (*this == ReferenceCells::Quadrilateral)
1548  {
1550  face,
1551  vertex,
1552  Utilities::get_bit(face_orientation, 0),
1553  Utilities::get_bit(face_orientation, 2),
1554  Utilities::get_bit(face_orientation, 1));
1555  }
1556  else if (*this == ReferenceCells::Tetrahedron)
1557  {
1558  static const ndarray<unsigned int, 4, 3> table = {
1559  {{{0, 1, 2}}, {{1, 0, 3}}, {{0, 2, 3}}, {{2, 1, 3}}}};
1560 
1561  return table[face][standard_to_real_face_vertex(
1562  vertex, face, face_orientation)];
1563  }
1564  else if (*this == ReferenceCells::Pyramid)
1565  {
1566  constexpr auto X = numbers::invalid_unsigned_int;
1567  static const ndarray<unsigned int, 5, 4> table = {{{{0, 1, 2, 3}},
1568  {{0, 2, 4, X}},
1569  {{3, 1, 4, X}},
1570  {{1, 0, 4, X}},
1571  {{2, 3, 4, X}}}};
1572 
1573  return table[face][standard_to_real_face_vertex(
1574  vertex, face, face_orientation)];
1575  }
1576  else if (*this == ReferenceCells::Wedge)
1577  {
1578  constexpr auto X = numbers::invalid_unsigned_int;
1579  static const ndarray<unsigned int, 6, 4> table = {{{{1, 0, 2, X}},
1580  {{3, 4, 5, X}},
1581  {{0, 1, 3, 4}},
1582  {{1, 2, 4, 5}},
1583  {{2, 0, 5, 3}}}};
1584 
1585  return table[face][standard_to_real_face_vertex(
1586  vertex, face, face_orientation)];
1587  }
1588  else if (*this == ReferenceCells::Hexahedron)
1589  {
1591  face,
1592  vertex,
1593  Utilities::get_bit(face_orientation, 0),
1594  Utilities::get_bit(face_orientation, 2),
1595  Utilities::get_bit(face_orientation, 1));
1596  }
1597 
1598  Assert(false, ExcNotImplemented());
1600 }
1601 
1602 
1603 
1604 inline unsigned int
1606  const unsigned int vertex,
1607  const unsigned int face,
1608  const unsigned char face_orientation) const
1609 {
1610  AssertIndexRange(face, n_faces());
1612 
1613  if (*this == ReferenceCells::Quadrilateral ||
1614  *this == ReferenceCells::Triangle)
1615  {
1616  static const ndarray<unsigned int, 2, 2> table = {{{{1, 0}}, {{0, 1}}}};
1617 
1618  return table[face_orientation][vertex];
1619  }
1620  else if (*this == ReferenceCells::Hexahedron)
1621  {
1622  static const ndarray<unsigned int, 8, 4> table = {{{{0, 2, 1, 3}},
1623  {{0, 1, 2, 3}},
1624  {{2, 3, 0, 1}},
1625  {{2, 0, 3, 1}},
1626  {{3, 1, 2, 0}},
1627  {{3, 2, 1, 0}},
1628  {{1, 0, 3, 2}},
1629  {{1, 3, 0, 2}}}};
1630  return table[face_orientation][vertex];
1631  }
1632  else if (*this == ReferenceCells::Tetrahedron)
1633  {
1634  static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 1}},
1635  {{0, 1, 2}},
1636  {{2, 1, 0}},
1637  {{1, 2, 0}},
1638  {{1, 0, 2}},
1639  {{2, 0, 1}}}};
1640 
1641  return table[face_orientation][vertex];
1642  }
1643  else if (*this == ReferenceCells::Pyramid)
1644  {
1645  if (face == 0) // The quadrilateral face
1646  {
1648  vertex,
1649  Utilities::get_bit(face_orientation, 0),
1650  Utilities::get_bit(face_orientation, 2),
1651  Utilities::get_bit(face_orientation, 1));
1652  }
1653  else // One of the triangular faces
1654  {
1655  static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 1}},
1656  {{0, 1, 2}},
1657  {{2, 1, 0}},
1658  {{1, 2, 0}},
1659  {{1, 0, 2}},
1660  {{2, 0, 1}}}};
1661 
1662  return table[face_orientation][vertex];
1663  }
1664  }
1665  else if (*this == ReferenceCells::Wedge)
1666  {
1667  if (face > 1) // One of the quadrilateral faces
1668  {
1670  vertex,
1671  Utilities::get_bit(face_orientation, 0),
1672  Utilities::get_bit(face_orientation, 2),
1673  Utilities::get_bit(face_orientation, 1));
1674  }
1675  else // One of the triangular faces
1676  {
1677  static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 1}},
1678  {{0, 1, 2}},
1679  {{2, 1, 0}},
1680  {{1, 2, 0}},
1681  {{1, 0, 2}},
1682  {{2, 0, 1}}}};
1683 
1684  return table[face_orientation][vertex];
1685  }
1686  }
1687  else if (*this == ReferenceCells::Vertex)
1688  {
1689  Assert(false, ExcNotImplemented());
1690  }
1691  else if (*this == ReferenceCells::Line)
1692  {
1693  Assert(false, ExcNotImplemented());
1694  }
1695 
1696  Assert(false, ExcNotImplemented());
1698 }
1699 
1700 
1701 
1702 inline unsigned int
1704  const unsigned int line,
1705  const unsigned int face,
1706  const unsigned char face_orientation) const
1707 {
1708  AssertIndexRange(face, n_faces());
1710 
1711  // start with the most common cases
1712  if (*this == ReferenceCells::Hexahedron)
1713  {
1714  static const ndarray<unsigned int, 8, 4> table = {{{{2, 3, 0, 1}},
1715  {{0, 1, 2, 3}},
1716  {{0, 1, 3, 2}},
1717  {{3, 2, 0, 1}},
1718  {{3, 2, 1, 0}},
1719  {{1, 0, 3, 2}},
1720  {{1, 0, 2, 3}},
1721  {{2, 3, 1, 0}}}};
1722  return table[face_orientation][line];
1723  }
1724  else if (*this == ReferenceCells::Tetrahedron)
1725  {
1726  static const ndarray<unsigned int, 6, 3> table = {{{{2, 1, 0}},
1727  {{0, 1, 2}},
1728  {{1, 0, 2}},
1729  {{1, 2, 0}},
1730  {{0, 2, 1}},
1731  {{2, 0, 1}}}};
1732 
1733  return table[face_orientation][line];
1734  }
1735  else if (*this == ReferenceCells::Pyramid)
1736  {
1737  if (face == 0) // The quadrilateral face
1738  {
1740  line,
1741  Utilities::get_bit(face_orientation, 0),
1742  Utilities::get_bit(face_orientation, 2),
1743  Utilities::get_bit(face_orientation, 1));
1744  }
1745  else // One of the triangular faces
1746  {
1747  static const ndarray<unsigned int, 6, 3> table = {{{{2, 1, 0}},
1748  {{0, 1, 2}},
1749  {{1, 0, 2}},
1750  {{1, 2, 0}},
1751  {{0, 2, 1}},
1752  {{2, 0, 1}}}};
1753 
1754  return table[face_orientation][line];
1755  }
1756  }
1757  else if (*this == ReferenceCells::Wedge)
1758  {
1759  if (face > 1) // One of the quadrilateral faces
1760  {
1762  line,
1763  Utilities::get_bit(face_orientation, 0),
1764  Utilities::get_bit(face_orientation, 2),
1765  Utilities::get_bit(face_orientation, 1));
1766  }
1767  else // One of the triangular faces
1768  {
1769  static const ndarray<unsigned int, 6, 3> table = {{{{2, 1, 0}},
1770  {{0, 1, 2}},
1771  {{1, 0, 2}},
1772  {{1, 2, 0}},
1773  {{0, 2, 1}},
1774  {{2, 0, 1}}}};
1775 
1776  return table[face_orientation][line];
1777  }
1778  }
1779  else if (*this == ReferenceCells::Vertex)
1780  {
1781  Assert(false, ExcNotImplemented());
1782  }
1783  else if (*this == ReferenceCells::Line)
1784  {
1785  Assert(false, ExcNotImplemented());
1786  }
1787  else if (*this == ReferenceCells::Triangle)
1788  {
1789  Assert(false, ExcNotImplemented());
1790  }
1791  else if (*this == ReferenceCells::Quadrilateral)
1792  {
1793  Assert(false, ExcNotImplemented());
1794  }
1795 
1796  Assert(false, ExcNotImplemented());
1798 }
1799 
1800 
1801 
1802 namespace ReferenceCells
1803 {
1804  template <int dim>
1805  inline constexpr const ReferenceCell &
1807  {
1808  switch (dim)
1809  {
1810  case 0:
1811  return ReferenceCells::Vertex;
1812  case 1:
1813  return ReferenceCells::Line;
1814  case 2:
1815  return ReferenceCells::Triangle;
1816  case 3:
1818  default:
1819  Assert(false, ExcNotImplemented());
1820  return ReferenceCells::Invalid;
1821  }
1822  }
1823 
1824 
1825 
1826  template <int dim>
1827  inline constexpr const ReferenceCell &
1829  {
1830  switch (dim)
1831  {
1832  case 0:
1833  return ReferenceCells::Vertex;
1834  case 1:
1835  return ReferenceCells::Line;
1836  case 2:
1838  case 3:
1840  default:
1841  Assert(false, ExcNotImplemented());
1842  return ReferenceCells::Invalid;
1843  }
1844  }
1845 } // namespace ReferenceCells
1846 
1847 
1848 inline ReferenceCell
1849 ReferenceCell::n_vertices_to_type(const int dim, const unsigned int n_vertices)
1850 {
1851  AssertIndexRange(dim, 4);
1853 
1854  const auto X = ReferenceCells::Invalid;
1855  static const ndarray<ReferenceCell, 4, 9> table = {
1856  {// dim 0
1857  {{X, ReferenceCells::Vertex, X, X, X, X, X, X, X}},
1858  // dim 1
1859  {{X, X, ReferenceCells::Line, X, X, X, X, X, X}},
1860  // dim 2
1861  {{X,
1862  X,
1863  X,
1866  X,
1867  X,
1868  X,
1869  X}},
1870  // dim 3
1871  {{X,
1872  X,
1873  X,
1874  X,
1878  X,
1880  Assert(table[dim][n_vertices] != ReferenceCells::Invalid,
1881  ExcMessage("The combination of dim = " + std::to_string(dim) +
1882  " and n_vertices = " + std::to_string(n_vertices) +
1883  " does not correspond to a known reference cell type."));
1884  return table[dim][n_vertices];
1885 }
1886 
1887 
1888 
1889 template <int dim>
1890 inline double
1892  const unsigned int i) const
1893 {
1895  if (*this == ReferenceCells::get_hypercube<dim>())
1897 
1898  if (*this ==
1899  ReferenceCells::Triangle) // see also
1900  // BarycentricPolynomials<2>::compute_value
1901  {
1902  switch (i)
1903  {
1904  case 0:
1905  return 1.0 - xi[std::min(0, dim - 1)] - xi[std::min(1, dim - 1)];
1906  case 1:
1907  return xi[std::min(0, dim - 1)];
1908  case 2:
1909  return xi[std::min(1, dim - 1)];
1910  }
1911  }
1912 
1913  if (*this ==
1914  ReferenceCells::Tetrahedron) // see also
1915  // BarycentricPolynomials<3>::compute_value
1916  {
1917  switch (i)
1918  {
1919  case 0:
1920  return 1.0 - xi[std::min(0, dim - 1)] - xi[std::min(1, dim - 1)] -
1921  xi[std::min(2, dim - 1)];
1922  case 1:
1923  return xi[std::min(0, dim - 1)];
1924  case 2:
1925  return xi[std::min(1, dim - 1)];
1926  case 3:
1927  return xi[std::min(2, dim - 1)];
1928  }
1929  }
1930 
1931  if (*this ==
1932  ReferenceCells::Wedge) // see also
1933  // ScalarLagrangePolynomialWedge::compute_value
1934  {
1936  .d_linear_shape_function<2>(Point<2>(xi[std::min(0, dim - 1)],
1937  xi[std::min(1, dim - 1)]),
1938  i % 3) *
1940  .d_linear_shape_function<1>(Point<1>(xi[std::min(2, dim - 1)]),
1941  i / 3);
1942  }
1943 
1944  if (*this ==
1945  ReferenceCells::Pyramid) // see also
1946  // ScalarLagrangePolynomialPyramid::compute_value
1947  {
1948  const double Q14 = 0.25;
1949  double ration;
1950 
1951  const double r = xi[std::min(0, dim - 1)];
1952  const double s = xi[std::min(1, dim - 1)];
1953  const double t = xi[std::min(2, dim - 1)];
1954 
1955  if (fabs(t - 1.0) > 1.0e-14)
1956  {
1957  ration = (r * s * t) / (1.0 - t);
1958  }
1959  else
1960  {
1961  ration = 0.0;
1962  }
1963 
1964  if (i == 0)
1965  return Q14 * ((1.0 - r) * (1.0 - s) - t + ration);
1966  if (i == 1)
1967  return Q14 * ((1.0 + r) * (1.0 - s) - t - ration);
1968  if (i == 2)
1969  return Q14 * ((1.0 - r) * (1.0 + s) - t - ration);
1970  if (i == 3)
1971  return Q14 * ((1.0 + r) * (1.0 + s) - t + ration);
1972  else
1973  return t;
1974  }
1975 
1976  Assert(false, ExcNotImplemented());
1977 
1978  return 0.0;
1979 }
1980 
1981 
1982 
1983 template <int dim>
1984 inline Tensor<1, dim>
1986  const unsigned int i) const
1987 {
1989  if (*this == ReferenceCells::get_hypercube<dim>())
1991 
1992  if (*this ==
1993  ReferenceCells::Triangle) // see also
1994  // BarycentricPolynomials<2>::compute_grad
1995  {
1996  switch (i)
1997  {
1998  case 0:
1999  return Point<dim>(-1.0, -1.0);
2000  case 1:
2001  return Point<dim>(+1.0, +0.0);
2002  case 2:
2003  return Point<dim>(+0.0, +1.0);
2004  }
2005  }
2006 
2007  Assert(false, ExcNotImplemented());
2008 
2009  return Point<dim>(+0.0, +0.0, +0.0);
2010 }
2011 
2012 
2013 
2014 inline double
2016 {
2017  if (*this == ReferenceCells::Vertex)
2018  return 0;
2019  else if (*this == ReferenceCells::Line)
2020  return 1;
2021  else if (*this == ReferenceCells::Triangle)
2022  return 1. / 2.;
2023  else if (*this == ReferenceCells::Quadrilateral)
2024  return 1;
2025  else if (*this == ReferenceCells::Tetrahedron)
2026  return 1. / 6.;
2027  else if (*this == ReferenceCells::Wedge)
2028  return 1. / 2.;
2029  else if (*this == ReferenceCells::Pyramid)
2030  return 4. / 3.;
2031  else if (*this == ReferenceCells::Hexahedron)
2032  return 1;
2033 
2034  Assert(false, ExcNotImplemented());
2035  return 0.0;
2036 }
2037 
2038 
2039 
2040 template <int dim>
2041 inline Point<dim>
2043 {
2045 
2046  if (*this == ReferenceCells::Vertex)
2047  return Point<dim>();
2048  else if (*this == ReferenceCells::Line)
2049  return Point<dim>(1. / 2.);
2050  else if (*this == ReferenceCells::Triangle)
2051  return Point<dim>(1. / 3., 1. / 3.);
2052  else if (*this == ReferenceCells::Quadrilateral)
2053  return Point<dim>(1. / 2., 1. / 2.);
2054  else if (*this == ReferenceCells::Tetrahedron)
2055  return Point<dim>(1. / 4., 1. / 4., 1. / 4.);
2056  else if (*this == ReferenceCells::Wedge)
2057  return Point<dim>(1. / 3, 1. / 3, 1. / 2.);
2058  else if (*this == ReferenceCells::Pyramid)
2059  return Point<dim>(0, 0, 1. / 4.);
2060  else if (*this == ReferenceCells::Hexahedron)
2061  return Point<dim>(1. / 2., 1. / 2., 1. / 2.);
2062 
2063  Assert(false, ExcNotImplemented());
2064  return Point<dim>();
2065 }
2066 
2067 
2068 
2069 template <int dim>
2070 inline bool
2071 ReferenceCell::contains_point(const Point<dim> &p, const double tolerance) const
2072 {
2074 
2075  if (*this == ReferenceCells::Vertex)
2076  {
2077  // Vertices are special cases in that they do not actually
2078  // have coordinates. Error out if this function is called
2079  // with a vertex:
2080  Assert(false,
2081  ExcMessage("Vertices are zero-dimensional objects and "
2082  "as a consequence have no coordinates. You "
2083  "cannot meaningfully ask whether a point is "
2084  "inside a vertex (within a certain tolerance) "
2085  "without coordinate values."));
2086  return false;
2087  }
2088  else if (*this == ReferenceCells::get_hypercube<dim>())
2089  {
2090  for (unsigned int d = 0; d < dim; ++d)
2091  if ((p[d] < -tolerance) || (p[d] > 1 + tolerance))
2092  return false;
2093  return true;
2094  }
2095  else if (*this == ReferenceCells::get_simplex<dim>())
2096  {
2097  // First make sure that we are in the first quadrant or octant
2098  for (unsigned int d = 0; d < dim; ++d)
2099  if (p[d] < -tolerance)
2100  return false;
2101 
2102  // Now we also need to make sure that we are below the diagonal line
2103  // or plane that delineates the simplex. This diagonal is given by
2104  // sum(p[d])<=1, and a diagonal a distance eps away is given by
2105  // sum(p[d])<=1+eps*sqrt(d). (For example, the point at (1,1) is a
2106  // distance of 1/sqrt(2) away from the diagonal. That is, its
2107  // sum satisfies
2108  // sum(p[d]) = 2 <= 1 + (1/sqrt(2)) * sqrt(2)
2109  // in other words, it satisfies the predicate with eps=1/sqrt(2).)
2110  double sum = 0;
2111  for (unsigned int d = 0; d < dim; ++d)
2112  sum += p[d];
2113  return (sum <= 1 + tolerance * std::sqrt(1. * dim));
2114  }
2115  else if (*this == ReferenceCells::Wedge)
2116  {
2117  Assert(false, ExcNotImplemented());
2118  }
2119  else if (*this == ReferenceCells::Pyramid)
2120  {
2121  Assert(false, ExcNotImplemented());
2122  }
2123 
2124  Assert(false, ExcNotImplemented());
2125 
2126  return false;
2127 }
2128 
2129 
2130 
2131 template <int dim>
2132 inline Tensor<1, dim>
2133 ReferenceCell::unit_tangential_vectors(const unsigned int face_no,
2134  const unsigned int i) const
2135 {
2137  AssertIndexRange(i, dim - 1);
2138 
2139  if (*this == ReferenceCells::get_hypercube<dim>())
2140  {
2142  return GeometryInfo<dim>::unit_tangential_vectors[face_no][i];
2143  }
2144  else if (*this == ReferenceCells::Triangle)
2145  {
2146  AssertIndexRange(face_no, 3);
2147  static const std::array<Tensor<1, dim>, 3> table = {
2148  {Point<dim>(1, 0),
2149  Point<dim>(-std::sqrt(0.5), +std::sqrt(0.5)),
2150  Point<dim>(0, -1)}};
2151 
2152  return table[face_no];
2153  }
2154  else if (*this == ReferenceCells::Tetrahedron)
2155  {
2156  AssertIndexRange(face_no, 4);
2157  static const ndarray<Tensor<1, dim>, 4, 2> table = {
2158  {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
2159  {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
2160  {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
2161  {{Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
2162  +std::pow(1.0 / 3.0, 1.0 / 4.0),
2163  0),
2164  Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
2165  0,
2166  +std::pow(1.0 / 3.0, 1.0 / 4.0))}}}};
2167 
2168  return table[face_no][i];
2169  }
2170  else if (*this == ReferenceCells::Wedge)
2171  {
2172  AssertIndexRange(face_no, 5);
2173  static const ndarray<Tensor<1, dim>, 5, 2> table = {
2174  {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
2175  {{Point<dim>(1, 0, 0), Point<dim>(0, 1, 0)}},
2176  {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
2177  {{Point<dim>(-1 / std::sqrt(2.0), +1 / std::sqrt(2.0), 0),
2178  Point<dim>(0, 0, 1)}},
2179  {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}}}};
2180 
2181  return table[face_no][i];
2182  }
2183  else if (*this == ReferenceCells::Pyramid)
2184  {
2185  AssertIndexRange(face_no, 5);
2186  static const ndarray<Tensor<1, dim>, 5, 2> table = {
2187  {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
2188  {{Point<dim>(+1.0 / sqrt(2.0), 0, +1.0 / sqrt(2.0)),
2189  Point<dim>(0, 1, 0)}},
2190  {{Point<dim>(+1.0 / sqrt(2.0), 0, -1.0 / sqrt(2.0)),
2191  Point<dim>(0, 1, 0)}},
2192  {{Point<dim>(1, 0, 0),
2193  Point<dim>(0, +1.0 / sqrt(2.0), +1.0 / sqrt(2.0))}},
2194  {{Point<dim>(1, 0, 0),
2195  Point<dim>(0, +1.0 / sqrt(2.0), -1.0 / sqrt(2.0))}}}};
2196 
2197  return table[face_no][i];
2198  }
2199 
2200  Assert(false, ExcNotImplemented());
2201 
2202  return {};
2203 }
2204 
2205 
2206 
2207 template <int dim>
2208 inline Tensor<1, dim>
2209 ReferenceCell::unit_normal_vectors(const unsigned int face_no) const
2210 {
2211  AssertDimension(dim, this->get_dimension());
2212 
2213  if (is_hyper_cube())
2214  {
2216  return GeometryInfo<dim>::unit_normal_vector[face_no];
2217  }
2218  else if (dim == 2)
2219  {
2221 
2222  // Return the rotated vector
2223  return cross_product_2d(unit_tangential_vectors<dim>(face_no, 0));
2224  }
2225  else if (dim == 3)
2226  {
2227  return cross_product_3d(unit_tangential_vectors<dim>(face_no, 0),
2228  unit_tangential_vectors<dim>(face_no, 1));
2229  }
2230 
2231  Assert(false, ExcNotImplemented());
2232 
2233  return {};
2234 }
2235 
2236 
2237 
2238 inline bool
2240  const unsigned int line,
2241  const unsigned char face_orientation_raw,
2242  const bool line_orientation) const
2243 {
2244  if (*this == ReferenceCells::Hexahedron)
2245  {
2246  static const bool bool_table[2][2][2][2] = {
2247  {{{true, false}, // lines 0/1, face_orientation=false,
2248  // face_flip=false, face_rotation=false and true
2249  {false, true}}, // lines 0/1, face_orientation=false,
2250  // face_flip=true, face_rotation=false and true
2251  {{true, true}, // lines 0/1, face_orientation=true,
2252  // face_flip=false, face_rotation=false and true
2253  {false, false}}}, // lines 0/1, face_orientation=true,
2254  // face_flip=true, face_rotation=false and true
2255 
2256  {{{true, true}, // lines 2/3 ...
2257  {false, false}},
2258  {{true, false}, {false, true}}}};
2259 
2260  const bool face_orientation = Utilities::get_bit(face_orientation_raw, 0);
2261  const bool face_flip = Utilities::get_bit(face_orientation_raw, 2);
2262  const bool face_rotation = Utilities::get_bit(face_orientation_raw, 1);
2263 
2264  return (line_orientation ==
2265  bool_table[line / 2][face_orientation][face_flip][face_rotation]);
2266  }
2267  else
2268  // TODO: This might actually be wrong for some of the other
2269  // kinds of objects. We should check this
2270  return true;
2271 }
2272 
2273 
2274 
2275 namespace internal
2276 {
2277  template <typename T, std::size_t N>
2279  {
2280  public:
2284  NoPermutation(const ::ReferenceCell &entity_type,
2285  const std::array<T, N> & vertices_0,
2286  const std::array<T, N> & vertices_1)
2290  {}
2291 
2295  virtual ~NoPermutation() noexcept override = default;
2296 
2300  virtual void
2301  print_info(std::ostream &out) const override
2302  {
2303  out << '[';
2304 
2305  const unsigned int n_vertices = entity_type.n_vertices();
2306 
2307  for (unsigned int i = 0; i < n_vertices; ++i)
2308  {
2309  out << vertices_0[i];
2310  if (i + 1 != n_vertices)
2311  out << ',';
2312  }
2313 
2314  out << "] is not a permutation of [";
2315 
2316  for (unsigned int i = 0; i < n_vertices; ++i)
2317  {
2318  out << vertices_1[i];
2319  if (i + 1 != n_vertices)
2320  out << ',';
2321  }
2322 
2323  out << "]." << std::endl;
2324  }
2325 
2329  const ::ReferenceCell entity_type;
2330 
2334  const std::array<T, N> vertices_0;
2335 
2339  const std::array<T, N> vertices_1;
2340  };
2341 } // namespace internal
2342 
2343 
2344 
2345 template <typename T, std::size_t N>
2346 inline unsigned char
2347 ReferenceCell::compute_orientation(const std::array<T, N> &vertices_0,
2348  const std::array<T, N> &vertices_1) const
2349 {
2350  AssertIndexRange(n_vertices(), N + 1);
2351  if (*this == ReferenceCells::Line)
2352  {
2353  const std::array<T, 2> i{{vertices_0[0], vertices_0[1]}};
2354  const std::array<T, 2> j{{vertices_1[0], vertices_1[1]}};
2355 
2356  // line_orientation=true
2357  if (i == std::array<T, 2>{{j[0], j[1]}})
2358  return 1;
2359 
2360  // line_orientation=false
2361  if (i == std::array<T, 2>{{j[1], j[0]}})
2362  return 0;
2363  }
2364  else if (*this == ReferenceCells::Triangle)
2365  {
2366  const std::array<T, 3> i{{vertices_0[0], vertices_0[1], vertices_0[2]}};
2367  const std::array<T, 3> j{{vertices_1[0], vertices_1[1], vertices_1[2]}};
2368 
2369  // face_orientation=true, face_rotation=false, face_flip=false
2370  if (i == std::array<T, 3>{{j[0], j[1], j[2]}})
2371  return 1;
2372 
2373  // face_orientation=true, face_rotation=true, face_flip=false
2374  if (i == std::array<T, 3>{{j[1], j[2], j[0]}})
2375  return 3;
2376 
2377  // face_orientation=true, face_rotation=false, face_flip=true
2378  if (i == std::array<T, 3>{{j[2], j[0], j[1]}})
2379  return 5;
2380 
2381  // face_orientation=false, face_rotation=false, face_flip=false
2382  if (i == std::array<T, 3>{{j[0], j[2], j[1]}})
2383  return 0;
2384 
2385  // face_orientation=false, face_rotation=true, face_flip=false
2386  if (i == std::array<T, 3>{{j[2], j[1], j[0]}})
2387  return 2;
2388 
2389  // face_orientation=false, face_rotation=false, face_flip=true
2390  if (i == std::array<T, 3>{{j[1], j[0], j[2]}})
2391  return 4;
2392  }
2393  else if (*this == ReferenceCells::Quadrilateral)
2394  {
2395  const std::array<T, 4> i{
2396  {vertices_0[0], vertices_0[1], vertices_0[2], vertices_0[3]}};
2397  const std::array<T, 4> j{
2398  {vertices_1[0], vertices_1[1], vertices_1[2], vertices_1[3]}};
2399 
2400  // face_orientation=true, face_rotation=false, face_flip=false
2401  if (i == std::array<T, 4>{{j[0], j[1], j[2], j[3]}})
2402  return 1;
2403 
2404  // face_orientation=true, face_rotation=true, face_flip=false
2405  if (i == std::array<T, 4>{{j[2], j[0], j[3], j[1]}})
2406  return 3;
2407 
2408  // face_orientation=true, face_rotation=false, face_flip=true
2409  if (i == std::array<T, 4>{{j[3], j[2], j[1], j[0]}})
2410  return 5;
2411 
2412  // face_orientation=true, face_rotation=true, face_flip=true
2413  if (i == std::array<T, 4>{{j[1], j[3], j[0], j[2]}})
2414  return 7;
2415 
2416  // face_orientation=false, face_rotation=false, face_flip=false
2417  if (i == std::array<T, 4>{{j[0], j[2], j[1], j[3]}})
2418  return 0;
2419 
2420  // face_orientation=false, face_rotation=true, face_flip=false
2421  if (i == std::array<T, 4>{{j[2], j[3], j[0], j[1]}})
2422  return 2;
2423 
2424  // face_orientation=false, face_rotation=false, face_flip=true
2425  if (i == std::array<T, 4>{{j[3], j[1], j[2], j[0]}})
2426  return 4;
2427 
2428  // face_orientation=false, face_rotation=true, face_flip=true
2429  if (i == std::array<T, 4>{{j[1], j[0], j[3], j[2]}})
2430  return 6;
2431  }
2432 
2433  Assert(false, (internal::NoPermutation<T, N>(*this, vertices_0, vertices_1)));
2434 
2435  return -1;
2436 }
2437 
2438 
2439 
2440 template <typename T, std::size_t N>
2441 inline std::array<T, N>
2443  const std::array<T, N> &vertices,
2444  const unsigned int orientation) const
2445 {
2446  std::array<T, 4> temp;
2447 
2448  if (*this == ReferenceCells::Line)
2449  {
2450  switch (orientation)
2451  {
2452  case 1:
2453  temp = {{vertices[0], vertices[1]}};
2454  break;
2455  case 0:
2456  temp = {{vertices[1], vertices[0]}};
2457  break;
2458  default:
2459  Assert(false, ExcNotImplemented());
2460  }
2461  }
2462  else if (*this == ReferenceCells::Triangle)
2463  {
2464  switch (orientation)
2465  {
2466  case 1:
2467  temp = {{vertices[0], vertices[1], vertices[2]}};
2468  break;
2469  case 3:
2470  temp = {{vertices[1], vertices[2], vertices[0]}};
2471  break;
2472  case 5:
2473  temp = {{vertices[2], vertices[0], vertices[1]}};
2474  break;
2475  case 0:
2476  temp = {{vertices[0], vertices[2], vertices[1]}};
2477  break;
2478  case 2:
2479  temp = {{vertices[2], vertices[1], vertices[0]}};
2480  break;
2481  case 4:
2482  temp = {{vertices[1], vertices[0], vertices[2]}};
2483  break;
2484  default:
2485  Assert(false, ExcNotImplemented());
2486  }
2487  }
2488  else if (*this == ReferenceCells::Quadrilateral)
2489  {
2490  switch (orientation)
2491  {
2492  case 1:
2493  temp = {{vertices[0], vertices[1], vertices[2], vertices[3]}};
2494  break;
2495  case 3:
2496  temp = {{vertices[2], vertices[0], vertices[3], vertices[1]}};
2497  break;
2498  case 5:
2499  temp = {{vertices[3], vertices[2], vertices[1], vertices[0]}};
2500  break;
2501  case 7:
2502  temp = {{vertices[1], vertices[3], vertices[0], vertices[2]}};
2503  break;
2504  case 0:
2505  temp = {{vertices[0], vertices[2], vertices[1], vertices[3]}};
2506  break;
2507  case 2:
2508  temp = {{vertices[2], vertices[3], vertices[0], vertices[1]}};
2509  break;
2510  case 4:
2511  temp = {{vertices[3], vertices[1], vertices[2], vertices[0]}};
2512  break;
2513  case 6:
2514  temp = {{vertices[1], vertices[0], vertices[3], vertices[2]}};
2515  break;
2516  default:
2517  Assert(false, ExcNotImplemented());
2518  }
2519  }
2520  else
2521  {
2522  AssertThrow(false, ExcNotImplemented());
2523  }
2524 
2525  std::array<T, N> temp_;
2526  std::copy_n(temp.begin(), N, temp_.begin());
2527 
2528  return temp_;
2529 }
2530 
2531 
2533 
2534 #endif
OutputOperator< VectorType > & operator<<(OutputOperator< VectorType > &out, unsigned int step)
Definition: operator.h:165
Abstract base class for mapping classes.
Definition: mapping.h:311
Definition: point.h:111
std::istream & operator>>(std::istream &in, Point< dim, Number > &p)
Definition: point.h:688
static Point< dim, Number > unit_vector(const unsigned int i)
ArrayView< const unsigned int > faces_for_given_vertex(const unsigned int vertex_index) const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices() const
unsigned int n_vertices() 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
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
unsigned int child_cell_on_face(const unsigned int face, const unsigned int subface, const unsigned char face_orientation=1) const
Tensor< 1, dim > unit_tangential_vectors(const unsigned int face_no, const unsigned int i) 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
double volume() const
unsigned int n_faces() 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
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
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
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
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
Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i) const
const ::ReferenceCell entity_type
const std::array< T, N > vertices_0
const std::array< T, N > vertices_1
virtual void print_info(std::ostream &out) const override
virtual ~NoPermutation() noexcept override=default
NoPermutation(const ::ReferenceCell &entity_type, const std::array< T, N > &vertices_0, const std::array< T, N > &vertices_1)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_CONSTEXPR
Definition: config.h:177
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Point< 3 > vertices[4]
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1473
std::string to_string(const T &t)
Definition: patterns.h:2403
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
Expression fabs(const Expression &x)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
static const char U
static const char N
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()
T sum(const T &t, const MPI_Comm &mpi_communicator)
bool get_bit(const unsigned char number, const unsigned int n)
Definition: utilities.h:1707
constexpr ::ReferenceCell make_reference_cell_from_int(const std::uint8_t kind)
static const unsigned int invalid_unsigned_int
Definition: types.h:201
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 standard_to_real_face_vertex(const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
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)
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2635
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2660