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