Reference documentation for deal.II version Git 500a7ed831 2022-01-17 20:04:17 -0700
\(\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 - 2021 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 
91 {
92 public:
100  static ReferenceCell
101  n_vertices_to_type(const int dim, const unsigned int n_vertices);
102 
112  ReferenceCell();
113 
124  bool
125  is_hyper_cube() const;
126 
130  bool
131  is_simplex() const;
132 
137  unsigned int
138  get_dimension() const;
139 
153  template <int dim>
154  double
155  d_linear_shape_function(const Point<dim> &xi, const unsigned int i) const;
156 
161  template <int dim>
163  d_linear_shape_function_gradient(const Point<dim> & xi,
164  const unsigned int i) const;
165 
175  template <int dim, int spacedim = dim>
176  std::unique_ptr<Mapping<dim, spacedim>>
177  get_default_mapping(const unsigned int degree) const;
178 
190  template <int dim, int spacedim = dim>
191  const Mapping<dim, spacedim> &
193 
202  template <int dim>
204  get_gauss_type_quadrature(const unsigned n_points_1D) const;
205 
212  template <int dim>
213  const Quadrature<dim> &
214  get_nodal_type_quadrature() const;
215 
230  unsigned int
231  n_vertices() const;
232 
238  vertex_indices() const;
239 
247  template <int dim>
248  Point<dim>
249  vertex(const unsigned int v) const;
250 
256  unsigned int
257  n_lines() const;
258 
264  line_indices() const;
265 
271  unsigned int
272  n_faces() const;
273 
279  face_indices() const;
280 
293  ReferenceCell
294  face_reference_cell(const unsigned int face_no) const;
295 
346  unsigned int
347  child_cell_on_face(const unsigned int face_n,
348  const unsigned int subface_n,
349  const unsigned char face_orientation = 1) const;
350 
360  std::array<unsigned int, 2>
361  standard_vertex_to_face_and_vertex_index(const unsigned int vertex) const;
362 
372  std::array<unsigned int, 2>
373  standard_line_to_face_and_line_index(const unsigned int line) const;
374 
378  unsigned int
379  face_to_cell_lines(const unsigned int face,
380  const unsigned int line,
381  const unsigned char face_orientation) const;
382 
386  unsigned int
387  face_to_cell_vertices(const unsigned int face,
388  const unsigned int vertex,
389  const unsigned char face_orientation) const;
390 
394  unsigned int
395  standard_to_real_face_vertex(const unsigned int vertex,
396  const unsigned int face,
397  const unsigned char face_orientation) const;
398 
402  unsigned int
403  standard_to_real_face_line(const unsigned int line,
404  const unsigned int face,
405  const unsigned char face_orientation) const;
406 
417  bool
418  standard_vs_true_line_orientation(const unsigned int line,
419  const unsigned char face_orientation,
420  const bool line_orientation) const;
421 
432  /*
433  * Return @f$i@f$-th unit tangential vector of a face of the reference cell.
434  * The vectors are arranged such that the
435  * cross product between the two vectors returns the unit normal vector.
436  *
437  * @pre @f$i@f$ must be between zero and `dim-1`.
438  */
439  template <int dim>
441  unit_tangential_vectors(const unsigned int face_no,
442  const unsigned int i) const;
443 
447  template <int dim>
449  unit_normal_vectors(const unsigned int face_no) const;
450 
455  template <typename T, std::size_t N>
456  unsigned char
457  compute_orientation(const std::array<T, N> &vertices_0,
458  const std::array<T, N> &vertices_1) const;
459 
463  template <typename T, std::size_t N>
464  std::array<T, N>
465  permute_according_orientation(const std::array<T, N> &vertices,
466  const unsigned int orientation) const;
467 
472  faces_for_given_vertex(const unsigned int vertex_index) const;
473 
486  unsigned int
487  exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const;
488 
492  unsigned int
493  exodusii_face_to_deal_face(const unsigned int face_n) const;
494 
498  unsigned int
499  unv_vertex_to_deal_vertex(const unsigned int vertex_n) const;
500 
504  unsigned int
505  vtk_linear_type() const;
506 
511  unsigned int
512  vtk_quadratic_type() const;
513 
518  unsigned int
519  vtk_lagrange_type() const;
520 
524  unsigned int
525  gmsh_element_type() const;
526 
540  std::string
541  to_string() const;
542 
546  constexpr operator std::uint8_t() const;
547 
551  constexpr bool
552  operator==(const ReferenceCell &type) const;
553 
557  constexpr bool
558  operator!=(const ReferenceCell &type) const;
559 
565  template <class Archive>
566  void
567  serialize(Archive &archive, const unsigned int /*version*/);
568 
573 private:
577  std::uint8_t kind;
578 
585  constexpr ReferenceCell(const std::uint8_t kind);
586 
591  friend DEAL_II_CONSTEXPR ReferenceCell
593 
594  friend std::ostream &
595  operator<<(std::ostream &out, const ReferenceCell &reference_cell);
596 
597  friend std::istream &
598  operator>>(std::istream &in, ReferenceCell &reference_cell);
599 };
600 
601 
609 std::ostream &
610 operator<<(std::ostream &out, const ReferenceCell &reference_cell);
611 
619 std::istream &
620 operator>>(std::istream &in, ReferenceCell &reference_cell);
621 
622 
623 inline constexpr ReferenceCell::ReferenceCell(const std::uint8_t kind)
624  : kind(kind)
625 {}
626 
627 
628 
629 inline constexpr ReferenceCell::operator std::uint8_t() const
630 {
631  return kind;
632 }
633 
634 
635 
636 inline constexpr bool
638 {
639  return kind == type.kind;
640 }
641 
642 
643 
644 inline constexpr bool
646 {
647  return kind != type.kind;
648 }
649 
650 
651 
652 namespace internal
653 {
654  namespace ReferenceCell
655  {
656  inline DEAL_II_CONSTEXPR ::ReferenceCell
657  make_reference_cell_from_int(const std::uint8_t kind)
658  {
659  // Make sure these are the only indices from which objects can be
660  // created.
661  Assert((kind == static_cast<std::uint8_t>(-1)) || (kind < 8),
662  ExcInternalError());
663 
664  // Call the private constructor, which we can from here because this
665  // function is a 'friend'.
666  return {kind};
667  }
668  } // namespace ReferenceCell
669 } // namespace internal
670 
671 
672 
680 namespace ReferenceCells
681 {
700  static_cast<std::uint8_t>(-1));
701 
707  template <int dim>
708  constexpr const ReferenceCell &
709  get_simplex();
710 
716  template <int dim>
717  constexpr const ReferenceCell &
718  get_hypercube();
719 } // namespace ReferenceCells
720 
721 
722 
723 inline DEAL_II_CONSTEXPR
725  : ReferenceCell(ReferenceCells::Invalid)
726 {}
727 
728 
729 
730 template <class Archive>
731 inline void
732 ReferenceCell::serialize(Archive &archive, const unsigned int /*version*/)
733 {
734  archive &kind;
735 }
736 
737 
738 
741 {
742  if (*this == ReferenceCells::Line)
743  {
745  return {&GeometryInfo<2>::vertex_to_face[vertex][0], 1};
746  }
747  else if (*this == ReferenceCells::Quadrilateral)
748  {
750  return {&GeometryInfo<2>::vertex_to_face[vertex][0], 2};
751  }
752  else if (*this == ReferenceCells::Hexahedron)
753  {
755  return {&GeometryInfo<3>::vertex_to_face[vertex][0], 3};
756  }
757  else if (*this == ReferenceCells::Triangle)
758  {
759  AssertIndexRange(vertex, 3);
760  static const ndarray<unsigned int, 3, 2> table = {
761  {{{0, 2}}, {{0, 1}}, {{1, 2}}}};
762 
763  return table[vertex];
764  }
765  else if (*this == ReferenceCells::Tetrahedron)
766  {
767  AssertIndexRange(vertex, 4);
768  static const ndarray<unsigned int, 4, 3> table = {
769  {{{0, 1, 2}}, {{0, 1, 3}}, {{0, 2, 3}}, {{1, 2, 3}}}};
770 
771  return table[vertex];
772  }
773  else if (*this == ReferenceCells::Wedge)
774  {
775  AssertIndexRange(vertex, 6);
776  static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 4}},
777  {{0, 2, 3}},
778  {{0, 3, 4}},
779  {{1, 2, 4}},
780  {{1, 2, 3}},
781  {{1, 3, 4}}}};
782 
783  return table[vertex];
784  }
785  else if (*this == ReferenceCells::Pyramid)
786  {
787  AssertIndexRange(vertex, 5);
788  static const unsigned int X = numbers::invalid_unsigned_int;
789  static const ndarray<unsigned int, 5, 4> table = {{{{0, 1, 3, X}},
790  {{0, 2, 3, X}},
791  {{0, 1, 4, X}},
792  {{0, 2, 4, X}},
793  {{1, 2, 3, 4}}}};
794 
795  return {&table[vertex][0], vertex == 4 ? 4u : 3u};
796  }
797 
798  Assert(false, ExcNotImplemented());
799 
800  return {};
801 }
802 
803 
804 
805 inline bool
807 {
808  return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
810  *this == ReferenceCells::Hexahedron);
811 }
812 
813 
814 
815 inline bool
817 {
818  return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
819  *this == ReferenceCells::Triangle ||
820  *this == ReferenceCells::Tetrahedron);
821 }
822 
823 
824 
825 inline unsigned int
827 {
828  if (*this == ReferenceCells::Vertex)
829  return 0;
830  else if (*this == ReferenceCells::Line)
831  return 1;
832  else if ((*this == ReferenceCells::Triangle) ||
834  return 2;
835  else if ((*this == ReferenceCells::Tetrahedron) ||
836  (*this == ReferenceCells::Pyramid) ||
837  (*this == ReferenceCells::Wedge) ||
838  (*this == ReferenceCells::Hexahedron))
839  return 3;
840 
841  Assert(false, ExcNotImplemented());
843 }
844 
845 
846 
847 inline unsigned int
849 {
850  if (*this == ReferenceCells::Vertex)
851  return 1;
852  else if (*this == ReferenceCells::Line)
853  return 2;
854  else if (*this == ReferenceCells::Triangle)
855  return 3;
856  else if (*this == ReferenceCells::Quadrilateral)
857  return 4;
858  else if (*this == ReferenceCells::Tetrahedron)
859  return 4;
860  else if (*this == ReferenceCells::Pyramid)
861  return 5;
862  else if (*this == ReferenceCells::Wedge)
863  return 6;
864  else if (*this == ReferenceCells::Hexahedron)
865  return 8;
866 
867  Assert(false, ExcNotImplemented());
869 }
870 
871 
872 
873 inline unsigned int
875 {
876  if (*this == ReferenceCells::Vertex)
877  return 0;
878  else if (*this == ReferenceCells::Line)
879  return 1;
880  else if (*this == ReferenceCells::Triangle)
881  return 3;
882  else if (*this == ReferenceCells::Quadrilateral)
883  return 4;
884  else if (*this == ReferenceCells::Tetrahedron)
885  return 6;
886  else if (*this == ReferenceCells::Pyramid)
887  return 7;
888  else if (*this == ReferenceCells::Wedge)
889  return 9;
890  else if (*this == ReferenceCells::Hexahedron)
891  return 12;
892 
893  Assert(false, ExcNotImplemented());
895 }
896 
897 
898 
899 template <int dim>
901 ReferenceCell::vertex(const unsigned int v) const
902 {
905 
906  if ((dim == 0) && (*this == ReferenceCells::Vertex))
907  {
908  return Point<dim>(0);
909  }
910  else if ((dim == 1) && (*this == ReferenceCells::Line))
911  {
912  static const Point<dim> vertices[2] = {
913  Point<dim>(), // the origin
914  Point<dim>::unit_vector(0) // unit point along x-axis
915  };
916  return vertices[v];
917  }
918  else if ((dim == 2) && (*this == ReferenceCells::Quadrilateral))
919  {
920  static const Point<dim> vertices[4] = {
921  // First the two points on the x-axis
922  Point<dim>(),
924  // Then these two points shifted in the y-direction
927  return vertices[v];
928  }
929  else if ((dim == 3) && (*this == ReferenceCells::Hexahedron))
930  {
931  static const Point<dim> vertices[8] = {
932  // First the two points on the x-axis
933  Point<dim>(),
935  // Then these two points shifted in the y-direction
938  // And now all four points shifted in the z-direction
944  return vertices[v];
945  }
946  else if ((dim == 2) && (*this == ReferenceCells::Triangle))
947  {
948  static const Point<dim> vertices[3] = {
949  Point<dim>(), // the origin
950  Point<dim>::unit_vector(0), // unit point along x-axis
951  Point<dim>::unit_vector(1) // unit point along y-axis
952  };
953  return vertices[v];
954  }
955  else if ((dim == 3) && (*this == ReferenceCells::Tetrahedron))
956  {
957  static const Point<dim> vertices[4] = {
958  Point<dim>(), // the origin
959  Point<dim>::unit_vector(0), // unit point along x-axis
960  Point<dim>::unit_vector(1), // unit point along y-axis
961  Point<dim>::unit_vector(2) // unit point along z-axis
962  };
963  return vertices[v];
964  }
965  else if ((dim == 3) && (*this == ReferenceCells::Pyramid))
966  {
967  static const Point<dim> vertices[5] = {Point<dim>{-1.0, -1.0, 0.0},
968  Point<dim>{+1.0, -1.0, 0.0},
969  Point<dim>{-1.0, +1.0, 0.0},
970  Point<dim>{+1.0, +1.0, 0.0},
971  Point<dim>{+0.0, +0.0, 1.0}};
972  return vertices[v];
973  }
974  else if ((dim == 3) && (*this == ReferenceCells::Wedge))
975  {
976  static const Point<dim> vertices[6] = {
977  // First the three points on the triangular base of the wedge:
978  Point<dim>(),
981  // And now everything shifted in the z-direction again
985  return vertices[v];
986  }
987  else
988  {
989  Assert(false, ExcNotImplemented());
990  return Point<dim>();
991  }
992 }
993 
994 
995 inline unsigned int
997 {
998  if (*this == ReferenceCells::Vertex)
999  return 0;
1000  else if (*this == ReferenceCells::Line)
1001  return 2;
1002  else if (*this == ReferenceCells::Triangle)
1003  return 3;
1004  else if (*this == ReferenceCells::Quadrilateral)
1005  return 4;
1006  else if (*this == ReferenceCells::Tetrahedron)
1007  return 4;
1008  else if (*this == ReferenceCells::Pyramid)
1009  return 5;
1010  else if (*this == ReferenceCells::Wedge)
1011  return 5;
1012  else if (*this == ReferenceCells::Hexahedron)
1013  return 6;
1014 
1015  Assert(false, ExcNotImplemented());
1017 }
1018 
1019 
1020 
1023 {
1024  return {0U, n_vertices()};
1025 }
1026 
1027 
1028 
1031 {
1032  return {0U, n_lines()};
1033 }
1034 
1035 
1036 
1039 {
1040  return {0U, n_faces()};
1041 }
1042 
1043 
1044 
1045 inline ReferenceCell
1046 ReferenceCell::face_reference_cell(const unsigned int face_no) const
1047 {
1048  AssertIndexRange(face_no, n_faces());
1049 
1050  if (*this == ReferenceCells::Vertex)
1051  return ReferenceCells::Invalid;
1052  else if (*this == ReferenceCells::Line)
1053  return ReferenceCells::Vertex;
1054  else if (*this == ReferenceCells::Triangle)
1055  return ReferenceCells::Line;
1056  else if (*this == ReferenceCells::Quadrilateral)
1057  return ReferenceCells::Line;
1058  else if (*this == ReferenceCells::Tetrahedron)
1059  return ReferenceCells::Triangle;
1060  else if (*this == ReferenceCells::Pyramid)
1061  {
1062  if (face_no == 0)
1064  else
1065  return ReferenceCells::Triangle;
1066  }
1067  else if (*this == ReferenceCells::Wedge)
1068  {
1069  if (face_no > 1)
1071  else
1072  return ReferenceCells::Triangle;
1073  }
1074  else if (*this == ReferenceCells::Hexahedron)
1076 
1077  Assert(false, ExcNotImplemented());
1078  return ReferenceCells::Invalid;
1079 }
1080 
1081 
1082 
1083 inline unsigned int
1085  const unsigned int face,
1086  const unsigned int subface,
1087  const unsigned char face_orientation_raw) const
1088 {
1089  AssertIndexRange(face, n_faces());
1090 
1091  if (*this == ReferenceCells::Vertex)
1092  {
1093  Assert(false, ExcNotImplemented());
1094  }
1095  else if (*this == ReferenceCells::Line)
1096  {
1097  Assert(false, ExcNotImplemented());
1098  }
1099  else if (*this == ReferenceCells::Triangle)
1100  {
1101  static const ndarray<unsigned int, 3, 2> subcells = {
1102  {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
1103 
1104  return subcells[face][subface];
1105  }
1106  else if (*this == ReferenceCells::Quadrilateral)
1107  {
1108  const bool face_orientation = Utilities::get_bit(face_orientation_raw, 0);
1109  const bool face_flip = Utilities::get_bit(face_orientation_raw, 2);
1110  const bool face_rotation = Utilities::get_bit(face_orientation_raw, 1);
1111 
1114  face,
1115  subface,
1116  face_orientation,
1117  face_flip,
1118  face_rotation);
1119  }
1120  else if (*this == ReferenceCells::Tetrahedron)
1121  {
1122  Assert(false, ExcNotImplemented());
1123  }
1124  else if (*this == ReferenceCells::Pyramid)
1125  {
1126  Assert(false, ExcNotImplemented());
1127  }
1128  else if (*this == ReferenceCells::Wedge)
1129  {
1130  Assert(false, ExcNotImplemented());
1131  }
1132  else if (*this == ReferenceCells::Hexahedron)
1133  {
1134  const bool face_orientation = Utilities::get_bit(face_orientation_raw, 0);
1135  const bool face_flip = Utilities::get_bit(face_orientation_raw, 2);
1136  const bool face_rotation = Utilities::get_bit(face_orientation_raw, 1);
1137 
1140  face,
1141  subface,
1142  face_orientation,
1143  face_flip,
1144  face_rotation);
1145  }
1146 
1147  Assert(false, ExcNotImplemented());
1148  return {};
1149 }
1150 
1151 
1152 
1153 inline std::array<unsigned int, 2>
1155  const unsigned int vertex) const
1156 {
1157  AssertIndexRange(vertex, n_vertices());
1158 
1159  if (*this == ReferenceCells::Vertex)
1160  {
1161  Assert(false, ExcNotImplemented());
1162  }
1163  else if (*this == ReferenceCells::Line)
1164  {
1165  Assert(false, ExcNotImplemented());
1166  }
1167  else if (*this == ReferenceCells::Triangle)
1168  {
1169  static const ndarray<unsigned int, 3, 2> table = {
1170  {{{0, 0}}, {{0, 1}}, {{1, 1}}}};
1171 
1172  return table[vertex];
1173  }
1174  else if (*this == ReferenceCells::Quadrilateral)
1175  {
1177  }
1178  else if (*this == ReferenceCells::Tetrahedron)
1179  {
1180  static const ndarray<unsigned int, 4, 2> table = {
1181  {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 2}}}};
1182 
1183  return table[vertex];
1184  }
1185  else if (*this == ReferenceCells::Pyramid)
1186  {
1187  static const ndarray<unsigned int, 5, 2> table = {
1188  {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{0, 3}}, {{1, 2}}}};
1189 
1190  return table[vertex];
1191  }
1192  else if (*this == ReferenceCells::Wedge)
1193  {
1194  static const ndarray<unsigned int, 6, 2> table = {
1195  {{{0, 1}}, {{0, 0}}, {{0, 2}}, {{1, 0}}, {{1, 1}}, {{1, 2}}}};
1196 
1197  return table[vertex];
1198  }
1199  else if (*this == ReferenceCells::Hexahedron)
1200  {
1202  }
1203 
1204  Assert(false, ExcNotImplemented());
1205  return {};
1206 }
1207 
1208 
1209 
1210 inline std::array<unsigned int, 2>
1212  const unsigned int line) const
1213 {
1214  AssertIndexRange(line, n_lines());
1215 
1216  if (*this == ReferenceCells::Vertex)
1217  {
1218  Assert(false, ExcNotImplemented());
1219  }
1220  else if (*this == ReferenceCells::Line)
1221  {
1222  Assert(false, ExcNotImplemented());
1223  }
1224  else if (*this == ReferenceCells::Triangle)
1225  {
1226  Assert(false, ExcNotImplemented());
1227  }
1228  else if (*this == ReferenceCells::Quadrilateral)
1229  {
1230  Assert(false, ExcNotImplemented());
1231  }
1232  else if (*this == ReferenceCells::Tetrahedron)
1233  {
1234  static const std::array<unsigned int, 2> table[6] = {
1235  {{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}};
1236 
1237  return table[line];
1238  }
1239  else if (*this == ReferenceCells::Pyramid)
1240  {
1241  static const std::array<unsigned int, 2> table[8] = {{{0, 0}},
1242  {{0, 1}},
1243  {{0, 2}},
1244  {{0, 3}},
1245  {{1, 2}},
1246  {{2, 1}},
1247  {{1, 1}},
1248  {{2, 2}}};
1249 
1250  return table[line];
1251  }
1252  else if (*this == ReferenceCells::Wedge)
1253  {
1254  static const std::array<unsigned int, 2> table[9] = {{{0, 0}},
1255  {{0, 2}},
1256  {{0, 1}},
1257  {{1, 0}},
1258  {{1, 1}},
1259  {{1, 2}},
1260  {{2, 0}},
1261  {{2, 1}},
1262  {{3, 1}}};
1263 
1264  return table[line];
1265  }
1266  else if (*this == ReferenceCells::Hexahedron)
1267  {
1269  }
1270 
1271  Assert(false, ExcNotImplemented());
1272  return {};
1273 }
1274 
1275 
1276 
1277 inline unsigned int
1278 ReferenceCell::face_to_cell_lines(const unsigned int face,
1279  const unsigned int line,
1280  const unsigned char face_orientation) const
1281 {
1282  AssertIndexRange(face, n_faces());
1284 
1285  if (*this == ReferenceCells::Vertex)
1286  {
1287  Assert(false, ExcNotImplemented());
1288  }
1289  else if (*this == ReferenceCells::Line)
1290  {
1292  face,
1293  line,
1294  Utilities::get_bit(face_orientation, 0),
1295  Utilities::get_bit(face_orientation, 2),
1296  Utilities::get_bit(face_orientation, 1));
1297  }
1298  else if (*this == ReferenceCells::Triangle)
1299  {
1300  return face;
1301  }
1302  else if (*this == ReferenceCells::Quadrilateral)
1303  {
1305  face,
1306  line,
1307  Utilities::get_bit(face_orientation, 0),
1308  Utilities::get_bit(face_orientation, 2),
1309  Utilities::get_bit(face_orientation, 1));
1310  }
1311  else if (*this == ReferenceCells::Tetrahedron)
1312  {
1313  const static ndarray<unsigned int, 4, 3> table = {
1314  {{{0, 1, 2}}, {{0, 3, 4}}, {{2, 5, 3}}, {{1, 4, 5}}}};
1315 
1316  return table[face]
1317  [standard_to_real_face_line(line, face, face_orientation)];
1318  }
1319  else if (*this == ReferenceCells::Pyramid)
1320  {
1321  Assert(false, ExcNotImplemented());
1322  }
1323  else if (*this == ReferenceCells::Wedge)
1324  {
1325  Assert(false, ExcNotImplemented());
1326  }
1327  else if (*this == ReferenceCells::Hexahedron)
1328  {
1330  face,
1331  line,
1332  Utilities::get_bit(face_orientation, 0),
1333  Utilities::get_bit(face_orientation, 2),
1334  Utilities::get_bit(face_orientation, 1));
1335  }
1336 
1337  Assert(false, ExcNotImplemented());
1339 }
1340 
1341 
1342 
1343 inline unsigned int
1345  const unsigned int vertex,
1346  const unsigned char face_orientation) const
1347 {
1348  AssertIndexRange(face, n_faces());
1350 
1351  if (*this == ReferenceCells::Vertex)
1352  {
1353  Assert(false, ExcNotImplemented());
1354  }
1355  else if (*this == ReferenceCells::Line)
1356  {
1358  face,
1359  vertex,
1360  Utilities::get_bit(face_orientation, 0),
1361  Utilities::get_bit(face_orientation, 2),
1362  Utilities::get_bit(face_orientation, 1));
1363  }
1364  else if (*this == ReferenceCells::Triangle)
1365  {
1366  static const ndarray<unsigned int, 3, 2> table = {
1367  {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
1368 
1369  return table[face][face_orientation != 0u ? vertex : (1 - vertex)];
1370  }
1371  else if (*this == ReferenceCells::Quadrilateral)
1372  {
1374  face,
1375  vertex,
1376  Utilities::get_bit(face_orientation, 0),
1377  Utilities::get_bit(face_orientation, 2),
1378  Utilities::get_bit(face_orientation, 1));
1379  }
1380  else if (*this == ReferenceCells::Tetrahedron)
1381  {
1382  static const ndarray<unsigned int, 4, 3> table = {
1383  {{{0, 1, 2}}, {{1, 0, 3}}, {{0, 2, 3}}, {{2, 1, 3}}}};
1384 
1385  return table[face][standard_to_real_face_vertex(
1386  vertex, face, face_orientation)];
1387  }
1388  else if (*this == ReferenceCells::Pyramid)
1389  {
1390  constexpr auto X = numbers::invalid_unsigned_int;
1391  static const ndarray<unsigned int, 5, 4> table = {{{{0, 1, 2, 3}},
1392  {{0, 2, 4, X}},
1393  {{3, 1, 4, X}},
1394  {{1, 0, 4, X}},
1395  {{2, 3, 4, X}}}};
1396 
1397  return table[face][standard_to_real_face_vertex(
1398  vertex, face, face_orientation)];
1399  }
1400  else if (*this == ReferenceCells::Wedge)
1401  {
1402  constexpr auto X = numbers::invalid_unsigned_int;
1403  static const ndarray<unsigned int, 6, 4> table = {{{{1, 0, 2, X}},
1404  {{3, 4, 5, X}},
1405  {{0, 1, 3, 4}},
1406  {{1, 2, 4, 5}},
1407  {{2, 0, 5, 3}}}};
1408 
1409  return table[face][standard_to_real_face_vertex(
1410  vertex, face, face_orientation)];
1411  }
1412  else if (*this == ReferenceCells::Hexahedron)
1413  {
1415  face,
1416  vertex,
1417  Utilities::get_bit(face_orientation, 0),
1418  Utilities::get_bit(face_orientation, 2),
1419  Utilities::get_bit(face_orientation, 1));
1420  }
1421 
1422  Assert(false, ExcNotImplemented());
1424 }
1425 
1426 
1427 
1428 inline unsigned int
1430  const unsigned int vertex,
1431  const unsigned int face,
1432  const unsigned char face_orientation) const
1433 {
1434  AssertIndexRange(face, n_faces());
1436 
1437  if (*this == ReferenceCells::Vertex)
1438  {
1439  Assert(false, ExcNotImplemented());
1440  }
1441  else if (*this == ReferenceCells::Line)
1442  {
1443  Assert(false, ExcNotImplemented());
1444  }
1445  else if (*this == ReferenceCells::Triangle)
1446  {
1447  static const ndarray<unsigned int, 2, 2> table = {{{{1, 0}}, {{0, 1}}}};
1448 
1449  return table[face_orientation][vertex];
1450  }
1451  else if (*this == ReferenceCells::Quadrilateral)
1452  {
1454  face_orientation !=
1455  0u);
1456  }
1457  else if (*this == ReferenceCells::Tetrahedron)
1458  {
1459  static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 1}},
1460  {{0, 1, 2}},
1461  {{2, 1, 0}},
1462  {{1, 2, 0}},
1463  {{1, 0, 2}},
1464  {{2, 0, 1}}}};
1465 
1466  return table[face_orientation][vertex];
1467  }
1468  else if (*this == ReferenceCells::Pyramid)
1469  {
1470  if (face == 0) // The quadrilateral face
1471  {
1473  vertex,
1474  Utilities::get_bit(face_orientation, 0),
1475  Utilities::get_bit(face_orientation, 2),
1476  Utilities::get_bit(face_orientation, 1));
1477  }
1478  else // One of the triangular faces
1479  {
1480  static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 1}},
1481  {{0, 1, 2}},
1482  {{2, 1, 0}},
1483  {{1, 2, 0}},
1484  {{1, 0, 2}},
1485  {{2, 0, 1}}}};
1486 
1487  return table[face_orientation][vertex];
1488  }
1489  }
1490  else if (*this == ReferenceCells::Wedge)
1491  {
1492  if (face > 1) // One of the quadrilateral faces
1493  {
1495  vertex,
1496  Utilities::get_bit(face_orientation, 0),
1497  Utilities::get_bit(face_orientation, 2),
1498  Utilities::get_bit(face_orientation, 1));
1499  }
1500  else // One of the triangular faces
1501  {
1502  static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 1}},
1503  {{0, 1, 2}},
1504  {{2, 1, 0}},
1505  {{1, 2, 0}},
1506  {{1, 0, 2}},
1507  {{2, 0, 1}}}};
1508 
1509  return table[face_orientation][vertex];
1510  }
1511  }
1512  else if (*this == ReferenceCells::Hexahedron)
1513  {
1515  vertex,
1516  Utilities::get_bit(face_orientation, 0),
1517  Utilities::get_bit(face_orientation, 2),
1518  Utilities::get_bit(face_orientation, 1));
1519  }
1520 
1521  Assert(false, ExcNotImplemented());
1523 }
1524 
1525 
1526 
1527 inline unsigned int
1529  const unsigned int line,
1530  const unsigned int face,
1531  const unsigned char face_orientation) const
1532 {
1533  AssertIndexRange(face, n_faces());
1535 
1536  if (*this == ReferenceCells::Vertex)
1537  {
1538  Assert(false, ExcNotImplemented());
1539  }
1540  else if (*this == ReferenceCells::Line)
1541  {
1542  Assert(false, ExcNotImplemented());
1543  }
1544  else if (*this == ReferenceCells::Triangle)
1545  {
1546  Assert(false, ExcNotImplemented());
1547  }
1548  else if (*this == ReferenceCells::Quadrilateral)
1549  {
1550  Assert(false, ExcNotImplemented());
1551  }
1552  else if (*this == ReferenceCells::Tetrahedron)
1553  {
1554  static const ndarray<unsigned int, 6, 3> table = {{{{2, 1, 0}},
1555  {{0, 1, 2}},
1556  {{1, 0, 2}},
1557  {{1, 2, 0}},
1558  {{0, 2, 1}},
1559  {{2, 0, 1}}}};
1560 
1561  return table[face_orientation][line];
1562  }
1563  else if (*this == ReferenceCells::Pyramid)
1564  {
1565  if (face == 0) // The quadrilateral face
1566  {
1568  line,
1569  Utilities::get_bit(face_orientation, 0),
1570  Utilities::get_bit(face_orientation, 2),
1571  Utilities::get_bit(face_orientation, 1));
1572  }
1573  else // One of the triangular faces
1574  {
1575  static const ndarray<unsigned int, 6, 3> table = {{{{2, 1, 0}},
1576  {{0, 1, 2}},
1577  {{1, 0, 2}},
1578  {{1, 2, 0}},
1579  {{0, 2, 1}},
1580  {{2, 0, 1}}}};
1581 
1582  return table[face_orientation][line];
1583  }
1584  }
1585  else if (*this == ReferenceCells::Wedge)
1586  {
1587  if (face > 1) // One of the quadrilateral faces
1588  {
1590  line,
1591  Utilities::get_bit(face_orientation, 0),
1592  Utilities::get_bit(face_orientation, 2),
1593  Utilities::get_bit(face_orientation, 1));
1594  }
1595  else // One of the triangular faces
1596  {
1597  static const ndarray<unsigned int, 6, 3> table = {{{{2, 1, 0}},
1598  {{0, 1, 2}},
1599  {{1, 0, 2}},
1600  {{1, 2, 0}},
1601  {{0, 2, 1}},
1602  {{2, 0, 1}}}};
1603 
1604  return table[face_orientation][line];
1605  }
1606  }
1607  else if (*this == ReferenceCells::Hexahedron)
1608  {
1610  line,
1611  Utilities::get_bit(face_orientation, 0),
1612  Utilities::get_bit(face_orientation, 2),
1613  Utilities::get_bit(face_orientation, 1));
1614  }
1615 
1616  Assert(false, ExcNotImplemented());
1618 }
1619 
1620 
1621 
1622 namespace ReferenceCells
1623 {
1624  template <int dim>
1625  inline constexpr const ReferenceCell &
1627  {
1628  switch (dim)
1629  {
1630  case 0:
1631  return ReferenceCells::Vertex;
1632  case 1:
1633  return ReferenceCells::Line;
1634  case 2:
1635  return ReferenceCells::Triangle;
1636  case 3:
1638  default:
1639  Assert(false, ExcNotImplemented());
1640  return ReferenceCells::Invalid;
1641  }
1642  }
1643 
1644 
1645 
1646  template <int dim>
1647  inline constexpr const ReferenceCell &
1649  {
1650  switch (dim)
1651  {
1652  case 0:
1653  return ReferenceCells::Vertex;
1654  case 1:
1655  return ReferenceCells::Line;
1656  case 2:
1658  case 3:
1660  default:
1661  Assert(false, ExcNotImplemented());
1662  return ReferenceCells::Invalid;
1663  }
1664  }
1665 } // namespace ReferenceCells
1666 
1667 
1668 inline ReferenceCell
1669 ReferenceCell::n_vertices_to_type(const int dim, const unsigned int n_vertices)
1670 {
1671  AssertIndexRange(dim, 4);
1672  AssertIndexRange(n_vertices, 9);
1673 
1674  const auto X = ReferenceCells::Invalid;
1675  static const ndarray<ReferenceCell, 4, 9> table = {
1676  {// dim 0
1677  {{X, ReferenceCells::Vertex, X, X, X, X, X, X, X}},
1678  // dim 1
1679  {{X, X, ReferenceCells::Line, X, X, X, X, X, X}},
1680  // dim 2
1681  {{X,
1682  X,
1683  X,
1686  X,
1687  X,
1688  X,
1689  X}},
1690  // dim 3
1691  {{X,
1692  X,
1693  X,
1694  X,
1698  X,
1700  Assert(table[dim][n_vertices] != ReferenceCells::Invalid,
1701  ExcMessage("The combination of dim = " + std::to_string(dim) +
1702  " and n_vertices = " + std::to_string(n_vertices) +
1703  " does not correspond to a known reference cell type."));
1704  return table[dim][n_vertices];
1705 }
1706 
1707 
1708 
1709 template <int dim>
1710 inline double
1712  const unsigned int i) const
1713 {
1715  if (*this == ReferenceCells::get_hypercube<dim>())
1717 
1718  if (*this ==
1719  ReferenceCells::Triangle) // see also
1720  // BarycentricPolynomials<2>::compute_value
1721  {
1722  switch (i)
1723  {
1724  case 0:
1725  return 1.0 - xi[std::min(0, dim - 1)] - xi[std::min(1, dim - 1)];
1726  case 1:
1727  return xi[std::min(0, dim - 1)];
1728  case 2:
1729  return xi[std::min(1, dim - 1)];
1730  }
1731  }
1732 
1733  if (*this ==
1734  ReferenceCells::Tetrahedron) // see also
1735  // BarycentricPolynomials<3>::compute_value
1736  {
1737  switch (i)
1738  {
1739  case 0:
1740  return 1.0 - xi[std::min(0, dim - 1)] - xi[std::min(1, dim - 1)] -
1741  xi[std::min(2, dim - 1)];
1742  case 1:
1743  return xi[std::min(0, dim - 1)];
1744  case 2:
1745  return xi[std::min(1, dim - 1)];
1746  case 3:
1747  return xi[std::min(2, dim - 1)];
1748  }
1749  }
1750 
1751  if (*this ==
1752  ReferenceCells::Wedge) // see also
1753  // ScalarLagrangePolynomialWedge::compute_value
1754  {
1756  .d_linear_shape_function<2>(Point<2>(xi[std::min(0, dim - 1)],
1757  xi[std::min(1, dim - 1)]),
1758  i % 3) *
1760  .d_linear_shape_function<1>(Point<1>(xi[std::min(2, dim - 1)]),
1761  i / 3);
1762  }
1763 
1764  if (*this ==
1765  ReferenceCells::Pyramid) // see also
1766  // ScalarLagrangePolynomialPyramid::compute_value
1767  {
1768  const double Q14 = 0.25;
1769  double ration;
1770 
1771  const double r = xi[std::min(0, dim - 1)];
1772  const double s = xi[std::min(1, dim - 1)];
1773  const double t = xi[std::min(2, dim - 1)];
1774 
1775  if (fabs(t - 1.0) > 1.0e-14)
1776  {
1777  ration = (r * s * t) / (1.0 - t);
1778  }
1779  else
1780  {
1781  ration = 0.0;
1782  }
1783 
1784  if (i == 0)
1785  return Q14 * ((1.0 - r) * (1.0 - s) - t + ration);
1786  if (i == 1)
1787  return Q14 * ((1.0 + r) * (1.0 - s) - t - ration);
1788  if (i == 2)
1789  return Q14 * ((1.0 - r) * (1.0 + s) - t - ration);
1790  if (i == 3)
1791  return Q14 * ((1.0 + r) * (1.0 + s) - t + ration);
1792  else
1793  return t;
1794  }
1795 
1796  Assert(false, ExcNotImplemented());
1797 
1798  return 0.0;
1799 }
1800 
1801 
1802 
1803 template <int dim>
1804 inline Tensor<1, dim>
1806  const unsigned int i) const
1807 {
1809  if (*this == ReferenceCells::get_hypercube<dim>())
1811 
1812  if (*this ==
1813  ReferenceCells::Triangle) // see also
1814  // BarycentricPolynomials<2>::compute_grad
1815  {
1816  switch (i)
1817  {
1818  case 0:
1819  return Point<dim>(-1.0, -1.0);
1820  case 1:
1821  return Point<dim>(+1.0, +0.0);
1822  case 2:
1823  return Point<dim>(+0.0, +1.0);
1824  }
1825  }
1826 
1827  Assert(false, ExcNotImplemented());
1828 
1829  return Point<dim>(+0.0, +0.0, +0.0);
1830 }
1831 
1832 
1833 template <int dim>
1834 inline Tensor<1, dim>
1835 ReferenceCell::unit_tangential_vectors(const unsigned int face_no,
1836  const unsigned int i) const
1837 {
1839  AssertIndexRange(i, dim - 1);
1840 
1841  if (*this == ReferenceCells::get_hypercube<dim>())
1842  {
1844  return GeometryInfo<dim>::unit_tangential_vectors[face_no][i];
1845  }
1846  else if (*this == ReferenceCells::Triangle)
1847  {
1848  AssertIndexRange(face_no, 3);
1849  static const std::array<Tensor<1, dim>, 3> table = {
1850  {Point<dim>(1, 0),
1851  Point<dim>(-std::sqrt(0.5), +std::sqrt(0.5)),
1852  Point<dim>(0, -1)}};
1853 
1854  return table[face_no];
1855  }
1856  else if (*this == ReferenceCells::Tetrahedron)
1857  {
1858  AssertIndexRange(face_no, 4);
1859  static const ndarray<Tensor<1, dim>, 4, 2> table = {
1860  {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
1861  {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
1862  {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
1863  {{Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
1864  +std::pow(1.0 / 3.0, 1.0 / 4.0),
1865  0),
1866  Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
1867  0,
1868  +std::pow(1.0 / 3.0, 1.0 / 4.0))}}}};
1869 
1870  return table[face_no][i];
1871  }
1872  else if (*this == ReferenceCells::Wedge)
1873  {
1874  AssertIndexRange(face_no, 5);
1875  static const ndarray<Tensor<1, dim>, 5, 2> table = {
1876  {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
1877  {{Point<dim>(1, 0, 0), Point<dim>(0, 1, 0)}},
1878  {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
1879  {{Point<dim>(-1 / std::sqrt(2.0), +1 / std::sqrt(2.0), 0),
1880  Point<dim>(0, 0, 1)}},
1881  {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}}}};
1882 
1883  return table[face_no][i];
1884  }
1885  else if (*this == ReferenceCells::Pyramid)
1886  {
1887  AssertIndexRange(face_no, 5);
1888  static const ndarray<Tensor<1, dim>, 5, 2> table = {
1889  {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
1890  {{Point<dim>(+1.0 / sqrt(2.0), 0, +1.0 / sqrt(2.0)),
1891  Point<dim>(0, 1, 0)}},
1892  {{Point<dim>(+1.0 / sqrt(2.0), 0, -1.0 / sqrt(2.0)),
1893  Point<dim>(0, 1, 0)}},
1894  {{Point<dim>(1, 0, 0),
1895  Point<dim>(0, +1.0 / sqrt(2.0), +1.0 / sqrt(2.0))}},
1896  {{Point<dim>(1, 0, 0),
1897  Point<dim>(0, +1.0 / sqrt(2.0), -1.0 / sqrt(2.0))}}}};
1898 
1899  return table[face_no][i];
1900  }
1901 
1902  Assert(false, ExcNotImplemented());
1903 
1904  return {};
1905 }
1906 
1907 
1908 
1909 template <int dim>
1910 inline Tensor<1, dim>
1911 ReferenceCell::unit_normal_vectors(const unsigned int face_no) const
1912 {
1913  AssertDimension(dim, this->get_dimension());
1914 
1915  if (is_hyper_cube())
1916  {
1918  return GeometryInfo<dim>::unit_normal_vector[face_no];
1919  }
1920  else if (dim == 2)
1921  {
1923 
1924  // Return the rotated vector
1925  return cross_product_2d(unit_tangential_vectors<dim>(face_no, 0));
1926  }
1927  else if (dim == 3)
1928  {
1929  return cross_product_3d(unit_tangential_vectors<dim>(face_no, 0),
1930  unit_tangential_vectors<dim>(face_no, 1));
1931  }
1932 
1933  Assert(false, ExcNotImplemented());
1934 
1935  return {};
1936 }
1937 
1938 
1939 
1940 inline bool
1942  const unsigned int line,
1943  const unsigned char face_orientation_raw,
1944  const bool line_orientation) const
1945 {
1946  if (*this == ReferenceCells::Hexahedron)
1947  {
1948  static const bool bool_table[2][2][2][2] = {
1949  {{{true, false}, // lines 0/1, face_orientation=false,
1950  // face_flip=false, face_rotation=false and true
1951  {false, true}}, // lines 0/1, face_orientation=false,
1952  // face_flip=true, face_rotation=false and true
1953  {{true, true}, // lines 0/1, face_orientation=true,
1954  // face_flip=false, face_rotation=false and true
1955  {false, false}}}, // lines 0/1, face_orientation=true,
1956  // face_flip=true, face_rotation=false and true
1957 
1958  {{{true, true}, // lines 2/3 ...
1959  {false, false}},
1960  {{true, false}, {false, true}}}};
1961 
1962  const bool face_orientation = Utilities::get_bit(face_orientation_raw, 0);
1963  const bool face_flip = Utilities::get_bit(face_orientation_raw, 2);
1964  const bool face_rotation = Utilities::get_bit(face_orientation_raw, 1);
1965 
1966  return (line_orientation ==
1967  bool_table[line / 2][face_orientation][face_flip][face_rotation]);
1968  }
1969  else
1970  // TODO: This might actually be wrong for some of the other
1971  // kinds of objects. We should check this
1972  return true;
1973 }
1974 
1975 
1976 
1977 namespace internal
1978 {
1979  template <typename T, std::size_t N>
1981  {
1982  public:
1986  NoPermutation(const ::ReferenceCell &entity_type,
1987  const std::array<T, N> & vertices_0,
1988  const std::array<T, N> & vertices_1)
1989  : entity_type(entity_type)
1990  , vertices_0(vertices_0)
1991  , vertices_1(vertices_1)
1992  {}
1993 
1997  virtual ~NoPermutation() noexcept override = default;
1998 
2002  virtual void
2003  print_info(std::ostream &out) const override
2004  {
2005  out << "[";
2006 
2007  const unsigned int n_vertices = entity_type.n_vertices();
2008 
2009  for (unsigned int i = 0; i < n_vertices; ++i)
2010  {
2011  out << vertices_0[i];
2012  if (i + 1 != n_vertices)
2013  out << ",";
2014  }
2015 
2016  out << "] is not a permutation of [";
2017 
2018  for (unsigned int i = 0; i < n_vertices; ++i)
2019  {
2020  out << vertices_1[i];
2021  if (i + 1 != n_vertices)
2022  out << ",";
2023  }
2024 
2025  out << "]." << std::endl;
2026  }
2027 
2031  const ::ReferenceCell entity_type;
2032 
2036  const std::array<T, N> vertices_0;
2037 
2041  const std::array<T, N> vertices_1;
2042  };
2043 } // namespace internal
2044 
2045 
2046 
2047 template <typename T, std::size_t N>
2048 inline unsigned char
2049 ReferenceCell::compute_orientation(const std::array<T, N> &vertices_0,
2050  const std::array<T, N> &vertices_1) const
2051 {
2052  AssertIndexRange(n_vertices(), N + 1);
2053  if (*this == ReferenceCells::Line)
2054  {
2055  const std::array<T, 2> i{{vertices_0[0], vertices_0[1]}};
2056  const std::array<T, 2> j{{vertices_1[0], vertices_1[1]}};
2057 
2058  // line_orientation=true
2059  if (i == std::array<T, 2>{{j[0], j[1]}})
2060  return 1;
2061 
2062  // line_orientation=false
2063  if (i == std::array<T, 2>{{j[1], j[0]}})
2064  return 0;
2065  }
2066  else if (*this == ReferenceCells::Triangle)
2067  {
2068  const std::array<T, 3> i{{vertices_0[0], vertices_0[1], vertices_0[2]}};
2069  const std::array<T, 3> j{{vertices_1[0], vertices_1[1], vertices_1[2]}};
2070 
2071  // face_orientation=true, face_rotation=false, face_flip=false
2072  if (i == std::array<T, 3>{{j[0], j[1], j[2]}})
2073  return 1;
2074 
2075  // face_orientation=true, face_rotation=true, face_flip=false
2076  if (i == std::array<T, 3>{{j[1], j[2], j[0]}})
2077  return 3;
2078 
2079  // face_orientation=true, face_rotation=false, face_flip=true
2080  if (i == std::array<T, 3>{{j[2], j[0], j[1]}})
2081  return 5;
2082 
2083  // face_orientation=false, face_rotation=false, face_flip=false
2084  if (i == std::array<T, 3>{{j[0], j[2], j[1]}})
2085  return 0;
2086 
2087  // face_orientation=false, face_rotation=true, face_flip=false
2088  if (i == std::array<T, 3>{{j[2], j[1], j[0]}})
2089  return 2;
2090 
2091  // face_orientation=false, face_rotation=false, face_flip=true
2092  if (i == std::array<T, 3>{{j[1], j[0], j[2]}})
2093  return 4;
2094  }
2095  else if (*this == ReferenceCells::Quadrilateral)
2096  {
2097  const std::array<T, 4> i{
2098  {vertices_0[0], vertices_0[1], vertices_0[2], vertices_0[3]}};
2099  const std::array<T, 4> j{
2100  {vertices_1[0], vertices_1[1], vertices_1[2], vertices_1[3]}};
2101 
2102  // face_orientation=true, face_rotation=false, face_flip=false
2103  if (i == std::array<T, 4>{{j[0], j[1], j[2], j[3]}})
2104  return 1;
2105 
2106  // face_orientation=true, face_rotation=true, face_flip=false
2107  if (i == std::array<T, 4>{{j[2], j[0], j[3], j[1]}})
2108  return 3;
2109 
2110  // face_orientation=true, face_rotation=false, face_flip=true
2111  if (i == std::array<T, 4>{{j[3], j[2], j[1], j[0]}})
2112  return 5;
2113 
2114  // face_orientation=true, face_rotation=true, face_flip=true
2115  if (i == std::array<T, 4>{{j[1], j[3], j[0], j[2]}})
2116  return 7;
2117 
2118  // face_orientation=false, face_rotation=false, face_flip=false
2119  if (i == std::array<T, 4>{{j[0], j[2], j[1], j[3]}})
2120  return 0;
2121 
2122  // face_orientation=false, face_rotation=true, face_flip=false
2123  if (i == std::array<T, 4>{{j[2], j[3], j[0], j[1]}})
2124  return 2;
2125 
2126  // face_orientation=false, face_rotation=false, face_flip=true
2127  if (i == std::array<T, 4>{{j[3], j[1], j[2], j[0]}})
2128  return 4;
2129 
2130  // face_orientation=false, face_rotation=true, face_flip=true
2131  if (i == std::array<T, 4>{{j[1], j[0], j[3], j[2]}})
2132  return 6;
2133  }
2134 
2135  Assert(false, (internal::NoPermutation<T, N>(*this, vertices_0, vertices_1)));
2136 
2137  return -1;
2138 }
2139 
2140 
2141 
2142 template <typename T, std::size_t N>
2143 inline std::array<T, N>
2145  const std::array<T, N> &vertices,
2146  const unsigned int orientation) const
2147 {
2148  std::array<T, 4> temp;
2149 
2150  if (*this == ReferenceCells::Line)
2151  {
2152  switch (orientation)
2153  {
2154  case 1:
2155  temp = {{vertices[0], vertices[1]}};
2156  break;
2157  case 0:
2158  temp = {{vertices[1], vertices[0]}};
2159  break;
2160  default:
2161  Assert(false, ExcNotImplemented());
2162  }
2163  }
2164  else if (*this == ReferenceCells::Triangle)
2165  {
2166  switch (orientation)
2167  {
2168  case 1:
2169  temp = {{vertices[0], vertices[1], vertices[2]}};
2170  break;
2171  case 3:
2172  temp = {{vertices[1], vertices[2], vertices[0]}};
2173  break;
2174  case 5:
2175  temp = {{vertices[2], vertices[0], vertices[1]}};
2176  break;
2177  case 0:
2178  temp = {{vertices[0], vertices[2], vertices[1]}};
2179  break;
2180  case 2:
2181  temp = {{vertices[2], vertices[1], vertices[0]}};
2182  break;
2183  case 4:
2184  temp = {{vertices[1], vertices[0], vertices[2]}};
2185  break;
2186  default:
2187  Assert(false, ExcNotImplemented());
2188  }
2189  }
2190  else if (*this == ReferenceCells::Quadrilateral)
2191  {
2192  switch (orientation)
2193  {
2194  case 1:
2195  temp = {{vertices[0], vertices[1], vertices[2], vertices[3]}};
2196  break;
2197  case 3:
2198  temp = {{vertices[2], vertices[0], vertices[3], vertices[1]}};
2199  break;
2200  case 5:
2201  temp = {{vertices[3], vertices[2], vertices[1], vertices[0]}};
2202  break;
2203  case 7:
2204  temp = {{vertices[1], vertices[3], vertices[0], vertices[2]}};
2205  break;
2206  case 0:
2207  temp = {{vertices[0], vertices[2], vertices[1], vertices[3]}};
2208  break;
2209  case 2:
2210  temp = {{vertices[2], vertices[3], vertices[0], vertices[1]}};
2211  break;
2212  case 4:
2213  temp = {{vertices[3], vertices[1], vertices[2], vertices[0]}};
2214  break;
2215  case 6:
2216  temp = {{vertices[1], vertices[0], vertices[3], vertices[2]}};
2217  break;
2218  default:
2219  Assert(false, ExcNotImplemented());
2220  }
2221  }
2222  else
2223  {
2224  AssertThrow(false, ExcNotImplemented());
2225  }
2226 
2227  std::array<T, N> temp_;
2228  std::copy_n(temp.begin(), N, temp_.begin());
2229 
2230  return temp_;
2231 }
2232 
2233 
2235 
2236 #endif
unsigned int standard_to_real_face_line(const unsigned int line, const unsigned int face, const unsigned char face_orientation) const
constexpr const ReferenceCell Pyramid
static const unsigned int invalid_unsigned_int
Definition: types.h:196
std::ostream & operator<<(std::ostream &out, const ReferenceCell &reference_cell)
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)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
bool standard_vs_true_line_orientation(const unsigned int line, const unsigned char face_orientation, const bool line_orientation) const
bool is_hyper_cube() const
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)
ArrayView< const unsigned int > faces_for_given_vertex(const unsigned int vertex_index) const
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)
constexpr const ReferenceCell & get_hypercube()
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
unsigned int standard_to_real_face_vertex(const unsigned int vertex, const unsigned int face, const unsigned char face_orientation) const
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition: mapping.cc:260
bool is_simplex() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1720
constexpr const ReferenceCell Triangle
static Point< dim, Number > unit_vector(const unsigned int i)
static ReferenceCell n_vertices_to_type(const int dim, const unsigned int n_vertices)
constexpr ::ReferenceCell make_reference_cell_from_int(const std::uint8_t kind)
static const char U
Tensor< 1, dim > unit_normal_vectors(const unsigned int face_no) const
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
double d_linear_shape_function(const Point< dim > &xi, const unsigned int i) const
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
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:2664
void serialize(Archive &archive, const unsigned int)
constexpr bool operator==(const ReferenceCell &type) const
constexpr const ReferenceCell Line
static std::array< unsigned int, 2 > standard_quad_vertex_to_line_vertex_index(const unsigned int vertex)
constexpr const ReferenceCell Wedge
Tensor< 1, dim > unit_tangential_vectors(const unsigned int face_no, const unsigned int i) const
static ::ExceptionBase & ExcMessage(std::string arg1)
NoPermutation(const ::ReferenceCell &entity_type, const std::array< T, N > &vertices_0, const std::array< T, N > &vertices_1)
unsigned int child_cell_on_face(const unsigned int face_n, const unsigned int subface_n, const unsigned char face_orientation=1) const
constexpr const ReferenceCell Tetrahedron
#define Assert(cond, exc)
Definition: exceptions.h:1461
std::uint8_t kind
constexpr bool operator!=(const ReferenceCell &type) const
Abstract base class for mapping classes.
Definition: mapping.h:310
boost::integer_range< IncrementableType > iota_view
Definition: iota_view.h:46
unsigned int vertex_indices[2]
Definition: grid_tools.cc:1256
std::array< unsigned int, 2 > standard_line_to_face_and_line_index(const unsigned int line) const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:407
std::array< T, N > permute_according_orientation(const std::array< T, N > &vertices, const unsigned int orientation) const
static unsigned int standard_to_real_line_vertex(const unsigned int vertex, const bool line_orientation=true)
std::string to_string(const T &t)
Definition: patterns.h:2330
constexpr const ReferenceCell Hexahedron
std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices() const
Point< 3 > vertices[4]
unsigned int n_vertices() const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > line_indices() const
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)
Expression fabs(const Expression &x)
unsigned char compute_orientation(const std::array< T, N > &vertices_0, const std::array< T, N > &vertices_1) const
Point< dim > vertex(const unsigned int v) const
virtual void print_info(std::ostream &out) const override
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
unsigned int get_dimension() const
unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const unsigned char face_orientation) const
unsigned int n_lines() const
unsigned int face_to_cell_lines(const unsigned int face, const unsigned int line, const unsigned char face_orientation) const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices() const
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
constexpr const ReferenceCell Invalid
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)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:406
constexpr const ReferenceCell Quadrilateral
static const char N
constexpr const ReferenceCell Vertex
std::array< unsigned int, 2 > standard_vertex_to_face_and_vertex_index(const unsigned int vertex) const
static std::array< unsigned int, 2 > standard_hex_vertex_to_quad_vertex_index(const unsigned int vertex)
static ::ExceptionBase & ExcNotImplemented()
Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i) const
unsigned int n_faces() const
std::istream & operator>>(std::istream &in, ReferenceCell &reference_cell)
constexpr ReferenceCell()
const std::array< T, N > vertices_1
const ::ReferenceCell entity_type
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2639
constexpr const ReferenceCell & get_simplex()
ReferenceCell face_reference_cell(const unsigned int face_no) const
static std::array< unsigned int, 2 > standard_hex_line_to_quad_line_index(const unsigned int line)
const std::array< T, N > vertices_0
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108
#define DEAL_II_CONSTEXPR
Definition: config.h:176
bool get_bit(const unsigned char number, const unsigned int n)
Definition: utilities.h:1381
static ::ExceptionBase & ExcInternalError()