Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-2684-gc61376a70f 2025-02-22 15:30:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
reference_cell.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2020 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_tria_reference_cell_h
16#define dealii_tria_reference_cell_h
17
18#include <deal.II/base/config.h>
19
23#include <deal.II/base/point.h>
24#include <deal.II/base/tensor.h>
25#include <deal.II/base/types.h>
27
29
30#include <boost/container/small_vector.hpp>
31
32#include <iosfwd>
33#include <string>
34
36
37// Forward declarations
38#ifndef DOXYGEN
39template <int dim, int spacedim>
40class Mapping;
41
42template <int dim>
43class Quadrature;
44
45class ReferenceCell;
46#endif
47
48
49namespace internal
50{
63 constexpr ReferenceCell
64 make_reference_cell_from_int(const std::uint8_t kind);
65} // namespace internal
66
67
68
75enum class IsotropicRefinementChoice : std::uint8_t
76{
84 cut_tet_68 = 1,
88 cut_tet_57 = 2,
92 cut_tet_49 = 3,
93};
94
95
96
128{
129public:
137 static ReferenceCell
138 n_vertices_to_type(const int dim, const unsigned int n_vertices);
139
148 constexpr ReferenceCell();
149
160 bool
161 is_hyper_cube() const;
162
166 bool
167 is_simplex() const;
168
173 unsigned int
174 get_dimension() const;
175
193 template <int dim>
194 double
195 d_linear_shape_function(const Point<dim> &xi, const unsigned int i) const;
196
205 template <int dim>
208 const unsigned int i) const;
209
219 template <int dim, int spacedim = dim>
220 std::unique_ptr<Mapping<dim, spacedim>>
221 get_default_mapping(const unsigned int degree) const;
222
234 template <int dim, int spacedim = dim>
237
250 template <int dim>
252 get_gauss_type_quadrature(const unsigned n_points_1d) const;
253
267 template <int dim>
270
288 template <int dim>
289 const Quadrature<dim> &
291
306 unsigned int
307 n_vertices() const;
308
314 vertex_indices() const;
315
327 template <int dim>
329 vertex(const unsigned int v) const;
330
336 unsigned int
337 n_lines() const;
338
344 line_indices() const;
345
351 unsigned int
352 n_faces() const;
353
359 face_indices() const;
360
364 unsigned int
366
371 get_isotropic_refinement_choice(const unsigned int ref_choice) const;
372
384 unsigned int
385 n_isotropic_children() const;
386
395 template <int dim>
396 unsigned int
397 n_children(const RefinementCase<dim> ref_case =
399
406
420 face_reference_cell(const unsigned int face_no) const;
421
437 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
438 "Use numbers::default_geometric_orientation instead.")
439 static constexpr types::geometric_orientation
441
450 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
451 "Use numbers::reverse_line_orientation instead.")
452 static constexpr types::geometric_orientation
454
462 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
463 "Use the version of this function which takes the orientation as the third "
464 "parameter instead.")
465 unsigned int
466 child_cell_on_face(const unsigned int face, const unsigned int subface) const;
467
509 unsigned int
510 child_cell_on_face(const unsigned int face,
511 const unsigned int subface,
512 const types::geometric_orientation face_orientation) const;
513
523 std::array<unsigned int, 2>
524 standard_vertex_to_face_and_vertex_index(const unsigned int vertex) const;
525
535 std::array<unsigned int, 2>
536 standard_line_to_face_and_line_index(const unsigned int line) const;
537
546 unsigned int
547 line_to_cell_vertices(const unsigned int line,
548 const unsigned int vertex) const;
549
553 unsigned int
554 face_to_cell_lines(const unsigned int face,
555 const unsigned int line,
556 const types::geometric_orientation face_orientation) const;
557
561 unsigned int
563 const unsigned int face,
564 const unsigned int vertex,
565 const types::geometric_orientation face_orientation) const;
566
590 template <int dim>
592 face_vertex_location(const unsigned int face,
593 const unsigned int vertex) const;
594
598 unsigned int
600 const unsigned int vertex,
601 const unsigned int face,
602 const types::geometric_orientation face_orientation) const;
603
607 unsigned int
609 const unsigned int line,
610 const unsigned int face,
611 const types::geometric_orientation face_orientation) const;
612
620 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
621 "Use face_to_cell_line_orientation() instead.")
624 const unsigned int line,
625 const unsigned int face,
626 const types::geometric_orientation face_orientation,
627 const types::geometric_orientation line_orientation) const;
628
653 const unsigned int face_line_no,
654 const unsigned int face_no,
655 const types::geometric_orientation face_orientation,
656 const types::geometric_orientation line_orientation) const;
657
681 double
682 volume() const;
683
700 template <int dim>
702 barycenter() const;
703
727 template <int dim>
728 bool
729 contains_point(const Point<dim> &p, const double tolerance = 0) const;
730
739 template <int dim>
741 closest_point(const Point<dim> &p) const;
742
754 template <int dim>
756 face_tangent_vector(const unsigned int face_no, const unsigned int i) const;
757
771 template <int dim>
772 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT("Use face_tangent_vector() instead.")
773 Tensor<1, dim> unit_tangential_vectors(const unsigned int face_no,
774 const unsigned int i) const;
775
783 template <int dim>
785 face_normal_vector(const unsigned int face_no) const;
786
796 template <int dim>
797 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT("Use face_normal_vector() instead.")
798 Tensor<1, dim> unit_normal_vectors(const unsigned int face_no) const;
799
805 unsigned int
806 n_face_orientations(const unsigned int face_no) const;
807
824 template <typename T, std::size_t N>
826 compute_orientation(const std::array<T, N> &vertices_0,
827 const std::array<T, N> &vertices_1) const;
828
854 template <typename T>
857 const ArrayView<const T> &vertices_1) const;
858
859
873 template <typename T, std::size_t N>
874 DEAL_II_DEPRECATED std::array<T, N>
875 permute_according_orientation(const std::array<T, N> &vertices,
876 const unsigned int orientation) const;
877
889 template <typename T>
890 boost::container::small_vector<T, 8>
892 const ArrayView<const T> &vertices,
893 const types::geometric_orientation orientation) const;
894
902 const types::geometric_orientation orientation) const;
903
908 faces_for_given_vertex(const unsigned int vertex_index) const;
909
914 constexpr ::ndarray<unsigned int, 12, 4>
915 new_isotropic_child_face_lines(const unsigned int refinement_choice) const;
916
921 constexpr ::ndarray<unsigned int, 12, 4, 2>
923 const unsigned int refinement_choice) const;
924
929 constexpr ::ndarray<unsigned int, 8, 6>
930 new_isotropic_child_cell_faces(const unsigned int refinement_choice) const;
931
936 constexpr ::ndarray<unsigned int, 8, 4>
937 new_isotropic_child_cell_vertices(const unsigned int refinement_choice) const;
938
951 unsigned int
952 exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const;
953
957 unsigned int
958 exodusii_face_to_deal_face(const unsigned int face_n) const;
959
963 unsigned int
964 unv_vertex_to_deal_vertex(const unsigned int vertex_n) const;
965
969 unsigned int
970 vtk_linear_type() const;
971
976 unsigned int
977 vtk_quadratic_type() const;
978
983 unsigned int
984 vtk_lagrange_type() const;
985
1002 template <int dim>
1003 unsigned int
1005 const std::array<unsigned, dim> &node_indices,
1006 const std::array<unsigned, dim> &nodes_per_direction,
1007 const bool legacy_format) const;
1008
1012 unsigned int
1013 vtk_vertex_to_deal_vertex(const unsigned int vertex_index) const;
1014
1018 unsigned int
1019 gmsh_element_type() const;
1020
1034 std::string
1035 to_string() const;
1036
1040 constexpr operator std::uint8_t() const;
1041
1045 constexpr bool
1046 operator==(const ReferenceCell &type) const;
1047
1051 constexpr bool
1052 operator!=(const ReferenceCell &type) const;
1053
1059 template <class Archive>
1060 void
1061 serialize(Archive &archive, const unsigned int /*version*/);
1062
1066 static constexpr std::size_t
1068
1073private:
1077 std::uint8_t kind;
1078
1085 constexpr ReferenceCell(const std::uint8_t kind);
1086
1093 {{{0}}}};
1094
1099 {{{1, 0}}, {{0, 1}}}};
1100
1105 {{{0, 2, 1}},
1106 {{0, 1, 2}},
1107 {{2, 1, 0}},
1108 {{2, 0, 1}},
1109 {{1, 0, 2}},
1110 {{1, 2, 0}}}};
1111
1115 static constexpr ndarray<unsigned int, 8, 4>
1117 {{0, 1, 2, 3}},
1118 {{2, 3, 0, 1}},
1119 {{2, 0, 3, 1}},
1120 {{3, 1, 2, 0}},
1121 {{3, 2, 1, 0}},
1122 {{1, 0, 3, 2}},
1123 {{1, 3, 0, 2}}}};
1124
1129 friend constexpr ReferenceCell
1131
1132 friend std::ostream &
1133 operator<<(std::ostream &out, const ReferenceCell &reference_cell);
1134
1135 friend std::istream &
1136 operator>>(std::istream &in, ReferenceCell &reference_cell);
1137};
1138
1139
1147std::ostream &
1148operator<<(std::ostream &out, const ReferenceCell &reference_cell);
1149
1157std::istream &
1158operator>>(std::istream &in, ReferenceCell &reference_cell);
1159
1160
1161inline constexpr ReferenceCell::ReferenceCell(const std::uint8_t kind)
1162 : kind(kind)
1163{}
1164
1165
1166
1167inline constexpr ReferenceCell::operator std::uint8_t() const
1168{
1169 return kind;
1170}
1171
1172
1173
1174inline constexpr bool
1176{
1177 return kind == type.kind;
1178}
1179
1180
1181
1182inline constexpr bool
1184{
1185 return kind != type.kind;
1186}
1187
1188
1189
1190namespace internal
1191{
1192 inline constexpr ReferenceCell
1193 make_reference_cell_from_int(const std::uint8_t kind)
1194 {
1195#ifndef DEAL_II_CXX14_CONSTEXPR_BUG
1196 // Make sure these are the only indices from which objects can be
1197 // created.
1198 Assert((kind == std::numeric_limits<std::uint8_t>::max()) || (kind < 8),
1200#endif
1201
1202 // Call the private constructor, which we can from here because this
1203 // function is a 'friend'.
1204 return {kind};
1205 }
1206} // namespace internal
1207
1208
1209
1218{
1219 constexpr const ReferenceCell Vertex =
1221 constexpr const ReferenceCell Line =
1223 constexpr const ReferenceCell Triangle =
1227 constexpr const ReferenceCell Tetrahedron =
1229 constexpr const ReferenceCell Pyramid =
1231 constexpr const ReferenceCell Wedge =
1233 constexpr const ReferenceCell Hexahedron =
1235 constexpr const ReferenceCell Invalid =
1237 std::numeric_limits<std::uint8_t>::max());
1238
1244 template <int dim>
1245 constexpr const ReferenceCell &
1246 get_simplex();
1247
1253 template <int dim>
1254 constexpr const ReferenceCell &
1255 get_hypercube();
1256
1264 template <int structdim>
1265 constexpr unsigned int
1267
1275 template <int structdim>
1276 constexpr unsigned int
1277 max_n_lines();
1278
1294 template <int structdim>
1295 constexpr unsigned int
1296 max_n_faces();
1297} // namespace ReferenceCells
1298
1299
1300
1302 : ReferenceCell(ReferenceCells::Invalid)
1303{}
1304
1305
1306
1307template <class Archive>
1308inline void
1309ReferenceCell::serialize(Archive &archive, const unsigned int /*version*/)
1310{
1311 archive &kind;
1312}
1313
1314
1315
1316inline constexpr std::size_t
1321
1322
1323
1325ReferenceCell::faces_for_given_vertex(const unsigned int vertex) const
1326{
1328 switch (this->kind)
1329 {
1331 return {&GeometryInfo<2>::vertex_to_face[vertex][0], 1};
1333 return {&GeometryInfo<2>::vertex_to_face[vertex][0], 2};
1335 {
1336 static constexpr ndarray<unsigned int, 3, 2> table = {
1337 {{{0, 2}}, {{0, 1}}, {{1, 2}}}};
1338 return table[vertex];
1339 }
1341 {
1342 static constexpr ndarray<unsigned int, 4, 3> table = {
1343 {{{0, 1, 2}}, {{0, 1, 3}}, {{0, 2, 3}}, {{1, 2, 3}}}};
1344
1345 return table[vertex];
1346 }
1348 {
1349 static constexpr unsigned int X = numbers::invalid_unsigned_int;
1350 static constexpr ndarray<unsigned int, 5, 4> table = {
1351 {{{0, 1, 3, X}},
1352 {{0, 2, 3, X}},
1353 {{0, 1, 4, X}},
1354 {{0, 2, 4, X}},
1355 {{1, 2, 3, 4}}}};
1356
1357 return {&table[vertex][0], vertex == 4 ? 4u : 3u};
1358 }
1360 {
1362 static constexpr ndarray<unsigned int, 6, 3> table = {{{{0, 2, 4}},
1363 {{0, 2, 3}},
1364 {{0, 3, 4}},
1365 {{1, 2, 4}},
1366 {{1, 2, 3}},
1367 {{1, 3, 4}}}};
1368
1369 return table[vertex];
1370 }
1372 return {&GeometryInfo<3>::vertex_to_face[vertex][0], 3};
1373 default:
1375 }
1376
1377 return {};
1378}
1379
1380
1381constexpr ::ndarray<unsigned int, 12, 4>
1383 const unsigned int refienement_choice) const
1384{
1385 // AssertIndexRange(refienement_choice, n_isotropic_refinement_choices());
1386 const unsigned int X = numbers::invalid_unsigned_int;
1387 switch (this->kind)
1388 {
1390 {
1391 constexpr ::ndarray<unsigned int, 3, 12, 4> new_quad_lines_tet =
1392 {{// new line is (6,8)
1393 {{{{2, 3, 8, X}},
1394 {{0, 9, 5, X}},
1395 {{1, 6, 11, X}},
1396 {{4, 10, 7, X}},
1397 {{2, 12, 5, X}},
1398 {{1, 9, 12, X}},
1399 {{4, 8, 12, X}},
1400 {{6, 12, 10, X}},
1401 {{X, X, X, X}},
1402 {{X, X, X, X}},
1403 {{X, X, X, X}},
1404 {{X, X, X, X}}}},
1405 // new line is (5,7)
1406 {{{{2, 3, 8, X}},
1407 {{0, 9, 5, X}},
1408 {{1, 6, 11, X}},
1409 {{4, 10, 7, X}},
1410 {{0, 3, 12, X}},
1411 {{1, 12, 8, X}},
1412 {{4, 12, 9, X}},
1413 {{7, 11, 12, X}},
1414 {{X, X, X, X}},
1415 {{X, X, X, X}},
1416 {{X, X, X, X}},
1417 {{X, X, X, X}}}},
1418 // new line is (4,9)
1419 {{{{2, 3, 8, X}},
1420 {{0, 9, 5, X}},
1421 {{1, 6, 11, X}},
1422 {{4, 10, 7, X}},
1423 {{0, 12, 11, X}},
1424 {{2, 6, 12, X}},
1425 {{3, 12, 7, X}},
1426 {{5, 10, 12, X}},
1427 {{X, X, X, X}},
1428 {{X, X, X, X}},
1429 {{X, X, X, X}},
1430 {{X, X, X, X}}}}}};
1431 return new_quad_lines_tet[refienement_choice];
1432 }
1433
1435 {
1436 constexpr ::ndarray<unsigned int, 12, 4> new_quad_lines_hex = {
1437 {{{10, 28, 16, 24}},
1438 {{28, 14, 17, 25}},
1439 {{11, 29, 24, 20}},
1440 {{29, 15, 25, 21}},
1441 {{18, 26, 0, 28}},
1442 {{26, 22, 1, 29}},
1443 {{19, 27, 28, 4}},
1444 {{27, 23, 29, 5}},
1445 {{2, 24, 8, 26}},
1446 {{24, 6, 9, 27}},
1447 {{3, 25, 26, 12}},
1448 {{25, 7, 27, 13}}}};
1449 return new_quad_lines_hex;
1450 }
1451 default:
1453 }
1454
1455 return {};
1456}
1457
1458
1459
1460constexpr ::ndarray<unsigned int, 12, 4, 2>
1462 const unsigned int refienement_choice) const
1463{
1464 // AssertIndexRange(refienement_choice, n_isotropic_refinement_choices());
1465 const unsigned int X = numbers::invalid_unsigned_int;
1466 switch (this->kind)
1467 {
1469 {
1470 constexpr ::ndarray<unsigned int, 3, 12, 4, 2> table_tet = {
1471 {// new line is (6, 8)
1472 {{{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
1473 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
1474 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
1475 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
1476 {{{{4, 6}}, {{6, 8}}, {{8, 4}}, {{X, X}}}},
1477 {{{{6, 5}}, {{5, 8}}, {{8, 6}}, {{X, X}}}},
1478 {{{{8, 7}}, {{7, 6}}, {{6, 8}}, {{X, X}}}},
1479 {{{{9, 6}}, {{6, 8}}, {{8, 9}}, {{X, X}}}},
1480 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1481 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1482 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1483 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}},
1484 // new line is (5, 7)
1485 {{{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
1486 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
1487 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
1488 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
1489 {{{{5, 4}}, {{4, 7}}, {{7, 5}}, {{X, X}}}},
1490 {{{{6, 5}}, {{5, 7}}, {{7, 6}}, {{X, X}}}},
1491 {{{{8, 7}}, {{7, 5}}, {{5, 8}}, {{X, X}}}},
1492 {{{{7, 9}}, {{9, 5}}, {{5, 7}}, {{X, X}}}},
1493 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1494 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1495 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1496 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}},
1497 // new line is (4, 9)
1498 {{{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
1499 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
1500 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
1501 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
1502 {{{{5, 4}}, {{4, 9}}, {{9, 5}}, {{X, X}}}},
1503 {{{{4, 6}}, {{6, 9}}, {{9, 4}}, {{X, X}}}},
1504 {{{{7, 4}}, {{4, 9}}, {{9, 7}}, {{X, X}}}},
1505 {{{{4, 8}}, {{8, 9}}, {{9, 4}}, {{X, X}}}},
1506 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1507 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1508 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1509 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}}}};
1510 return table_tet[refienement_choice];
1511 }
1512
1514 {
1515 constexpr ::ndarray<unsigned int, 12, 4, 2> table_hex = {
1516 {{{{{10, 22}}, {{24, 26}}, {{10, 24}}, {{22, 26}}}},
1517 {{{{24, 26}}, {{11, 23}}, {{24, 11}}, {{26, 23}}}},
1518 {{{{22, 14}}, {{26, 25}}, {{22, 26}}, {{14, 25}}}},
1519 {{{{26, 25}}, {{23, 15}}, {{26, 23}}, {{25, 15}}}},
1520 {{{{8, 24}}, {{20, 26}}, {{8, 20}}, {{24, 26}}}},
1521 {{{{20, 26}}, {{12, 25}}, {{20, 12}}, {{26, 25}}}},
1522 {{{{24, 9}}, {{26, 21}}, {{24, 26}}, {{9, 21}}}},
1523 {{{{26, 21}}, {{25, 13}}, {{26, 25}}, {{21, 13}}}},
1524 {{{{16, 20}}, {{22, 26}}, {{16, 22}}, {{20, 26}}}},
1525 {{{{22, 26}}, {{17, 21}}, {{22, 17}}, {{26, 21}}}},
1526 {{{{20, 18}}, {{26, 23}}, {{20, 26}}, {{18, 23}}}},
1527 {{{{26, 23}}, {{21, 19}}, {{26, 21}}, {{23, 19}}}}}};
1528 return table_hex;
1529 }
1530 default:
1532 }
1533
1534 return {};
1535}
1536
1537
1538
1539constexpr ::ndarray<unsigned int, 8, 6>
1541 const unsigned int refienement_choice) const
1542{
1543 // AssertIndexRange(refienement_choice, n_isotropic_refinement_choices());
1544 const unsigned int X = numbers::invalid_unsigned_int;
1545 switch (this->kind)
1546 {
1548 {
1549 constexpr ::ndarray<unsigned int, 3, 8, 6> cell_quads_tet = {
1550 {// new line is (6, 8)
1551 {{{{8, 13, 16, 0, X, X}},
1552 {{9, 12, 1, 21, X, X}},
1553 {{10, 2, 17, 20, X, X}},
1554 {{3, 14, 18, 22, X, X}},
1555 {{11, 1, 4, 5, X, X}},
1556 {{15, 0, 4, 6, X, X}},
1557 {{19, 7, 6, 3, X, X}},
1558 {{23, 5, 2, 7, X, X}}}},
1559 // new line is (5, 7)
1560 {{
1561 {{8, 13, 16, 0, X, X}},
1562 {{9, 12, 1, 21, X, X}},
1563 {{10, 2, 17, 20, X, X}},
1564 {{3, 14, 18, 22, X, X}},
1565 {{11, 4, 0, 5, X, X}},
1566 {{15, 4, 1, 6, X, X}},
1567 {{19, 2, 5, 7, X, X}},
1568 {{23, 6, 7, 3, X, X}},
1569 }},
1570 // new line is (4, 9)
1571 {{
1572 {{8, 13, 16, 0, X, X}},
1573 {{9, 12, 1, 21, X, X}},
1574 {{10, 2, 17, 20, X, X}},
1575 {{3, 14, 18, 22, X, X}},
1576 {{11, 4, 5, 2, X, X}},
1577 {{15, 6, 7, 3, X, X}},
1578 {{19, 5, 0, 6, X, X}},
1579 {{23, 1, 4, 7, X, X}},
1580 }}}};
1581 return cell_quads_tet[refienement_choice];
1582 }
1583
1585 {
1586 constexpr ::ndarray<unsigned int, 8, 6> cell_quads_hex = {{
1587 {{12, 0, 20, 4, 28, 8}}, // bottom children
1588 {{0, 16, 22, 6, 29, 9}}, //
1589 {{13, 1, 4, 24, 30, 10}}, //
1590 {{1, 17, 6, 26, 31, 11}}, //
1591 {{14, 2, 21, 5, 8, 32}}, // top children
1592 {{2, 18, 23, 7, 9, 33}}, //
1593 {{15, 3, 5, 25, 10, 34}}, //
1594 {{3, 19, 7, 27, 11, 35}} //
1595 }};
1596 return cell_quads_hex;
1597 }
1598 default:
1600 }
1601
1602 return {};
1603}
1604
1605
1606
1607constexpr ::ndarray<unsigned int, 8, 4>
1609 const unsigned int refienement_choice) const
1610{
1611 // AssertIndexRange(refienement_choice, n_isotropic_refinement_choices());
1612
1613 switch (this->kind)
1614 {
1616 {
1617 constexpr ::ndarray<unsigned int, 3, 8, 4> cell_vertices_tet = {
1618 {// new line is (6,8)
1619 {{{{0, 4, 6, 7}},
1620 {{4, 1, 5, 8}},
1621 {{6, 5, 2, 9}},
1622 {{7, 8, 9, 3}},
1623 {{4, 5, 6, 8}},
1624 {{4, 7, 8, 6}},
1625 {{6, 9, 7, 8}},
1626 {{5, 8, 9, 6}}}},
1627 // new line is (5,7)
1628 {{{{0, 4, 6, 7}},
1629 {{4, 1, 5, 8}},
1630 {{6, 5, 2, 9}},
1631 {{7, 8, 9, 3}},
1632 {{4, 5, 6, 7}},
1633 {{4, 7, 8, 5}},
1634 {{6, 9, 7, 5}},
1635 {{5, 8, 9, 7}}}},
1636 // new line is (4,9)
1637 {{{{0, 4, 6, 7}},
1638 {{4, 1, 5, 8}},
1639 {{6, 5, 2, 9}},
1640 {{7, 8, 9, 3}},
1641 {{4, 5, 6, 9}},
1642 {{4, 7, 8, 9}},
1643 {{6, 9, 7, 4}},
1644 {{5, 8, 9, 4}}}}}};
1645 return cell_vertices_tet[refienement_choice];
1646 }
1647 default:
1649 }
1650
1651 return {};
1652}
1653
1654
1655
1656inline bool
1658{
1659 return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
1662}
1663
1664
1665
1666inline bool
1668{
1669 return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
1670 *this == ReferenceCells::Triangle ||
1672}
1673
1674
1675
1676inline unsigned int
1678{
1679 switch (this->kind)
1680 {
1682 return 0;
1684 return 1;
1687 return 2;
1692 return 3;
1693 default:
1695 }
1696
1698}
1699
1700
1701
1702template <int dim>
1705{
1706 Assert(dim == get_dimension(),
1707 ExcMessage("You can only call this function with arguments for "
1708 "which 'dim' equals the dimension of the space in which "
1709 "the current reference cell lives. You have dim=" +
1710 std::to_string(dim) + " but the reference cell is " +
1711 std::to_string(get_dimension()) + " dimensional."));
1712
1713 return Quadrature<dim>(std::vector<Point<dim>>({barycenter<dim>()}),
1714 std::vector<double>({volume()}));
1715}
1716
1717
1718
1719inline unsigned int
1721{
1722 switch (this->kind)
1723 {
1725 return 1;
1727 return 2;
1729 return 3;
1731 return 4;
1733 return 4;
1735 return 5;
1737 return 6;
1739 return 8;
1740 default:
1742 }
1743
1745}
1746
1747
1748
1749inline unsigned int
1751{
1752 switch (this->kind)
1753 {
1755 return 0;
1757 return 1;
1759 return 3;
1761 return 4;
1763 return 6;
1765 return 8;
1767 return 9;
1769 return 12;
1770 default:
1772 }
1773
1775}
1776
1777
1778
1779template <int dim>
1781ReferenceCell::vertex(const unsigned int v) const
1782{
1783 Assert(dim == get_dimension(),
1784 ExcMessage("You can only call this function with arguments for "
1785 "which 'dim' equals the dimension of the space in which "
1786 "the current reference cell lives. You have dim=" +
1787 std::to_string(dim) + " but the reference cell is " +
1788 std::to_string(get_dimension()) + " dimensional."));
1790
1791 switch (dim)
1792 {
1793 case 0:
1794 {
1795 if (*this == ReferenceCells::Vertex)
1796 return Point<dim>();
1797 break;
1798 }
1799 case 1:
1800 {
1801 static const Point<dim> vertices[2] = {
1802 Point<dim>(), // the origin
1803 Point<dim>::unit_vector(0) // unit point along x-axis
1804 };
1805 if (*this == ReferenceCells::Line)
1806 return vertices[v];
1807 break;
1808 }
1809 case 2:
1810 {
1811 switch (this->kind)
1812 {
1814 {
1815 static const Point<dim> vertices[3] = {
1816 Point<dim>(), // the origin
1817 Point<dim>::unit_vector(0), // unit point along x-axis
1818 Point<dim>::unit_vector(1) // unit point along y-axis
1819 };
1820 return vertices[v];
1821 }
1823 {
1824 static const Point<dim> vertices[4] = {
1825 // First the two points on the x-axis
1826 Point<dim>(),
1828 // Then these two points shifted in the y-direction
1831 return vertices[v];
1832 }
1833 }
1834 break;
1835 }
1836 case 3:
1837 {
1838 switch (this->kind)
1839 {
1841 {
1842 static const Point<dim> vertices[4] = {
1843 Point<dim>(), // the origin
1844 Point<dim>::unit_vector(0), // unit point along x-axis
1845 Point<dim>::unit_vector(1), // unit point along y-axis
1846 Point<dim>::unit_vector(2) // unit point along z-axis
1847 };
1848 return vertices[v];
1849 }
1851 {
1852 static const Point<dim> vertices[5] = {
1853 Point<dim>{-1.0, -1.0, 0.0},
1854 Point<dim>{+1.0, -1.0, 0.0},
1855 Point<dim>{-1.0, +1.0, 0.0},
1856 Point<dim>{+1.0, +1.0, 0.0},
1857 Point<dim>{+0.0, +0.0, 1.0}};
1858 return vertices[v];
1859 }
1861 {
1862 static const Point<dim> vertices[6] = {
1863 // First the three points on the triangular base of the
1864 // wedge:
1865 Point<dim>(),
1868 // And now everything shifted in the z-direction again
1872 return vertices[v];
1873 }
1875 {
1876 static const Point<dim> vertices[8] = {
1877 // First the two points on the x-axis
1878 Point<dim>(),
1880 // Then these two points shifted in the y-direction
1883 // And now all four points shifted in the z-direction
1890 return vertices[v];
1891 }
1892 }
1893 break;
1894 }
1895 default:
1897 }
1898
1900 return Point<dim>();
1901}
1902
1903
1904inline unsigned int
1906{
1907 switch (this->kind)
1908 {
1910 return 0;
1912 return 2;
1914 return 3;
1916 return 4;
1918 return 4;
1920 return 5;
1922 return 5;
1924 return 6;
1925 default:
1927 }
1928
1930}
1931
1932
1933
1940
1941
1942
1943inline unsigned int
1945{
1946 switch (this->kind)
1947 {
1949 return 1;
1951 return 1;
1953 return 1;
1955 return 1;
1957 return 3;
1959 return 1;
1961 return 1;
1963 return 1;
1964 default:
1966 }
1967
1969}
1970
1971
1972
1975 const unsigned int ref_choice) const
1976{
1978 switch (this->kind)
1979 {
1981 {
1983 isotropic_ref_choices = {{IsotropicRefinementChoice::cut_tet_68,
1986 return isotropic_ref_choices[ref_choice];
1987 }
1988 default:
1989 {
1991 }
1992 }
1993
1995}
1996
1997
1998
1999inline unsigned int
2001{
2002 switch (this->kind)
2003 {
2005 return 0;
2007 return 2;
2009 return 4;
2011 return 4;
2013 return 8;
2015 // We haven't yet decided how to refine pyramids. Update this when we
2016 // have
2017 return 0;
2019 return 8;
2021 return 8;
2022 default:
2024 }
2025
2027}
2028
2029
2030
2031template <int dim>
2032unsigned int
2034{
2035 Assert(dim == get_dimension(),
2036 ExcMessage("You can only call this function with arguments for "
2037 "which 'dim' equals the dimension of the space in which "
2038 "the current reference cell lives. You have dim=" +
2039 std::to_string(dim) + " but the reference cell is " +
2040 std::to_string(get_dimension()) + " dimensional."));
2041
2042 // Use GeometryInfo here to keep it the single source of truth
2043 if (this->is_hyper_cube())
2044 return GeometryInfo<dim>::n_children(ref_case);
2045 else
2046 return this->n_isotropic_children();
2047}
2048
2049
2050
2057
2058
2059
2066
2067
2068
2075
2076
2077
2078inline ReferenceCell
2079ReferenceCell::face_reference_cell(const unsigned int face_no) const
2080{
2081 AssertIndexRange(face_no, n_faces());
2082
2083 switch (this->kind)
2084 {
2091 return ReferenceCells::Line;
2095 if (face_no == 0)
2097 else
2100 if (face_no > 1)
2102 else
2106 default:
2108 }
2109
2111}
2112
2113
2114
2115inline constexpr types::geometric_orientation
2120
2121
2122
2123inline constexpr types::geometric_orientation
2125{
2126 // For a reversed line 'orientation' is false and neither flip nor rotate are
2127 // defined.
2129}
2130
2131
2132
2133inline unsigned int
2135 const unsigned int subface) const
2136{
2137 return child_cell_on_face(face,
2138 subface,
2140}
2141
2142
2143
2144inline unsigned int
2146 const unsigned int face,
2147 const unsigned int subface,
2148 const types::geometric_orientation combined_face_orientation) const
2149{
2150 AssertIndexRange(face, n_faces());
2152
2153 switch (this->kind)
2154 {
2157 {
2159 break;
2160 }
2162 {
2163 static constexpr ndarray<unsigned int, 3, 2> subcells = {
2164 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
2165
2166 Assert(combined_face_orientation ==
2168 combined_face_orientation ==
2171 return subcells[face][combined_face_orientation ==
2173 subface :
2174 1 - subface];
2175 }
2177 {
2178 const auto [face_orientation, face_rotation, face_flip] =
2179 internal::split_face_orientation(combined_face_orientation);
2180
2183 face,
2184 subface,
2185 face_orientation,
2186 face_flip,
2187 face_rotation);
2188 }
2192 {
2194 break;
2195 }
2197 {
2198 const auto [face_orientation, face_rotation, face_flip] =
2199 internal::split_face_orientation(combined_face_orientation);
2200
2203 face,
2204 subface,
2205 face_orientation,
2206 face_flip,
2207 face_rotation);
2208 }
2209 default:
2211 }
2212
2214}
2215
2216
2217
2218inline std::array<unsigned int, 2>
2220 const unsigned int vertex) const
2221{
2223 // Work around a GCC warning at higher optimization levels by making all of
2224 // these tables the same size
2225 constexpr unsigned int X = numbers::invalid_unsigned_int;
2226
2227 switch (this->kind)
2228 {
2232 break;
2234 {
2235 static constexpr ndarray<unsigned int, 6, 2> table = {
2236 {{{0, 0}}, {{0, 1}}, {{1, 1}}, {{X, X}}, {{X, X}}, {{X, X}}}};
2237
2238 return table[vertex];
2239 }
2241 {
2243 vertex);
2244 }
2246 {
2247 static constexpr ndarray<unsigned int, 6, 2> table = {
2248 {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 2}}, {{X, X}}, {{X, X}}}};
2249
2250 return table[vertex];
2251 }
2253 {
2254 static constexpr ndarray<unsigned int, 6, 2> table = {
2255 {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{0, 3}}, {{1, 2}}, {{X, X}}}};
2256
2257 return table[vertex];
2258 }
2260 {
2261 static constexpr ndarray<unsigned int, 6, 2> table = {
2262 {{{0, 1}}, {{0, 0}}, {{0, 2}}, {{1, 0}}, {{1, 1}}, {{1, 2}}}};
2263
2264 return table[vertex];
2265 }
2267 {
2269 vertex);
2270 }
2271 default:
2273 }
2274
2275 return {};
2276}
2277
2278
2279
2280inline std::array<unsigned int, 2>
2282 const unsigned int line) const
2283{
2284 AssertIndexRange(line, n_lines());
2285
2286 switch (this->kind)
2287 {
2292 {
2294 break;
2295 }
2297 {
2298 static const std::array<unsigned int, 2> table[6] = {
2299 {{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}};
2300
2301 return table[line];
2302 }
2304 {
2305 static const std::array<unsigned int, 2> table[8] = {{{0, 0}},
2306 {{0, 1}},
2307 {{0, 2}},
2308 {{0, 3}},
2309 {{1, 2}},
2310 {{2, 1}},
2311 {{1, 1}},
2312 {{2, 2}}};
2313
2314 return table[line];
2315 }
2317 {
2318 static const std::array<unsigned int, 2> table[9] = {{{0, 0}},
2319 {{0, 2}},
2320 {{0, 1}},
2321 {{1, 0}},
2322 {{1, 1}},
2323 {{1, 2}},
2324 {{2, 0}},
2325 {{2, 1}},
2326 {{3, 1}}};
2327
2328 return table[line];
2329 }
2331 {
2333 }
2334 default:
2336 }
2337
2338 return {};
2339}
2340
2341
2342inline unsigned int
2344 const unsigned int vertex) const
2345{
2347 AssertIndexRange(line, n_lines());
2348
2349 switch (this->kind)
2350 {
2353 return vertex;
2355 {
2356 static constexpr ndarray<unsigned int, 3, 2> table = {
2357 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
2358 return table[line][vertex];
2359 }
2361 {
2362 static constexpr ndarray<unsigned int, 4, 2> table = {
2363 {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
2364 return table[line][vertex];
2365 }
2367 {
2368 static constexpr ndarray<unsigned int, 6, 2> table = {
2369 {{{0, 1}}, {{1, 2}}, {{2, 0}}, {{0, 3}}, {{1, 3}}, {{2, 3}}}};
2370 return table[line][vertex];
2371 }
2373 {
2374 static constexpr ndarray<unsigned int, 8, 2> table = {{{{0, 2}},
2375 {{1, 3}},
2376 {{0, 1}},
2377 {{2, 3}},
2378 {{4, 0}},
2379 {{1, 4}},
2380 {{2, 4}},
2381 {{4, 3}}}};
2382 return table[line][vertex];
2383 }
2385 {
2386 static constexpr ndarray<unsigned int, 9, 2> table = {{{{1, 0}},
2387 {{2, 1}},
2388 {{0, 2}},
2389 {{3, 4}},
2390 {{4, 5}},
2391 {{5, 3}},
2392 {{0, 3}},
2393 {{1, 4}},
2394 {{2, 5}}}};
2395 return table[line][vertex];
2396 }
2398 {
2399 // first four lines comprise the bottom face, next four are the top,
2400 // and the last four are 'bottom to top'
2401 static constexpr ndarray<unsigned int, 12, 2> table = {{{{0, 2}},
2402 {{1, 3}},
2403 {{0, 1}},
2404 {{2, 3}},
2405 {{4, 6}},
2406 {{5, 7}},
2407 {{4, 5}},
2408 {{6, 7}},
2409 {{0, 4}},
2410 {{1, 5}},
2411 {{2, 6}},
2412 {{3, 7}}}};
2413 return table[line][vertex];
2414 }
2415
2416 default:
2418 }
2419
2421}
2422
2423
2424inline unsigned int
2426 const unsigned int face,
2427 const unsigned int line,
2428 const types::geometric_orientation combined_face_orientation) const
2429{
2430 AssertIndexRange(face, n_faces());
2432
2433 static constexpr unsigned int X = numbers::invalid_unsigned_int;
2434
2435 switch (this->kind)
2436 {
2438 {
2440 break;
2441 }
2443 {
2444 const auto [face_orientation, face_rotation, face_flip] =
2445 internal::split_face_orientation(combined_face_orientation);
2446
2448 face, line, face_orientation, face_flip, face_rotation);
2449 }
2451 {
2452 return face;
2453 }
2455 {
2456 const auto [face_orientation, face_rotation, face_flip] =
2457 internal::split_face_orientation(combined_face_orientation);
2458
2460 face, line, face_orientation, face_flip, face_rotation);
2461 }
2463 {
2464 static constexpr ndarray<unsigned int, 4, 3> table = {
2465 {{{0, 1, 2}}, {{0, 3, 4}}, {{2, 5, 3}}, {{1, 4, 5}}}};
2466
2467 return table[face][standard_to_real_face_line(
2468 line, face, combined_face_orientation)];
2469 }
2471 {
2472 static constexpr ndarray<unsigned int, 5, 4> table = {
2473 {{{0, 1, 2, 3}},
2474 {{0, 6, 4, X}},
2475 {{1, 5, 7, X}},
2476 {{2, 4, 5, X}},
2477 {{3, 7, 6, X}}}};
2478
2479 return table[face][standard_to_real_face_line(
2480 line, face, combined_face_orientation)];
2481 }
2483 {
2484 static constexpr ndarray<unsigned int, 5, 4> table = {
2485 {{{0, 2, 1, X}},
2486 {{3, 4, 5, X}},
2487 {{6, 7, 0, 3}},
2488 {{7, 8, 1, 4}},
2489 {{8, 6, 5, 2}}}};
2490
2491 return table[face][standard_to_real_face_line(
2492 line, face, combined_face_orientation)];
2493 }
2495 {
2496 const auto [face_orientation, face_rotation, face_flip] =
2497 internal::split_face_orientation(combined_face_orientation);
2498
2500 face, line, face_orientation, face_flip, face_rotation);
2501 }
2502 default:
2504 }
2505
2507}
2508
2509
2510
2511inline unsigned int
2513 const unsigned int face,
2514 const unsigned int vertex,
2515 const types::geometric_orientation combined_face_orientation) const
2516{
2517 AssertIndexRange(face, n_faces());
2519 // TODO: once the default orientation is switched to 0 then we can remove this
2520 // special case for 1D.
2521 if (get_dimension() == 1)
2522 Assert(combined_face_orientation == numbers::default_geometric_orientation,
2523 ExcMessage("In 1D, all faces must have the default orientation."));
2524 else
2525 AssertIndexRange(combined_face_orientation, n_face_orientations(face));
2526
2527 switch (this->kind)
2528 {
2530 {
2532 break;
2533 }
2535 {
2536 const auto [face_orientation, face_rotation, face_flip] =
2537 internal::split_face_orientation(combined_face_orientation);
2538
2540 face, vertex, face_orientation, face_flip, face_rotation);
2541 }
2543 {
2544 static constexpr ndarray<unsigned int, 3, 2> table = {
2545 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
2546
2547 return table[face][combined_face_orientation ==
2549 vertex :
2550 (1 - vertex)];
2551 }
2553 {
2554 const auto [face_orientation, face_rotation, face_flip] =
2555 internal::split_face_orientation(combined_face_orientation);
2556
2558 face, vertex, face_orientation, face_flip, face_rotation);
2559 }
2561 {
2562 static constexpr ndarray<unsigned int, 4, 3> table = {
2563 {{{0, 1, 2}}, {{1, 0, 3}}, {{0, 2, 3}}, {{2, 1, 3}}}};
2564
2565 return table[face][standard_to_real_face_vertex(
2566 vertex, face, combined_face_orientation)];
2567 }
2569 {
2570 constexpr auto X = numbers::invalid_unsigned_int;
2571 static constexpr ndarray<unsigned int, 5, 4> table = {
2572 {{{0, 1, 2, 3}},
2573 {{0, 2, 4, X}},
2574 {{3, 1, 4, X}},
2575 {{1, 0, 4, X}},
2576 {{2, 3, 4, X}}}};
2577
2578 return table[face][standard_to_real_face_vertex(
2579 vertex, face, combined_face_orientation)];
2580 }
2582 {
2583 constexpr auto X = numbers::invalid_unsigned_int;
2584 static constexpr ndarray<unsigned int, 6, 4> table = {
2585 {{{1, 0, 2, X}},
2586 {{3, 4, 5, X}},
2587 {{0, 1, 3, 4}},
2588 {{1, 2, 4, 5}},
2589 {{2, 0, 5, 3}}}};
2590
2591 return table[face][standard_to_real_face_vertex(
2592 vertex, face, combined_face_orientation)];
2593 }
2595 {
2596 const auto [face_orientation, face_rotation, face_flip] =
2597 internal::split_face_orientation(combined_face_orientation);
2598
2600 face, vertex, face_orientation, face_flip, face_rotation);
2601 }
2602 default:
2604 }
2605
2607}
2608
2609
2610
2611template <int dim>
2614 const unsigned int vertex) const
2615{
2616 Assert(dim == get_dimension(),
2617 ExcMessage("You can only call this function with arguments for "
2618 "which 'dim' equals the dimension of the space in which "
2619 "the current reference cell lives. You have dim=" +
2620 std::to_string(dim) + " but the reference cell is " +
2621 std::to_string(get_dimension()) + " dimensional."));
2622 return this->template vertex<dim>(face_to_cell_vertices(
2624}
2625
2626
2627
2628inline unsigned int
2630 const unsigned int vertex,
2631 const unsigned int face,
2632 const types::geometric_orientation face_orientation) const
2633{
2634 AssertIndexRange(face, n_faces());
2636
2637 switch (this->kind)
2638 {
2641 break;
2644 ExcMessage(
2645 "In 1D, all faces must have the default orientation."));
2646 return vertex;
2649 return line_vertex_permutations[face_orientation][vertex];
2651 return triangle_vertex_permutations[face_orientation][vertex];
2653 // face 0 is a quadrilateral
2654 if (face == 0)
2655 return quadrilateral_vertex_permutations[face_orientation][vertex];
2656 else
2657 return triangle_vertex_permutations[face_orientation][vertex];
2659 // faces 0 and 1 are triangles
2660 if (face > 1)
2661 return quadrilateral_vertex_permutations[face_orientation][vertex];
2662 else
2663 return triangle_vertex_permutations[face_orientation][vertex];
2665 return quadrilateral_vertex_permutations[face_orientation][vertex];
2666 default:
2668 }
2669
2672}
2673
2674
2675
2676inline unsigned int
2678 const unsigned int line,
2679 const unsigned int face,
2680 const types::geometric_orientation combined_face_orientation) const
2681{
2682 AssertIndexRange(face, n_faces());
2684
2685 static constexpr ndarray<unsigned int, 6, 3> triangle_table = {{{{2, 1, 0}},
2686 {{0, 1, 2}},
2687 {{1, 0, 2}},
2688 {{2, 0, 1}},
2689 {{0, 2, 1}},
2690 {{1, 2, 0}}}};
2691
2692
2693 switch (this->kind)
2694 {
2700 break;
2702 return triangle_table[combined_face_orientation][line];
2704 if (face == 0) // The quadrilateral face
2705 {
2706 const auto [face_orientation, face_rotation, face_flip] =
2707 internal::split_face_orientation(combined_face_orientation);
2708
2710 face_orientation,
2711 face_flip,
2712 face_rotation);
2713 }
2714 else // One of the triangular faces
2715 {
2716 return triangle_table[combined_face_orientation][line];
2717 }
2719 if (face > 1) // One of the quadrilateral faces
2720 {
2721 const auto [face_orientation, face_rotation, face_flip] =
2722 internal::split_face_orientation(combined_face_orientation);
2723
2725 face_orientation,
2726 face_flip,
2727 face_rotation);
2728 }
2729 else // One of the triangular faces
2730 return triangle_table[combined_face_orientation][line];
2732 {
2733 static constexpr ndarray<unsigned int, 8, 4> table = {
2734 {{{2, 3, 0, 1}},
2735 {{0, 1, 2, 3}},
2736 {{0, 1, 3, 2}},
2737 {{3, 2, 0, 1}},
2738 {{3, 2, 1, 0}},
2739 {{1, 0, 3, 2}},
2740 {{1, 0, 2, 3}},
2741 {{2, 3, 1, 0}}}};
2742 return table[combined_face_orientation][line];
2743 }
2744 default:
2746 }
2747
2749}
2750
2751
2752
2753namespace ReferenceCells
2754{
2755 template <int dim>
2756 inline constexpr const ReferenceCell &
2758 {
2759 switch (dim)
2760 {
2761 case 0:
2763 case 1:
2764 return ReferenceCells::Line;
2765 case 2:
2767 case 3:
2769 default:
2771 }
2773 }
2774
2775
2776
2777 template <int dim>
2778 inline constexpr const ReferenceCell &
2780 {
2781 switch (dim)
2782 {
2783 case 0:
2785 case 1:
2786 return ReferenceCells::Line;
2787 case 2:
2789 case 3:
2791 default:
2793 }
2795 }
2796
2797
2798
2799 template <int structdim>
2800 inline constexpr unsigned int
2805
2806
2807
2808 template <int structdim>
2809 inline constexpr unsigned int
2814
2815
2816
2817 template <int structdim>
2818 inline constexpr unsigned int
2823} // namespace ReferenceCells
2824
2825
2826inline ReferenceCell
2827ReferenceCell::n_vertices_to_type(const int dim, const unsigned int n_vertices)
2828{
2829 AssertIndexRange(dim, 4);
2831
2832 const auto X = ReferenceCells::Invalid;
2833 static const ndarray<ReferenceCell, 4, 9> table = {
2834 {// dim 0
2835 {{X, ReferenceCells::Vertex, X, X, X, X, X, X, X}},
2836 // dim 1
2837 {{X, X, ReferenceCells::Line, X, X, X, X, X, X}},
2838 // dim 2
2839 {{X,
2840 X,
2841 X,
2844 X,
2845 X,
2846 X,
2847 X}},
2848 // dim 3
2849 {{X,
2850 X,
2851 X,
2852 X,
2856 X,
2859 ExcMessage("The combination of dim = " + std::to_string(dim) +
2860 " and n_vertices = " + std::to_string(n_vertices) +
2861 " does not correspond to a known reference cell type."));
2862 return table[dim][n_vertices];
2863}
2864
2865
2866
2867template <int dim>
2868inline double
2870 const unsigned int i) const
2871{
2872 Assert(dim == get_dimension(),
2873 ExcMessage("You can only call this function with arguments for "
2874 "which 'dim' equals the dimension of the space in which "
2875 "the current reference cell lives. You have dim=" +
2876 std::to_string(dim) + " but the reference cell is " +
2877 std::to_string(get_dimension()) + " dimensional."));
2879 switch (this->kind)
2880 {
2886 // see also BarycentricPolynomials<2>::compute_value
2888 {
2889 switch (i)
2890 {
2891 case 0:
2892 return 1.0 - xi[std::min(0, dim - 1)] -
2893 xi[std::min(1, dim - 1)];
2894 case 1:
2895 return xi[std::min(0, dim - 1)];
2896 case 2:
2897 return xi[std::min(1, dim - 1)];
2898 default:
2900 }
2901 }
2902 // see also BarycentricPolynomials<3>::compute_value
2904 {
2905 switch (i)
2906 {
2907 case 0:
2908 return 1.0 - xi[std::min(0, dim - 1)] -
2909 xi[std::min(1, dim - 1)] - xi[std::min(2, dim - 1)];
2910 case 1:
2911 return xi[std::min(0, dim - 1)];
2912 case 2:
2913 return xi[std::min(1, dim - 1)];
2914 case 3:
2915 return xi[std::min(2, dim - 1)];
2916 default:
2918 }
2919 }
2920 // see also ScalarLagrangePolynomialPyramid::compute_value()
2922 {
2923 const double Q14 = 0.25;
2924
2925 const double r = xi[std::min(0, dim - 1)];
2926 const double s = xi[std::min(1, dim - 1)];
2927 const double t = xi[std::min(2, dim - 1)];
2928
2929 const double ratio =
2930 (std::fabs(t - 1.0) > 1.0e-14 ? (r * s * t) / (1.0 - t) : 0.0);
2931
2932 if (i == 0)
2933 return Q14 * ((1.0 - r) * (1.0 - s) - t + ratio);
2934 if (i == 1)
2935 return Q14 * ((1.0 + r) * (1.0 - s) - t - ratio);
2936 if (i == 2)
2937 return Q14 * ((1.0 - r) * (1.0 + s) - t - ratio);
2938 if (i == 3)
2939 return Q14 * ((1.0 + r) * (1.0 + s) - t + ratio);
2940 else
2941 return t;
2942 }
2943 // see also ScalarLagrangePolynomialWedge::compute_value()
2946 .d_linear_shape_function<2>(Point<2>(xi[std::min(0, dim - 1)],
2947 xi[std::min(1, dim - 1)]),
2948 i % 3) *
2950 .d_linear_shape_function<1>(Point<1>(xi[std::min(2, dim - 1)]),
2951 i / 3);
2952 default:
2954 }
2955
2956 return 0.0;
2957}
2958
2959
2960
2961template <int dim>
2962inline Tensor<1, dim>
2964 const unsigned int i) const
2965{
2966 Assert(dim == get_dimension(),
2967 ExcMessage("You can only call this function with arguments for "
2968 "which 'dim' equals the dimension of the space in which "
2969 "the current reference cell lives. You have dim=" +
2970 std::to_string(dim) + " but the reference cell is " +
2971 std::to_string(get_dimension()) + " dimensional."));
2972 switch (this->kind)
2973 {
2979 // see also BarycentricPolynomials<2>::compute_grad()
2981 switch (i)
2982 {
2983 case 0:
2984 return Point<dim>(-1.0, -1.0);
2985 case 1:
2986 return Point<dim>(+1.0, +0.0);
2987 case 2:
2988 return Point<dim>(+0.0, +1.0);
2989 default:
2991 }
2992 default:
2994 }
2995
2996 return Point<dim>(+0.0, +0.0, +0.0);
2997}
2998
2999
3000
3001inline double
3003{
3004 switch (this->kind)
3005 {
3007 return 0;
3009 return 1;
3011 return 1. / 2.;
3013 return 1;
3015 return 1. / 6.;
3017 return 4. / 3.;
3019 return 1. / 2.;
3021 return 1;
3022 default:
3024 }
3025
3026 return 0.0;
3027}
3028
3029
3030
3031template <int dim>
3032inline Point<dim>
3034{
3035 Assert(dim == get_dimension(),
3036 ExcMessage("You can only call this function with arguments for "
3037 "which 'dim' equals the dimension of the space in which "
3038 "the current reference cell lives. You have dim=" +
3039 std::to_string(dim) + " but the reference cell is " +
3040 std::to_string(get_dimension()) + " dimensional."));
3041
3042 switch (this->kind)
3043 {
3045 return Point<dim>();
3047 return Point<dim>(1. / 2.);
3049 return Point<dim>(1. / 3., 1. / 3.);
3051 return Point<dim>(1. / 2., 1. / 2.);
3053 return Point<dim>(1. / 4., 1. / 4., 1. / 4.);
3055 return Point<dim>(0, 0, 1. / 4.);
3057 return Point<dim>(1. / 3, 1. / 3, 1. / 2.);
3059 return Point<dim>(1. / 2., 1. / 2., 1. / 2.);
3060 default:
3062 }
3063
3064 return Point<dim>();
3065}
3066
3067
3068
3069template <int dim>
3070inline bool
3071ReferenceCell::contains_point(const Point<dim> &p, const double tolerance) const
3072{
3073 Assert(dim == get_dimension(),
3074 ExcMessage("You can only call this function with arguments for "
3075 "which 'dim' equals the dimension of the space in which "
3076 "the current reference cell lives. You have dim=" +
3077 std::to_string(dim) + " but the reference cell is " +
3078 std::to_string(get_dimension()) + " dimensional."));
3079
3080 // Introduce abbreviations to silence compiler warnings about invalid
3081 // array accesses (that can't happen, of course, but the compiler
3082 // doesn't know that).
3083 constexpr unsigned int x_coordinate = 0;
3084 constexpr unsigned int y_coordinate = (dim >= 2 ? 1 : 0);
3085 constexpr unsigned int z_coordinate = (dim >= 3 ? 2 : 0);
3086
3087 switch (this->kind)
3088 {
3090 {
3091 // Vertices are special cases in that they do not actually
3092 // have coordinates. Error out if this function is called
3093 // with a vertex:
3094 Assert(false,
3095 ExcMessage("Vertices are zero-dimensional objects and "
3096 "as a consequence have no coordinates. You "
3097 "cannot meaningfully ask whether a point is "
3098 "inside a vertex (within a certain tolerance) "
3099 "without coordinate values."));
3100 return false;
3101 }
3105 {
3106 for (unsigned int d = 0; d < dim; ++d)
3107 if ((p[d] < -tolerance) || (p[d] > 1 + tolerance))
3108 return false;
3109 return true;
3110 }
3113 {
3114 // First make sure that we are in the first quadrant or octant
3115 for (unsigned int d = 0; d < dim; ++d)
3116 if (p[d] < -tolerance)
3117 return false;
3118
3119 // Now we also need to make sure that we are below the diagonal line
3120 // or plane that delineates the simplex. This diagonal is given by
3121 // sum(p[d])<=1, and a diagonal a distance eps away is given by
3122 // sum(p[d])<=1+eps*sqrt(d). (For example, the point at (1,1) is a
3123 // distance of 1/sqrt(2) away from the diagonal. That is, its
3124 // sum satisfies
3125 // sum(p[d]) = 2 <= 1 + (1/sqrt(2)) * sqrt(2)
3126 // in other words, it satisfies the predicate with eps=1/sqrt(2).)
3127 double sum = 0;
3128 for (unsigned int d = 0; d < dim; ++d)
3129 sum += p[d];
3130 return (sum <= 1 + tolerance * std::sqrt(1. * dim));
3131 }
3133 {
3134 // A pyramid only lives in the upper half-space:
3135 if (p[z_coordinate] < -tolerance)
3136 return false;
3137
3138 // It also only lives in the space below z=1:
3139 if (p[z_coordinate] > 1 + tolerance)
3140 return false;
3141
3142 // Within what's left of the space, a pyramid is a cone that tapers
3143 // towards the top. First compute the distance of the point to the
3144 // axis in the max norm (this is the right norm because the vertices
3145 // of the pyramid are at points +/-1, +/-1):
3146 const double distance_from_axis =
3147 std::max(std::fabs(p[x_coordinate]), std::fabs(p[y_coordinate]));
3148
3149 // We are inside the pyramid if the distance from the axis is less
3150 // than (1-z)
3151 return (distance_from_axis <= 1 + tolerance - p[z_coordinate]);
3152 }
3154 {
3155 // The wedge we use is a triangle extruded into the third
3156 // dimension by one unit. So we can use the same logic as for
3157 // triangles above (i.e., for the simplex above, using dim==2)
3158 // and then check the third dimension separately.
3159
3160 if ((p[x_coordinate] < -tolerance) || (p[y_coordinate] < -tolerance))
3161 return false;
3162
3163 const double sum = p[x_coordinate] + p[y_coordinate];
3164 if (sum > 1 + tolerance * std::sqrt(2.0))
3165 return false;
3166
3167 if (p[z_coordinate] < -tolerance)
3168 return false;
3169 if (p[z_coordinate] > 1 + tolerance)
3170 return false;
3171
3172 return true;
3173 }
3174 default:
3176 }
3177
3178 return false;
3179}
3180
3181
3182
3183template <int dim>
3184inline Tensor<1, dim>
3186 const unsigned int i) const
3187{
3188 return face_tangent_vector<dim>(face_no, i);
3189}
3190
3191
3192
3193template <int dim>
3194inline Tensor<1, dim>
3195ReferenceCell::face_tangent_vector(const unsigned int face_no,
3196 const unsigned int i) const
3197{
3198 Assert(dim == get_dimension(),
3199 ExcMessage("You can only call this function with arguments for "
3200 "which 'dim' equals the dimension of the space in which "
3201 "the current reference cell lives. You have dim=" +
3202 std::to_string(dim) + " but the reference cell is " +
3203 std::to_string(get_dimension()) + " dimensional."));
3204 AssertIndexRange(i, dim - 1);
3205
3206 switch (this->kind)
3207 {
3215 {
3216 AssertIndexRange(face_no, 3);
3217 static const std::array<Tensor<1, dim>, 3> table = {
3218 {Point<dim>(1, 0),
3219 Point<dim>(-std::sqrt(0.5), +std::sqrt(0.5)),
3220 Point<dim>(0, -1)}};
3221
3222 return table[face_no];
3223 }
3225 {
3226 AssertIndexRange(face_no, 4);
3227 static const ndarray<Tensor<1, dim>, 4, 2> table = {
3228 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
3229 {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
3230 {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
3231 {{Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
3232 +std::pow(1.0 / 3.0, 1.0 / 4.0),
3233 0),
3234 Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
3235 0,
3236 +std::pow(1.0 / 3.0, 1.0 / 4.0))}}}};
3237
3238 return table[face_no][i];
3239 }
3241 {
3242 AssertIndexRange(face_no, 5);
3243 static const ndarray<Tensor<1, dim>, 5, 2> table = {
3244 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
3245 {{Point<dim>(+1.0 / sqrt(2.0), 0, +1.0 / sqrt(2.0)),
3246 Point<dim>(0, 1, 0)}},
3247 {{Point<dim>(+1.0 / sqrt(2.0), 0, -1.0 / sqrt(2.0)),
3248 Point<dim>(0, 1, 0)}},
3249 {{Point<dim>(1, 0, 0),
3250 Point<dim>(0, +1.0 / sqrt(2.0), +1.0 / sqrt(2.0))}},
3251 {{Point<dim>(1, 0, 0),
3252 Point<dim>(0, +1.0 / sqrt(2.0), -1.0 / sqrt(2.0))}}}};
3253
3254 return table[face_no][i];
3255 }
3257 {
3258 AssertIndexRange(face_no, 5);
3259 static const ndarray<Tensor<1, dim>, 5, 2> table = {
3260 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
3261 {{Point<dim>(1, 0, 0), Point<dim>(0, 1, 0)}},
3262 {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
3263 {{Point<dim>(-1 / std::sqrt(2.0), +1 / std::sqrt(2.0), 0),
3264 Point<dim>(0, 0, 1)}},
3265 {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}}}};
3266
3267 return table[face_no][i];
3268 }
3269 default:
3271 }
3272
3273
3274 return {};
3275}
3276
3277
3278
3279template <int dim>
3280inline Tensor<1, dim>
3281ReferenceCell::unit_normal_vectors(const unsigned int face_no) const
3282{
3283 return face_normal_vector<dim>(face_no);
3284}
3285
3286
3287
3288template <int dim>
3289inline Tensor<1, dim>
3290ReferenceCell::face_normal_vector(const unsigned int face_no) const
3291{
3292 Assert(dim == get_dimension(),
3293 ExcMessage("You can only call this function with arguments for "
3294 "which 'dim' equals the dimension of the space in which "
3295 "the current reference cell lives. You have dim=" +
3296 std::to_string(dim) + " but the reference cell is " +
3297 std::to_string(get_dimension()) + " dimensional."));
3298
3299 if (is_hyper_cube())
3300 {
3303 }
3304 else if (dim == 2)
3305 {
3307
3308 // Return the rotated vector
3309 return cross_product_2d(face_tangent_vector<dim>(face_no, 0));
3310 }
3311 else if (dim == 3)
3312 {
3313 return cross_product_3d(face_tangent_vector<dim>(face_no, 0),
3314 face_tangent_vector<dim>(face_no, 1));
3315 }
3316
3318
3319 return {};
3320}
3321
3322
3323
3324inline unsigned int
3325ReferenceCell::n_face_orientations(const unsigned int face_no) const
3326{
3327 AssertIndexRange(face_no, n_faces());
3328 if (get_dimension() == 1)
3329 return 1;
3330 if (get_dimension() == 2)
3331 return 2;
3333 return 8;
3334 else if (face_reference_cell(face_no) == ReferenceCells::Triangle)
3335 return 6;
3336
3339}
3340
3341
3342
3345 const unsigned int line,
3346 const unsigned int face,
3347 const types::geometric_orientation combined_face_orientation,
3348 const types::geometric_orientation line_orientation) const
3349{
3351 face,
3352 combined_face_orientation,
3353 line_orientation);
3354}
3355
3356
3359 const unsigned int face_line_no,
3360 const unsigned int face_no,
3361 const types::geometric_orientation combined_face_orientation,
3362 const types::geometric_orientation line_orientation) const
3363{
3364 constexpr auto D = numbers::default_geometric_orientation;
3365 constexpr auto R = numbers::reverse_line_orientation;
3366 if (*this == ReferenceCells::Hexahedron)
3367 {
3368 static constexpr ::ndarray<types::geometric_orientation, 2, 8>
3369 table{{{{D, D, R, D, R, R, D, R}}, {{D, D, D, R, R, R, R, D}}}};
3370 // We use face_line_no / 2 here since lines i and i + 1 are parallel and,
3371 // on a given face, have the same relative orientations.
3372 const bool match =
3373 line_orientation == table[face_line_no / 2][combined_face_orientation];
3374
3377 }
3378 else if (*this == ReferenceCells::Tetrahedron)
3379 {
3380 static constexpr unsigned int X = numbers::invalid_unsigned_int;
3381 static constexpr ::ndarray<unsigned int, 4, 3> combined_lines{
3382 {{{0, 0, 0}}, {{X, 0, 1}}, {{X, 0, X}}, {{X, X, X}}}};
3383
3384 const auto combined_line = combined_lines[face_no][face_line_no];
3385
3386 Assert(combined_line != X,
3387 ExcMessage(
3388 "This function can only be called for following face-line "
3389 "combinations: (0,0), (0,1), (0,2), (1,1), (1,2), (2,1),"));
3390
3391 static constexpr ::ndarray<types::geometric_orientation, 2, 6>
3392 table{{{{R, D, R, D, R, D}}, {{D, R, D, R, D, R}}}};
3393
3394 const bool match =
3395 line_orientation == table[combined_line][combined_face_orientation];
3396
3399 }
3400 else
3401 // TODO: This might actually be wrong for some of the other
3402 // kinds of objects. We should check this
3404}
3405
3406
3407
3408namespace internal
3409{
3410 template <typename T>
3412 {
3413 public:
3427
3431 virtual ~NoPermutation() noexcept override = default;
3432
3436 virtual void
3437 print_info(std::ostream &out) const override
3438 {
3439 out << '[';
3440
3441 const unsigned int n_vertices = entity_type.n_vertices();
3442
3443 for (unsigned int i = 0; i < n_vertices; ++i)
3444 {
3445 out << vertices_0[i];
3446 if (i + 1 != n_vertices)
3447 out << ',';
3448 }
3449
3450 out << "] is not a valid permutation of [";
3451
3452 for (unsigned int i = 0; i < n_vertices; ++i)
3453 {
3454 out << vertices_1[i];
3455 if (i + 1 != n_vertices)
3456 out << ',';
3457 }
3458
3459 out << "]." << std::endl;
3460 }
3461
3466
3471
3476 };
3477
3489 << "The reference-cell type used on this cell (" << arg1.to_string()
3490 << ") does not match the reference-cell type of the finite element "
3491 << "associated with this cell (" << arg2.to_string() << "). "
3492 << "Did you accidentally use simplex elements on hypercube meshes "
3493 << "(or the other way around), or are you using a mixed mesh and "
3494 << "assigned a simplex element to a hypercube cell (or the other "
3495 << "way around) via the active_fe_index?");
3496} // namespace internal
3497
3498
3499
3500template <typename T, std::size_t N>
3502ReferenceCell::compute_orientation(const std::array<T, N> &vertices_0,
3503 const std::array<T, N> &vertices_1) const
3504{
3505 Assert(N >= n_vertices(),
3506 ExcMessage("The number of array elements must be equal to or "
3507 "greater than the number of vertices of the cell "
3508 "referenced by this object."));
3509
3510 // Call the non-deprecated function, taking care of calling it only with
3511 // those array elements that we actually care about (see the note
3512 // in the documentation about the arguments potentially being
3513 // larger arrays than necessary).
3515 make_array_view(vertices_0.begin(), vertices_0.begin() + n_vertices()),
3516 make_array_view(vertices_1.begin(), vertices_1.begin() + n_vertices()));
3517}
3518
3519
3520
3521template <typename T>
3524 const ArrayView<const T> &vertices_0,
3525 const ArrayView<const T> &vertices_1) const
3526{
3527 Assert(vertices_0.size() == n_vertices(),
3528 ExcMessage("The number of array elements must be equal to "
3529 "the number of vertices of the cell "
3530 "referenced by this object."));
3531 Assert(vertices_1.size() == n_vertices(),
3532 ExcMessage("The number of array elements must be equal to "
3533 "the number of vertices of the cell "
3534 "referenced by this object."));
3535
3536 auto compute_orientation = [&](const auto &table) {
3537 for (types::geometric_orientation o = 0; o < table.size(); ++o)
3538 {
3539 bool match = true;
3540 for (unsigned int j = 0; j < table[o].size(); ++j)
3541 match = (match && vertices_0[j] == vertices_1[table[o][j]]);
3542
3543 if (match)
3544 return o;
3545 }
3546
3547 // Do not use `this` in Assert because nvcc when using C++20 assumes that
3548 // `this` is an integer and we get the following error: invalid type
3549 // argument of unary '*' (have 'int')
3550 [[maybe_unused]] const auto &ref_cell = *this;
3551 Assert(false,
3552 (internal::NoPermutation<T>(ref_cell, vertices_0, vertices_1)));
3553 return std::numeric_limits<types::geometric_orientation>::max();
3554 };
3555
3556 switch (this->kind)
3557 {
3559 // TODO: we can get rid of this special-case and use
3560 // vertex_vertex_permutations once we make 0 the default orientation.
3562 if (vertices_0[0] == vertices_1[0])
3564 break;
3571 default:
3573 }
3574
3575 Assert(false, (internal::NoPermutation<T>(*this, vertices_0, vertices_1)));
3576 return std::numeric_limits<types::geometric_orientation>::max();
3577}
3578
3579
3580
3581template <typename T, std::size_t N>
3582inline std::array<T, N>
3584 const std::array<T, N> &vertices,
3585 const unsigned int orientation) const
3586{
3587 Assert(N >= n_vertices(),
3588 ExcMessage("The number of array elements must be equal to or "
3589 "greater than the number of vertices of the cell "
3590 "referenced by this object."));
3591
3592 // Call the non-deprecated function, taking care of calling it only with
3593 // those array elements that we actually care about (see the note
3594 // in the documentation about the arguments potentially being
3595 // larger arrays than necessary).
3596 const auto permutation = permute_by_combined_orientation(
3597 make_array_view(vertices.begin(), vertices.begin() + n_vertices()),
3598 orientation);
3599
3600 std::array<T, N> temp;
3601 std::copy(permutation.begin(), permutation.end(), temp.begin());
3602
3603 return temp;
3604}
3605
3606
3607
3608template <typename T>
3609boost::container::small_vector<T, 8>
3611 const ArrayView<const T> &vertices,
3612 const types::geometric_orientation orientation) const
3613{
3614 AssertDimension(vertices.size(), n_vertices());
3615 boost::container::small_vector<T, 8> result(n_vertices());
3616
3617 auto permute = [&](const auto &table) {
3618 AssertIndexRange(orientation, table.size());
3619 for (unsigned int j = 0; j < table[orientation].size(); ++j)
3620 result[j] = vertices[table[orientation][j]];
3621 };
3622
3623 switch (this->kind)
3624 {
3626 // TODO: we can get rid of this special-case and use
3627 // vertex_vertex_permutations once we make 0 the default orientation.
3628 std::copy(vertices.begin(), vertices.end(), result.begin());
3629 break;
3631 permute(line_vertex_permutations);
3632 break;
3635 break;
3638 break;
3639 default:
3641 }
3642
3643 return result;
3644}
3645
3646
3647
3650 const types::geometric_orientation orientation) const
3651{
3652 switch (this->kind)
3653 {
3655 // Things are always default-oriented in 1D
3656 return orientation;
3657
3659 // the 1d orientations are the identity and a flip: i.e., the identity
3660 // and an involutory mapping
3661 return orientation;
3662
3664 {
3665 AssertIndexRange(orientation, 6);
3666 constexpr std::array<types::geometric_orientation, 6> inverses{
3667 {0, 1, 2, 5, 4, 3}};
3668 return inverses[orientation];
3669 }
3671 {
3672 AssertIndexRange(orientation, 8);
3673 constexpr std::array<types::geometric_orientation, 8> inverses{
3674 {0, 1, 2, 7, 4, 5, 6, 3}};
3675 return inverses[orientation];
3676 }
3677 default:
3679 }
3680
3681 return std::numeric_limits<types::geometric_orientation>::max();
3682}
3683
3684
3685
3686template <>
3687unsigned int
3688ReferenceCell::vtk_lexicographic_to_node_index<0>(
3689 const std::array<unsigned, 0> &node_indices,
3690 const std::array<unsigned, 0> &nodes_per_direction,
3691 const bool legacy_format) const;
3692
3693template <>
3694unsigned int
3695ReferenceCell::vtk_lexicographic_to_node_index<1>(
3696 const std::array<unsigned, 1> &node_indices,
3697 const std::array<unsigned, 1> &nodes_per_direction,
3698 const bool legacy_format) const;
3699
3700template <>
3701unsigned int
3702ReferenceCell::vtk_lexicographic_to_node_index<2>(
3703 const std::array<unsigned, 2> &node_indices,
3704 const std::array<unsigned, 2> &nodes_per_direction,
3705 const bool legacy_format) const;
3706
3707template <>
3708unsigned int
3709ReferenceCell::vtk_lexicographic_to_node_index<3>(
3710 const std::array<unsigned, 3> &node_indices,
3711 const std::array<unsigned, 3> &nodes_per_direction,
3712 const bool legacy_format) const;
3713
3715
3716#endif
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
iterator begin() const
Definition array_view.h:702
iterator end() const
Definition array_view.h:711
std::size_t size() const
Definition array_view.h:684
Abstract base class for mapping classes.
Definition mapping.h:320
Definition point.h:113
static constexpr Point< dim, Number > unit_vector(const unsigned int i)
ArrayView< const unsigned int > faces_for_given_vertex(const unsigned int vertex_index) const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices() const
types::geometric_orientation face_to_cell_line_orientation(const unsigned int face_line_no, const unsigned int face_no, const types::geometric_orientation face_orientation, const types::geometric_orientation line_orientation) const
Convert a line orientation defined relative to a face to the canonical per-cell line orientation.
static constexpr types::geometric_orientation reversed_combined_line_orientation()
unsigned int n_vertices() const
Quadrature< dim > get_gauss_type_quadrature(const unsigned n_points_1d) const
bool is_hyper_cube() const
unsigned int vtk_linear_type() const
void serialize(Archive &archive, const unsigned int)
unsigned int n_isotropic_refinement_choices() const
unsigned int standard_to_real_face_vertex(const unsigned int vertex, const unsigned int face, const types::geometric_orientation face_orientation) const
types::geometric_orientation 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
Point< dim > vertex(const unsigned int v) const
const Quadrature< dim > & get_nodal_type_quadrature() const
Point< dim > face_vertex_location(const unsigned int face, const unsigned int vertex) const
constexpr ::ndarray< unsigned int, 12, 4 > new_isotropic_child_face_lines(const unsigned int refinement_choice) const
unsigned int face_to_cell_lines(const unsigned int face, const unsigned int line, const types::geometric_orientation face_orientation) const
static constexpr ndarray< unsigned int, 1, 1 > vertex_vertex_permutations
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
constexpr ::ndarray< unsigned int, 12, 4, 2 > new_isotropic_child_face_line_vertices(const unsigned int refinement_choice) const
double d_linear_shape_function(const Point< dim > &xi, const unsigned int i) const
types::geometric_orientation standard_vs_true_line_orientation(const unsigned int line, const unsigned int face, const types::geometric_orientation face_orientation, const types::geometric_orientation line_orientation) const
unsigned int child_cell_on_face(const unsigned int face, const unsigned int subface) const
types::geometric_orientation get_inverse_combined_orientation(const types::geometric_orientation orientation) const
Tensor< 1, dim > unit_tangential_vectors(const unsigned int face_no, const unsigned int i) const
static constexpr ndarray< unsigned int, 2, 2 > line_vertex_permutations
unsigned int vtk_lexicographic_to_node_index(const std::array< unsigned, dim > &node_indices, const std::array< unsigned, dim > &nodes_per_direction, const bool legacy_format) const
std::array< T, N > permute_according_orientation(const std::array< T, N > &vertices, const unsigned int orientation) const
constexpr ::ndarray< unsigned int, 8, 4 > new_isotropic_child_cell_vertices(const unsigned int refinement_choice) const
Tensor< 1, dim > face_normal_vector(const unsigned int face_no) const
bool contains_point(const Point< dim > &p, const double tolerance=0) const
constexpr ::ndarray< unsigned int, 8, 6 > new_isotropic_child_cell_faces(const unsigned int refinement_choice) const
unsigned int n_isotropic_children() const
Tensor< 1, dim > face_tangent_vector(const unsigned int face_no, const unsigned int i) const
boost::container::small_vector< T, 8 > permute_by_combined_orientation(const ArrayView< const T > &vertices, const types::geometric_orientation orientation) const
unsigned int gmsh_element_type() const
unsigned int n_face_orientations(const unsigned int face_no) const
static constexpr ndarray< unsigned int, 8, 4 > quadrilateral_vertex_permutations
double volume() const
unsigned int n_faces() const
unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex) const
std::uint8_t kind
constexpr ReferenceCell()
unsigned int exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const
unsigned int unv_vertex_to_deal_vertex(const unsigned int vertex_n) const
static constexpr types::geometric_orientation default_combined_face_orientation()
unsigned int vtk_vertex_to_deal_vertex(const unsigned int vertex_index) const
types::geometric_orientation get_combined_orientation(const ArrayView< const T > &vertices_0, const ArrayView< const T > &vertices_1) const
IsotropicRefinementChoice get_isotropic_refinement_choice(const unsigned int ref_choice) const
unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const types::geometric_orientation face_orientation) 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
unsigned int n_children(const RefinementCase< dim > ref_case=RefinementCase< dim >::isotropic_refinement) 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)
static constexpr std::size_t memory_consumption()
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
constexpr bool operator!=(const ReferenceCell &type) const
constexpr bool operator==(const ReferenceCell &type) const
Point< dim > barycenter() const
unsigned int standard_to_real_face_line(const unsigned int line, const unsigned int face, const types::geometric_orientation face_orientation) 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)
Point< dim > closest_point(const Point< dim > &p) const
static constexpr ndarray< unsigned int, 6, 3 > triangle_vertex_permutations
Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i) const
virtual void print_info(std::ostream &out) const override
const ReferenceCell entity_type
virtual ~NoPermutation() noexcept override=default
const ArrayView< const T > vertices_0
const ArrayView< const T > vertices_1
NoPermutation(const ReferenceCell &entity_type, const ArrayView< const T > &vertices_0, const ArrayView< const T > &vertices_1)
#define DEAL_II_DEPRECATED
Definition config.h:232
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:522
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:523
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcNonMatchingReferenceCellTypes(ReferenceCell arg1, ReferenceCell arg2)
#define Assert(cond, exc)
#define DeclException2(Exception2, type1, type2, outsequence)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
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 unsigned int max_n_lines()
constexpr const ReferenceCell & get_hypercube()
constexpr const ReferenceCell Vertex
constexpr unsigned int max_n_vertices()
constexpr unsigned int max_n_faces()
constexpr const ReferenceCell Line
constexpr const ReferenceCell & get_simplex()
constexpr ReferenceCell make_reference_cell_from_int(const std::uint8_t kind)
std::tuple< bool, bool, bool > split_face_orientation(const types::geometric_orientation combined_face_orientation)
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
constexpr types::geometric_orientation reverse_line_orientation
Definition types.h:359
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:346
boost::integer_range< IncrementableType > iota_view
Definition iota_view.h:45
STL namespace.
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(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)
unsigned char geometric_orientation
Definition types.h:40
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:107
IsotropicRefinementChoice
std::istream & operator>>(std::istream &in, ReferenceCell &reference_cell)
std::ostream & operator<<(std::ostream &out, const ReferenceCell &reference_cell)
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)
static unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static std::array< unsigned int, 2 > standard_hex_line_to_quad_line_index(const unsigned int line)
static unsigned int standard_to_real_face_line(const unsigned int line, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static 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 n_children(const RefinementCase< dim > &refinement_case)
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)