Reference documentation for deal.II version Git 1aa49a030b 2021-12-07 05:12:24 -0500
\(\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 unsigned char 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 ? 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  }
1456  else if (*this == ReferenceCells::Tetrahedron)
1457  {
1458  static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 1}},
1459  {{0, 1, 2}},
1460  {{2, 1, 0}},
1461  {{1, 2, 0}},
1462  {{1, 0, 2}},
1463  {{2, 0, 1}}}};
1464 
1465  return table[face_orientation][vertex];
1466  }
1467  else if (*this == ReferenceCells::Pyramid)
1468  {
1469  if (face == 0) // The quadrilateral face
1470  {
1472  vertex,
1473  Utilities::get_bit(face_orientation, 0),
1474  Utilities::get_bit(face_orientation, 2),
1475  Utilities::get_bit(face_orientation, 1));
1476  }
1477  else // One of the triangular faces
1478  {
1479  static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 1}},
1480  {{0, 1, 2}},
1481  {{2, 1, 0}},
1482  {{1, 2, 0}},
1483  {{1, 0, 2}},
1484  {{2, 0, 1}}}};
1485 
1486  return table[face_orientation][vertex];
1487  }
1488  }
1489  else if (*this == ReferenceCells::Wedge)
1490  {
1491  if (face > 1) // One of the quadrilateral faces
1492  {
1494  vertex,
1495  Utilities::get_bit(face_orientation, 0),
1496  Utilities::get_bit(face_orientation, 2),
1497  Utilities::get_bit(face_orientation, 1));
1498  }
1499  else // One of the triangular faces
1500  {
1501  static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 1}},
1502  {{0, 1, 2}},
1503  {{2, 1, 0}},
1504  {{1, 2, 0}},
1505  {{1, 0, 2}},
1506  {{2, 0, 1}}}};
1507 
1508  return table[face_orientation][vertex];
1509  }
1510  }
1511  else if (*this == ReferenceCells::Hexahedron)
1512  {
1514  vertex,
1515  Utilities::get_bit(face_orientation, 0),
1516  Utilities::get_bit(face_orientation, 2),
1517  Utilities::get_bit(face_orientation, 1));
1518  }
1519 
1520  Assert(false, ExcNotImplemented());
1522 }
1523 
1524 
1525 
1526 inline unsigned int
1528  const unsigned int line,
1529  const unsigned int face,
1530  const unsigned char face_orientation) const
1531 {
1532  AssertIndexRange(face, n_faces());
1534 
1535  if (*this == ReferenceCells::Vertex)
1536  {
1537  Assert(false, ExcNotImplemented());
1538  }
1539  else if (*this == ReferenceCells::Line)
1540  {
1541  Assert(false, ExcNotImplemented());
1542  }
1543  else if (*this == ReferenceCells::Triangle)
1544  {
1545  Assert(false, ExcNotImplemented());
1546  }
1547  else if (*this == ReferenceCells::Quadrilateral)
1548  {
1549  Assert(false, ExcNotImplemented());
1550  }
1551  else if (*this == ReferenceCells::Tetrahedron)
1552  {
1553  static const ndarray<unsigned int, 6, 3> table = {{{{2, 1, 0}},
1554  {{0, 1, 2}},
1555  {{1, 0, 2}},
1556  {{1, 2, 0}},
1557  {{0, 2, 1}},
1558  {{2, 0, 1}}}};
1559 
1560  return table[face_orientation][line];
1561  }
1562  else if (*this == ReferenceCells::Pyramid)
1563  {
1564  if (face == 0) // The quadrilateral face
1565  {
1567  line,
1568  Utilities::get_bit(face_orientation, 0),
1569  Utilities::get_bit(face_orientation, 2),
1570  Utilities::get_bit(face_orientation, 1));
1571  }
1572  else // One of the triangular faces
1573  {
1574  static const ndarray<unsigned int, 6, 3> table = {{{{2, 1, 0}},
1575  {{0, 1, 2}},
1576  {{1, 0, 2}},
1577  {{1, 2, 0}},
1578  {{0, 2, 1}},
1579  {{2, 0, 1}}}};
1580 
1581  return table[face_orientation][line];
1582  }
1583  }
1584  else if (*this == ReferenceCells::Wedge)
1585  {
1586  if (face > 1) // One of the quadrilateral faces
1587  {
1589  line,
1590  Utilities::get_bit(face_orientation, 0),
1591  Utilities::get_bit(face_orientation, 2),
1592  Utilities::get_bit(face_orientation, 1));
1593  }
1594  else // One of the triangular faces
1595  {
1596  static const ndarray<unsigned int, 6, 3> table = {{{{2, 1, 0}},
1597  {{0, 1, 2}},
1598  {{1, 0, 2}},
1599  {{1, 2, 0}},
1600  {{0, 2, 1}},
1601  {{2, 0, 1}}}};
1602 
1603  return table[face_orientation][line];
1604  }
1605  }
1606  else if (*this == ReferenceCells::Hexahedron)
1607  {
1609  line,
1610  Utilities::get_bit(face_orientation, 0),
1611  Utilities::get_bit(face_orientation, 2),
1612  Utilities::get_bit(face_orientation, 1));
1613  }
1614 
1615  Assert(false, ExcNotImplemented());
1617 }
1618 
1619 
1620 
1621 namespace ReferenceCells
1622 {
1623  template <int dim>
1624  inline constexpr const ReferenceCell &
1626  {
1627  switch (dim)
1628  {
1629  case 0:
1630  return ReferenceCells::Vertex;
1631  case 1:
1632  return ReferenceCells::Line;
1633  case 2:
1634  return ReferenceCells::Triangle;
1635  case 3:
1637  default:
1638  Assert(false, ExcNotImplemented());
1639  return ReferenceCells::Invalid;
1640  }
1641  }
1642 
1643 
1644 
1645  template <int dim>
1646  inline constexpr const ReferenceCell &
1648  {
1649  switch (dim)
1650  {
1651  case 0:
1652  return ReferenceCells::Vertex;
1653  case 1:
1654  return ReferenceCells::Line;
1655  case 2:
1657  case 3:
1659  default:
1660  Assert(false, ExcNotImplemented());
1661  return ReferenceCells::Invalid;
1662  }
1663  }
1664 } // namespace ReferenceCells
1665 
1666 
1667 inline ReferenceCell
1668 ReferenceCell::n_vertices_to_type(const int dim, const unsigned int n_vertices)
1669 {
1670  AssertIndexRange(dim, 4);
1671  AssertIndexRange(n_vertices, 9);
1672 
1673  const auto X = ReferenceCells::Invalid;
1674  static const ndarray<ReferenceCell, 4, 9> table = {
1675  {// dim 0
1676  {{X, ReferenceCells::Vertex, X, X, X, X, X, X, X}},
1677  // dim 1
1678  {{X, X, ReferenceCells::Line, X, X, X, X, X, X}},
1679  // dim 2
1680  {{X,
1681  X,
1682  X,
1685  X,
1686  X,
1687  X,
1688  X}},
1689  // dim 3
1690  {{X,
1691  X,
1692  X,
1693  X,
1697  X,
1699  Assert(table[dim][n_vertices] != ReferenceCells::Invalid,
1700  ExcMessage("The combination of dim = " + std::to_string(dim) +
1701  " and n_vertices = " + std::to_string(n_vertices) +
1702  " does not correspond to a known reference cell type."));
1703  return table[dim][n_vertices];
1704 }
1705 
1706 
1707 
1708 template <int dim>
1709 inline double
1711  const unsigned int i) const
1712 {
1714  if (*this == ReferenceCells::get_hypercube<dim>())
1716 
1717  if (*this ==
1718  ReferenceCells::Triangle) // see also
1719  // BarycentricPolynomials<2>::compute_value
1720  {
1721  switch (i)
1722  {
1723  case 0:
1724  return 1.0 - xi[std::min(0, dim - 1)] - xi[std::min(1, dim - 1)];
1725  case 1:
1726  return xi[std::min(0, dim - 1)];
1727  case 2:
1728  return xi[std::min(1, dim - 1)];
1729  }
1730  }
1731 
1732  if (*this ==
1733  ReferenceCells::Tetrahedron) // see also
1734  // BarycentricPolynomials<3>::compute_value
1735  {
1736  switch (i)
1737  {
1738  case 0:
1739  return 1.0 - xi[std::min(0, dim - 1)] - xi[std::min(1, dim - 1)] -
1740  xi[std::min(2, dim - 1)];
1741  case 1:
1742  return xi[std::min(0, dim - 1)];
1743  case 2:
1744  return xi[std::min(1, dim - 1)];
1745  case 3:
1746  return xi[std::min(2, dim - 1)];
1747  }
1748  }
1749 
1750  if (*this ==
1751  ReferenceCells::Wedge) // see also
1752  // ScalarLagrangePolynomialWedge::compute_value
1753  {
1755  .d_linear_shape_function<2>(Point<2>(xi[std::min(0, dim - 1)],
1756  xi[std::min(1, dim - 1)]),
1757  i % 3) *
1759  .d_linear_shape_function<1>(Point<1>(xi[std::min(2, dim - 1)]),
1760  i / 3);
1761  }
1762 
1763  if (*this ==
1764  ReferenceCells::Pyramid) // see also
1765  // ScalarLagrangePolynomialPyramid::compute_value
1766  {
1767  const double Q14 = 0.25;
1768  double ration;
1769 
1770  const double r = xi[std::min(0, dim - 1)];
1771  const double s = xi[std::min(1, dim - 1)];
1772  const double t = xi[std::min(2, dim - 1)];
1773 
1774  if (fabs(t - 1.0) > 1.0e-14)
1775  {
1776  ration = (r * s * t) / (1.0 - t);
1777  }
1778  else
1779  {
1780  ration = 0.0;
1781  }
1782 
1783  if (i == 0)
1784  return Q14 * ((1.0 - r) * (1.0 - s) - t + ration);
1785  if (i == 1)
1786  return Q14 * ((1.0 + r) * (1.0 - s) - t - ration);
1787  if (i == 2)
1788  return Q14 * ((1.0 - r) * (1.0 + s) - t - ration);
1789  if (i == 3)
1790  return Q14 * ((1.0 + r) * (1.0 + s) - t + ration);
1791  else
1792  return t;
1793  }
1794 
1795  Assert(false, ExcNotImplemented());
1796 
1797  return 0.0;
1798 }
1799 
1800 
1801 
1802 template <int dim>
1803 inline Tensor<1, dim>
1805  const unsigned int i) const
1806 {
1808  if (*this == ReferenceCells::get_hypercube<dim>())
1810 
1811  if (*this ==
1812  ReferenceCells::Triangle) // see also
1813  // BarycentricPolynomials<2>::compute_grad
1814  {
1815  switch (i)
1816  {
1817  case 0:
1818  return Point<dim>(-1.0, -1.0);
1819  case 1:
1820  return Point<dim>(+1.0, +0.0);
1821  case 2:
1822  return Point<dim>(+0.0, +1.0);
1823  }
1824  }
1825 
1826  Assert(false, ExcNotImplemented());
1827 
1828  return Point<dim>(+0.0, +0.0, +0.0);
1829 }
1830 
1831 
1832 template <int dim>
1833 inline Tensor<1, dim>
1834 ReferenceCell::unit_tangential_vectors(const unsigned int face_no,
1835  const unsigned int i) const
1836 {
1838  AssertIndexRange(i, dim - 1);
1839 
1840  if (*this == ReferenceCells::get_hypercube<dim>())
1841  {
1843  return GeometryInfo<dim>::unit_tangential_vectors[face_no][i];
1844  }
1845  else if (*this == ReferenceCells::Triangle)
1846  {
1847  AssertIndexRange(face_no, 3);
1848  static const std::array<Tensor<1, dim>, 3> table = {
1849  {Point<dim>(1, 0),
1850  Point<dim>(-std::sqrt(0.5), +std::sqrt(0.5)),
1851  Point<dim>(0, -1)}};
1852 
1853  return table[face_no];
1854  }
1855  else if (*this == ReferenceCells::Tetrahedron)
1856  {
1857  AssertIndexRange(face_no, 4);
1858  static const ndarray<Tensor<1, dim>, 4, 2> table = {
1859  {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
1860  {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
1861  {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
1862  {{Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
1863  +std::pow(1.0 / 3.0, 1.0 / 4.0),
1864  0),
1865  Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
1866  0,
1867  +std::pow(1.0 / 3.0, 1.0 / 4.0))}}}};
1868 
1869  return table[face_no][i];
1870  }
1871  else if (*this == ReferenceCells::Wedge)
1872  {
1873  AssertIndexRange(face_no, 5);
1874  static const ndarray<Tensor<1, dim>, 5, 2> table = {
1875  {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
1876  {{Point<dim>(1, 0, 0), Point<dim>(0, 1, 0)}},
1877  {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
1878  {{Point<dim>(-1 / std::sqrt(2.0), +1 / std::sqrt(2.0), 0),
1879  Point<dim>(0, 0, 1)}},
1880  {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}}}};
1881 
1882  return table[face_no][i];
1883  }
1884  else if (*this == ReferenceCells::Pyramid)
1885  {
1886  AssertIndexRange(face_no, 5);
1887  static const ndarray<Tensor<1, dim>, 5, 2> table = {
1888  {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
1889  {{Point<dim>(+1.0 / sqrt(2.0), 0, +1.0 / sqrt(2.0)),
1890  Point<dim>(0, 1, 0)}},
1891  {{Point<dim>(+1.0 / sqrt(2.0), 0, -1.0 / sqrt(2.0)),
1892  Point<dim>(0, 1, 0)}},
1893  {{Point<dim>(1, 0, 0),
1894  Point<dim>(0, +1.0 / sqrt(2.0), +1.0 / sqrt(2.0))}},
1895  {{Point<dim>(1, 0, 0),
1896  Point<dim>(0, +1.0 / sqrt(2.0), -1.0 / sqrt(2.0))}}}};
1897 
1898  return table[face_no][i];
1899  }
1900 
1901  Assert(false, ExcNotImplemented());
1902 
1903  return {};
1904 }
1905 
1906 
1907 
1908 template <int dim>
1909 inline Tensor<1, dim>
1910 ReferenceCell::unit_normal_vectors(const unsigned int face_no) const
1911 {
1912  AssertDimension(dim, this->get_dimension());
1913 
1914  if (is_hyper_cube())
1915  {
1917  return GeometryInfo<dim>::unit_normal_vector[face_no];
1918  }
1919  else if (dim == 2)
1920  {
1922 
1923  // Return the rotated vector
1924  return cross_product_2d(unit_tangential_vectors<dim>(face_no, 0));
1925  }
1926  else if (dim == 3)
1927  {
1928  return cross_product_3d(unit_tangential_vectors<dim>(face_no, 0),
1929  unit_tangential_vectors<dim>(face_no, 1));
1930  }
1931 
1932  Assert(false, ExcNotImplemented());
1933 
1934  return {};
1935 }
1936 
1937 
1938 
1939 inline bool
1941  const unsigned int line,
1942  const unsigned char face_orientation_raw,
1943  const unsigned char line_orientation) const
1944 {
1945  if (*this == ReferenceCells::Hexahedron)
1946  {
1947  static const bool bool_table[2][2][2][2] = {
1948  {{{true, false}, // lines 0/1, face_orientation=false,
1949  // face_flip=false, face_rotation=false and true
1950  {false, true}}, // lines 0/1, face_orientation=false,
1951  // face_flip=true, face_rotation=false and true
1952  {{true, true}, // lines 0/1, face_orientation=true,
1953  // face_flip=false, face_rotation=false and true
1954  {false, false}}}, // lines 0/1, face_orientation=true,
1955  // face_flip=true, face_rotation=false and true
1956 
1957  {{{true, true}, // lines 2/3 ...
1958  {false, false}},
1959  {{true, false}, {false, true}}}};
1960 
1961  const bool face_orientation = Utilities::get_bit(face_orientation_raw, 0);
1962  const bool face_flip = Utilities::get_bit(face_orientation_raw, 2);
1963  const bool face_rotation = Utilities::get_bit(face_orientation_raw, 1);
1964 
1965  return (static_cast<bool>(line_orientation) ==
1966  bool_table[line / 2][face_orientation][face_flip][face_rotation]);
1967  }
1968  else
1969  // TODO: This might actually be wrong for some of the other
1970  // kinds of objects. We should check this
1971  return true;
1972 }
1973 
1974 
1975 
1976 namespace internal
1977 {
1978  template <typename T, std::size_t N>
1980  {
1981  public:
1985  NoPermutation(const ::ReferenceCell &entity_type,
1986  const std::array<T, N> & vertices_0,
1987  const std::array<T, N> & vertices_1)
1988  : entity_type(entity_type)
1989  , vertices_0(vertices_0)
1990  , vertices_1(vertices_1)
1991  {}
1992 
1996  virtual ~NoPermutation() noexcept override = default;
1997 
2001  virtual void
2002  print_info(std::ostream &out) const override
2003  {
2004  out << "[";
2005 
2006  const unsigned int n_vertices = entity_type.n_vertices();
2007 
2008  for (unsigned int i = 0; i < n_vertices; ++i)
2009  {
2010  out << vertices_0[i];
2011  if (i + 1 != n_vertices)
2012  out << ",";
2013  }
2014 
2015  out << "] is not a permutation of [";
2016 
2017  for (unsigned int i = 0; i < n_vertices; ++i)
2018  {
2019  out << vertices_1[i];
2020  if (i + 1 != n_vertices)
2021  out << ",";
2022  }
2023 
2024  out << "]." << std::endl;
2025  }
2026 
2030  const ::ReferenceCell entity_type;
2031 
2035  const std::array<T, N> vertices_0;
2036 
2040  const std::array<T, N> vertices_1;
2041  };
2042 } // namespace internal
2043 
2044 
2045 
2046 template <typename T, std::size_t N>
2047 inline unsigned char
2048 ReferenceCell::compute_orientation(const std::array<T, N> &vertices_0,
2049  const std::array<T, N> &vertices_1) const
2050 {
2051  AssertIndexRange(n_vertices(), N + 1);
2052  if (*this == ReferenceCells::Line)
2053  {
2054  const std::array<T, 2> i{{vertices_0[0], vertices_0[1]}};
2055  const std::array<T, 2> j{{vertices_1[0], vertices_1[1]}};
2056 
2057  // line_orientation=true
2058  if (i == std::array<T, 2>{{j[0], j[1]}})
2059  return 1;
2060 
2061  // line_orientation=false
2062  if (i == std::array<T, 2>{{j[1], j[0]}})
2063  return 0;
2064  }
2065  else if (*this == ReferenceCells::Triangle)
2066  {
2067  const std::array<T, 3> i{{vertices_0[0], vertices_0[1], vertices_0[2]}};
2068  const std::array<T, 3> j{{vertices_1[0], vertices_1[1], vertices_1[2]}};
2069 
2070  // face_orientation=true, face_rotation=false, face_flip=false
2071  if (i == std::array<T, 3>{{j[0], j[1], j[2]}})
2072  return 1;
2073 
2074  // face_orientation=true, face_rotation=true, face_flip=false
2075  if (i == std::array<T, 3>{{j[1], j[2], j[0]}})
2076  return 3;
2077 
2078  // face_orientation=true, face_rotation=false, face_flip=true
2079  if (i == std::array<T, 3>{{j[2], j[0], j[1]}})
2080  return 5;
2081 
2082  // face_orientation=false, face_rotation=false, face_flip=false
2083  if (i == std::array<T, 3>{{j[0], j[2], j[1]}})
2084  return 0;
2085 
2086  // face_orientation=false, face_rotation=true, face_flip=false
2087  if (i == std::array<T, 3>{{j[2], j[1], j[0]}})
2088  return 2;
2089 
2090  // face_orientation=false, face_rotation=false, face_flip=true
2091  if (i == std::array<T, 3>{{j[1], j[0], j[2]}})
2092  return 4;
2093  }
2094  else if (*this == ReferenceCells::Quadrilateral)
2095  {
2096  const std::array<T, 4> i{
2097  {vertices_0[0], vertices_0[1], vertices_0[2], vertices_0[3]}};
2098  const std::array<T, 4> j{
2099  {vertices_1[0], vertices_1[1], vertices_1[2], vertices_1[3]}};
2100 
2101  // face_orientation=true, face_rotation=false, face_flip=false
2102  if (i == std::array<T, 4>{{j[0], j[1], j[2], j[3]}})
2103  return 1;
2104 
2105  // face_orientation=true, face_rotation=true, face_flip=false
2106  if (i == std::array<T, 4>{{j[2], j[0], j[3], j[1]}})
2107  return 3;
2108 
2109  // face_orientation=true, face_rotation=false, face_flip=true
2110  if (i == std::array<T, 4>{{j[3], j[2], j[1], j[0]}})
2111  return 5;
2112 
2113  // face_orientation=true, face_rotation=true, face_flip=true
2114  if (i == std::array<T, 4>{{j[1], j[3], j[0], j[2]}})
2115  return 7;
2116 
2117  // face_orientation=false, face_rotation=false, face_flip=false
2118  if (i == std::array<T, 4>{{j[0], j[2], j[1], j[3]}})
2119  return 0;
2120 
2121  // face_orientation=false, face_rotation=true, face_flip=false
2122  if (i == std::array<T, 4>{{j[2], j[3], j[0], j[1]}})
2123  return 2;
2124 
2125  // face_orientation=false, face_rotation=false, face_flip=true
2126  if (i == std::array<T, 4>{{j[3], j[1], j[2], j[0]}})
2127  return 4;
2128 
2129  // face_orientation=false, face_rotation=true, face_flip=true
2130  if (i == std::array<T, 4>{{j[1], j[0], j[3], j[2]}})
2131  return 6;
2132  }
2133 
2134  Assert(false, (internal::NoPermutation<T, N>(*this, vertices_0, vertices_1)));
2135 
2136  return -1;
2137 }
2138 
2139 
2140 
2141 template <typename T, std::size_t N>
2142 inline std::array<T, N>
2144  const std::array<T, N> &vertices,
2145  const unsigned int orientation) const
2146 {
2147  std::array<T, 4> temp;
2148 
2149  if (*this == ReferenceCells::Line)
2150  {
2151  switch (orientation)
2152  {
2153  case 1:
2154  temp = {{vertices[0], vertices[1]}};
2155  break;
2156  case 0:
2157  temp = {{vertices[1], vertices[0]}};
2158  break;
2159  default:
2160  Assert(false, ExcNotImplemented());
2161  }
2162  }
2163  else if (*this == ReferenceCells::Triangle)
2164  {
2165  switch (orientation)
2166  {
2167  case 1:
2168  temp = {{vertices[0], vertices[1], vertices[2]}};
2169  break;
2170  case 3:
2171  temp = {{vertices[1], vertices[2], vertices[0]}};
2172  break;
2173  case 5:
2174  temp = {{vertices[2], vertices[0], vertices[1]}};
2175  break;
2176  case 0:
2177  temp = {{vertices[0], vertices[2], vertices[1]}};
2178  break;
2179  case 2:
2180  temp = {{vertices[2], vertices[1], vertices[0]}};
2181  break;
2182  case 4:
2183  temp = {{vertices[1], vertices[0], vertices[2]}};
2184  break;
2185  default:
2186  Assert(false, ExcNotImplemented());
2187  }
2188  }
2189  else if (*this == ReferenceCells::Quadrilateral)
2190  {
2191  switch (orientation)
2192  {
2193  case 1:
2194  temp = {{vertices[0], vertices[1], vertices[2], vertices[3]}};
2195  break;
2196  case 3:
2197  temp = {{vertices[2], vertices[0], vertices[3], vertices[1]}};
2198  break;
2199  case 5:
2200  temp = {{vertices[3], vertices[2], vertices[1], vertices[0]}};
2201  break;
2202  case 7:
2203  temp = {{vertices[1], vertices[3], vertices[0], vertices[2]}};
2204  break;
2205  case 0:
2206  temp = {{vertices[0], vertices[2], vertices[1], vertices[3]}};
2207  break;
2208  case 2:
2209  temp = {{vertices[2], vertices[3], vertices[0], vertices[1]}};
2210  break;
2211  case 4:
2212  temp = {{vertices[3], vertices[1], vertices[2], vertices[0]}};
2213  break;
2214  case 6:
2215  temp = {{vertices[1], vertices[0], vertices[3], vertices[2]}};
2216  break;
2217  default:
2218  Assert(false, ExcNotImplemented());
2219  }
2220  }
2221  else
2222  {
2223  AssertThrow(false, ExcNotImplemented());
2224  }
2225 
2226  std::array<T, N> temp_;
2227  std::copy_n(temp.begin(), N, temp_.begin());
2228 
2229  return temp_;
2230 }
2231 
2232 
2234 
2235 #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 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:254
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:2650
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
bool standard_vs_true_line_orientation(const unsigned int line, const unsigned char face_orientation, const unsigned char line_orientation) 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:404
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:2329
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:403
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:2625
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:173
bool get_bit(const unsigned char number, const unsigned int n)
Definition: utilities.h:1381
static ::ExceptionBase & ExcInternalError()