deal.II version GIT relicensing-3128-g792dc1c45b 2025-04-23 21:00: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\}}\)
Loading...
Searching...
No Matches
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
199 template <int dim>
200 double
201 d_linear_shape_function(const Point<dim> &xi, const unsigned int i) const;
202
211 template <int dim>
214 const unsigned int i) const;
215
225 template <int dim, int spacedim = dim>
226 std::unique_ptr<Mapping<dim, spacedim>>
227 get_default_mapping(const unsigned int degree) const;
228
240 template <int dim, int spacedim = dim>
243
256 template <int dim>
258 get_gauss_type_quadrature(const unsigned n_points_1d) const;
259
273 template <int dim>
276
294 template <int dim>
295 const Quadrature<dim> &
297
312 unsigned int
313 n_vertices() const;
314
320 vertex_indices() const;
321
333 template <int dim>
335 vertex(const unsigned int v) const;
336
342 unsigned int
343 n_lines() const;
344
350 line_indices() const;
351
357 unsigned int
358 n_faces() const;
359
365 face_indices() const;
366
370 unsigned int
372
377 get_isotropic_refinement_choice(const unsigned int ref_choice) const;
378
390 unsigned int
391 n_isotropic_children() const;
392
401 template <int dim>
402 unsigned int
403 n_children(const RefinementCase<dim> ref_case =
405
412
426 face_reference_cell(const unsigned int face_no) const;
427
443 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
444 "Use numbers::default_geometric_orientation instead.")
445 static constexpr types::geometric_orientation
447
456 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
457 "Use numbers::reverse_line_orientation instead.")
458 static constexpr types::geometric_orientation
460
468 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
469 "Use the version of this function which takes the orientation as the third "
470 "parameter instead.")
471 unsigned int
472 child_cell_on_face(const unsigned int face, const unsigned int subface) const;
473
515 unsigned int
516 child_cell_on_face(const unsigned int face,
517 const unsigned int subface,
518 const types::geometric_orientation face_orientation) const;
519
529 std::array<unsigned int, 2>
530 standard_vertex_to_face_and_vertex_index(const unsigned int vertex) const;
531
541 std::array<unsigned int, 2>
542 standard_line_to_face_and_line_index(const unsigned int line) const;
543
552 unsigned int
553 line_to_cell_vertices(const unsigned int line,
554 const unsigned int vertex) const;
555
559 unsigned int
560 face_to_cell_lines(const unsigned int face,
561 const unsigned int line,
562 const types::geometric_orientation face_orientation) const;
563
567 unsigned int
569 const unsigned int face,
570 const unsigned int vertex,
571 const types::geometric_orientation face_orientation) const;
572
596 template <int dim>
598 face_vertex_location(const unsigned int face,
599 const unsigned int vertex) const;
600
604 unsigned int
606 const unsigned int vertex,
607 const unsigned int face,
608 const types::geometric_orientation face_orientation) const;
609
613 unsigned int
615 const unsigned int line,
616 const unsigned int face,
617 const types::geometric_orientation face_orientation) const;
618
626 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
627 "Use face_to_cell_line_orientation() instead.")
630 const unsigned int line,
631 const unsigned int face,
632 const types::geometric_orientation face_orientation,
633 const types::geometric_orientation line_orientation) const;
634
659 const unsigned int face_line_no,
660 const unsigned int face_no,
661 const types::geometric_orientation face_orientation,
662 const types::geometric_orientation line_orientation) const;
663
687 double
688 volume() const;
689
705 double
706 face_measure(const unsigned int face_no) const;
707
724 template <int dim>
726 barycenter() const;
727
751 template <int dim>
752 bool
753 contains_point(const Point<dim> &p, const double tolerance = 0) const;
754
763 template <int dim>
765 closest_point(const Point<dim> &p) const;
766
778 template <int dim>
780 face_tangent_vector(const unsigned int face_no, const unsigned int i) const;
781
795 template <int dim>
796 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT("Use face_tangent_vector() instead.")
797 Tensor<1, dim> unit_tangential_vectors(const unsigned int face_no,
798 const unsigned int i) const;
799
807 template <int dim>
809 face_normal_vector(const unsigned int face_no) const;
810
820 template <int dim>
821 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT("Use face_normal_vector() instead.")
822 Tensor<1, dim> unit_normal_vectors(const unsigned int face_no) const;
823
829 unsigned int
830 n_face_orientations(const unsigned int face_no) const;
831
848 template <typename T, std::size_t N>
850 compute_orientation(const std::array<T, N> &vertices_0,
851 const std::array<T, N> &vertices_1) const;
852
878 template <typename T>
881 const ArrayView<const T> &vertices_1) const;
882
883
897 template <typename T, std::size_t N>
898 DEAL_II_DEPRECATED std::array<T, N>
899 permute_according_orientation(const std::array<T, N> &vertices,
900 const unsigned int orientation) const;
901
913 template <typename T>
914 boost::container::small_vector<T, 8>
916 const ArrayView<const T> &vertices,
917 const types::geometric_orientation orientation) const;
918
926 const types::geometric_orientation orientation) const;
927
932 faces_for_given_vertex(const unsigned int vertex_index) const;
933
938 constexpr ::ndarray<unsigned int, 12, 4>
939 new_isotropic_child_face_lines(const unsigned int refinement_choice) const;
940
945 constexpr ::ndarray<unsigned int, 12, 4, 2>
947 const unsigned int refinement_choice) const;
948
953 constexpr ::ndarray<unsigned int, 8, 6>
954 new_isotropic_child_cell_faces(const unsigned int refinement_choice) const;
955
960 constexpr ::ndarray<unsigned int, 8, 4>
961 new_isotropic_child_cell_vertices(const unsigned int refinement_choice) const;
962
975 unsigned int
976 exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const;
977
981 unsigned int
982 exodusii_face_to_deal_face(const unsigned int face_n) const;
983
987 unsigned int
988 unv_vertex_to_deal_vertex(const unsigned int vertex_n) const;
989
993 unsigned int
994 vtk_linear_type() const;
995
1000 unsigned int
1001 vtk_quadratic_type() const;
1002
1007 unsigned int
1008 vtk_lagrange_type() const;
1009
1026 template <int dim>
1027 unsigned int
1029 const std::array<unsigned, dim> &node_indices,
1030 const std::array<unsigned, dim> &nodes_per_direction,
1031 const bool legacy_format) const;
1032
1036 unsigned int
1037 vtk_vertex_to_deal_vertex(const unsigned int vertex_index) const;
1038
1042 unsigned int
1043 gmsh_element_type() const;
1044
1058 std::string
1059 to_string() const;
1060
1064 constexpr operator std::uint8_t() const;
1065
1069 constexpr bool
1070 operator==(const ReferenceCell &type) const;
1071
1075 constexpr bool
1076 operator!=(const ReferenceCell &type) const;
1077
1083 template <class Archive>
1084 void
1085 serialize(Archive &archive, const unsigned int /*version*/);
1086
1090 static constexpr std::size_t
1092
1097private:
1101 std::uint8_t kind;
1102
1109 constexpr ReferenceCell(const std::uint8_t kind);
1110
1117 {{{0}}}};
1118
1123 {{{0, 1}}, {{1, 0}}}};
1124
1129 {{{0, 1, 2}},
1130 {{0, 2, 1}},
1131 {{2, 0, 1}},
1132 {{2, 1, 0}},
1133 {{1, 2, 0}},
1134 {{1, 0, 2}}}};
1135
1139 static constexpr ndarray<unsigned int, 8, 4>
1141 {{0, 1, 2, 3}},
1142 {{0, 2, 1, 3}},
1143 {{2, 0, 3, 1}},
1144 {{2, 3, 0, 1}},
1145 {{3, 2, 1, 0}},
1146 {{3, 1, 2, 0}},
1147 {{1, 3, 0, 2}},
1148 {{1, 0, 3, 2}},
1149 }};
1150
1155 friend constexpr ReferenceCell
1157
1158 friend std::ostream &
1159 operator<<(std::ostream &out, const ReferenceCell &reference_cell);
1160
1161 friend std::istream &
1162 operator>>(std::istream &in, ReferenceCell &reference_cell);
1163};
1164
1165
1173std::ostream &
1174operator<<(std::ostream &out, const ReferenceCell &reference_cell);
1175
1183std::istream &
1184operator>>(std::istream &in, ReferenceCell &reference_cell);
1185
1186
1187inline constexpr ReferenceCell::ReferenceCell(const std::uint8_t kind)
1188 : kind(kind)
1189{}
1190
1191
1192
1193inline constexpr ReferenceCell::operator std::uint8_t() const
1194{
1195 return kind;
1196}
1197
1198
1199
1200inline constexpr bool
1202{
1203 return kind == type.kind;
1204}
1205
1206
1207
1208inline constexpr bool
1210{
1211 return kind != type.kind;
1212}
1213
1214
1215
1216namespace internal
1217{
1218 inline constexpr ReferenceCell
1219 make_reference_cell_from_int(const std::uint8_t kind)
1220 {
1221#ifndef DEAL_II_CXX14_CONSTEXPR_BUG
1222 // Make sure these are the only indices from which objects can be
1223 // created.
1224 Assert((kind == std::numeric_limits<std::uint8_t>::max()) || (kind < 8),
1226#endif
1227
1228 // Call the private constructor, which we can from here because this
1229 // function is a 'friend'.
1230 return {kind};
1231 }
1232} // namespace internal
1233
1234
1235
1244{
1257 std::numeric_limits<std::uint8_t>::max());
1258
1264 template <int dim>
1265 constexpr const ReferenceCell &
1266 get_simplex();
1267
1273 template <int dim>
1274 constexpr const ReferenceCell &
1275 get_hypercube();
1276
1284 template <int structdim>
1285 constexpr unsigned int
1287
1295 template <int structdim>
1296 constexpr unsigned int
1297 max_n_lines();
1298
1314 template <int structdim>
1315 constexpr unsigned int
1316 max_n_faces();
1317} // namespace ReferenceCells
1318
1319
1320
1322 : ReferenceCell(ReferenceCells::Invalid)
1323{}
1324
1325
1326
1327template <class Archive>
1328inline void
1329ReferenceCell::serialize(Archive &archive, const unsigned int /*version*/)
1330{
1331 archive &kind;
1332}
1333
1334
1335
1336inline constexpr std::size_t
1341
1342
1343
1345ReferenceCell::faces_for_given_vertex(const unsigned int vertex) const
1346{
1348 switch (this->kind)
1349 {
1351 return {&GeometryInfo<2>::vertex_to_face[vertex][0], 1};
1353 return {&GeometryInfo<2>::vertex_to_face[vertex][0], 2};
1355 {
1356 static constexpr ndarray<unsigned int, 3, 2> table = {
1357 {{{0, 2}}, {{0, 1}}, {{1, 2}}}};
1358 return table[vertex];
1359 }
1361 {
1362 static constexpr ndarray<unsigned int, 4, 3> table = {
1363 {{{0, 1, 2}}, {{0, 1, 3}}, {{0, 2, 3}}, {{1, 2, 3}}}};
1364
1365 return table[vertex];
1366 }
1368 {
1369 static constexpr unsigned int X = numbers::invalid_unsigned_int;
1370 static constexpr ndarray<unsigned int, 5, 4> table = {
1371 {{{0, 1, 3, X}},
1372 {{0, 2, 3, X}},
1373 {{0, 1, 4, X}},
1374 {{0, 2, 4, X}},
1375 {{1, 2, 3, 4}}}};
1376
1377 return {&table[vertex][0], vertex == 4 ? 4u : 3u};
1378 }
1380 {
1382 static constexpr ndarray<unsigned int, 6, 3> table = {{{{0, 2, 4}},
1383 {{0, 2, 3}},
1384 {{0, 3, 4}},
1385 {{1, 2, 4}},
1386 {{1, 2, 3}},
1387 {{1, 3, 4}}}};
1388
1389 return table[vertex];
1390 }
1392 return {&GeometryInfo<3>::vertex_to_face[vertex][0], 3};
1393 default:
1395 }
1396
1397 return {};
1398}
1399
1400
1401constexpr ::ndarray<unsigned int, 12, 4>
1403 const unsigned int refienement_choice) const
1404{
1405 // AssertIndexRange(refienement_choice, n_isotropic_refinement_choices());
1406 const unsigned int X = numbers::invalid_unsigned_int;
1407 switch (this->kind)
1408 {
1410 {
1411 constexpr ::ndarray<unsigned int, 3, 12, 4> new_quad_lines_tet =
1412 {{// new line is (6,8)
1413 {{{{2, 3, 8, X}},
1414 {{0, 9, 5, X}},
1415 {{1, 6, 11, X}},
1416 {{4, 10, 7, X}},
1417 {{2, 12, 5, X}},
1418 {{1, 9, 12, X}},
1419 {{4, 8, 12, X}},
1420 {{6, 12, 10, X}},
1421 {{X, X, X, X}},
1422 {{X, X, X, X}},
1423 {{X, X, X, X}},
1424 {{X, X, X, X}}}},
1425 // new line is (5,7)
1426 {{{{2, 3, 8, X}},
1427 {{0, 9, 5, X}},
1428 {{1, 6, 11, X}},
1429 {{4, 10, 7, X}},
1430 {{0, 3, 12, X}},
1431 {{1, 12, 8, X}},
1432 {{4, 12, 9, X}},
1433 {{7, 11, 12, X}},
1434 {{X, X, X, X}},
1435 {{X, X, X, X}},
1436 {{X, X, X, X}},
1437 {{X, X, X, X}}}},
1438 // new line is (4,9)
1439 {{{{2, 3, 8, X}},
1440 {{0, 9, 5, X}},
1441 {{1, 6, 11, X}},
1442 {{4, 10, 7, X}},
1443 {{0, 12, 11, X}},
1444 {{2, 6, 12, X}},
1445 {{3, 12, 7, X}},
1446 {{5, 10, 12, X}},
1447 {{X, X, X, X}},
1448 {{X, X, X, X}},
1449 {{X, X, X, X}},
1450 {{X, X, X, X}}}}}};
1451 return new_quad_lines_tet[refienement_choice];
1452 }
1453
1455 {
1456 constexpr ::ndarray<unsigned int, 12, 4> new_quad_lines_hex = {
1457 {{{10, 28, 16, 24}},
1458 {{28, 14, 17, 25}},
1459 {{11, 29, 24, 20}},
1460 {{29, 15, 25, 21}},
1461 {{18, 26, 0, 28}},
1462 {{26, 22, 1, 29}},
1463 {{19, 27, 28, 4}},
1464 {{27, 23, 29, 5}},
1465 {{2, 24, 8, 26}},
1466 {{24, 6, 9, 27}},
1467 {{3, 25, 26, 12}},
1468 {{25, 7, 27, 13}}}};
1469 return new_quad_lines_hex;
1470 }
1471 default:
1473 }
1474
1475 return {};
1476}
1477
1478
1479
1480constexpr ::ndarray<unsigned int, 12, 4, 2>
1482 const unsigned int refienement_choice) const
1483{
1484 // AssertIndexRange(refienement_choice, n_isotropic_refinement_choices());
1485 const unsigned int X = numbers::invalid_unsigned_int;
1486 switch (this->kind)
1487 {
1489 {
1490 constexpr ::ndarray<unsigned int, 3, 12, 4, 2> table_tet = {
1491 {// new line is (6, 8)
1492 {{{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
1493 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
1494 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
1495 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
1496 {{{{4, 6}}, {{6, 8}}, {{8, 4}}, {{X, X}}}},
1497 {{{{6, 5}}, {{5, 8}}, {{8, 6}}, {{X, X}}}},
1498 {{{{8, 7}}, {{7, 6}}, {{6, 8}}, {{X, X}}}},
1499 {{{{9, 6}}, {{6, 8}}, {{8, 9}}, {{X, X}}}},
1500 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1501 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1502 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1503 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}},
1504 // new line is (5, 7)
1505 {{{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
1506 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
1507 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
1508 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
1509 {{{{5, 4}}, {{4, 7}}, {{7, 5}}, {{X, X}}}},
1510 {{{{6, 5}}, {{5, 7}}, {{7, 6}}, {{X, X}}}},
1511 {{{{8, 7}}, {{7, 5}}, {{5, 8}}, {{X, X}}}},
1512 {{{{7, 9}}, {{9, 5}}, {{5, 7}}, {{X, X}}}},
1513 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1514 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1515 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1516 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}},
1517 // new line is (4, 9)
1518 {{{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
1519 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
1520 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
1521 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
1522 {{{{5, 4}}, {{4, 9}}, {{9, 5}}, {{X, X}}}},
1523 {{{{4, 6}}, {{6, 9}}, {{9, 4}}, {{X, X}}}},
1524 {{{{7, 4}}, {{4, 9}}, {{9, 7}}, {{X, X}}}},
1525 {{{{4, 8}}, {{8, 9}}, {{9, 4}}, {{X, X}}}},
1526 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1527 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1528 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1529 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}}}};
1530 return table_tet[refienement_choice];
1531 }
1532
1534 {
1535 constexpr ::ndarray<unsigned int, 12, 4, 2> table_hex = {
1536 {{{{{10, 22}}, {{24, 26}}, {{10, 24}}, {{22, 26}}}},
1537 {{{{24, 26}}, {{11, 23}}, {{24, 11}}, {{26, 23}}}},
1538 {{{{22, 14}}, {{26, 25}}, {{22, 26}}, {{14, 25}}}},
1539 {{{{26, 25}}, {{23, 15}}, {{26, 23}}, {{25, 15}}}},
1540 {{{{8, 24}}, {{20, 26}}, {{8, 20}}, {{24, 26}}}},
1541 {{{{20, 26}}, {{12, 25}}, {{20, 12}}, {{26, 25}}}},
1542 {{{{24, 9}}, {{26, 21}}, {{24, 26}}, {{9, 21}}}},
1543 {{{{26, 21}}, {{25, 13}}, {{26, 25}}, {{21, 13}}}},
1544 {{{{16, 20}}, {{22, 26}}, {{16, 22}}, {{20, 26}}}},
1545 {{{{22, 26}}, {{17, 21}}, {{22, 17}}, {{26, 21}}}},
1546 {{{{20, 18}}, {{26, 23}}, {{20, 26}}, {{18, 23}}}},
1547 {{{{26, 23}}, {{21, 19}}, {{26, 21}}, {{23, 19}}}}}};
1548 return table_hex;
1549 }
1550 default:
1552 }
1553
1554 return {};
1555}
1556
1557
1558
1559constexpr ::ndarray<unsigned int, 8, 6>
1561 const unsigned int refienement_choice) const
1562{
1563 // AssertIndexRange(refienement_choice, n_isotropic_refinement_choices());
1564 const unsigned int X = numbers::invalid_unsigned_int;
1565 switch (this->kind)
1566 {
1568 {
1569 constexpr ::ndarray<unsigned int, 3, 8, 6> cell_quads_tet = {
1570 {// new line is (6, 8)
1571 {{{{8, 13, 16, 0, X, X}},
1572 {{9, 12, 1, 21, X, X}},
1573 {{10, 2, 17, 20, X, X}},
1574 {{3, 14, 18, 22, X, X}},
1575 {{11, 1, 4, 5, X, X}},
1576 {{15, 0, 4, 6, X, X}},
1577 {{19, 7, 6, 3, X, X}},
1578 {{23, 5, 2, 7, X, X}}}},
1579 // new line is (5, 7)
1580 {{
1581 {{8, 13, 16, 0, X, X}},
1582 {{9, 12, 1, 21, X, X}},
1583 {{10, 2, 17, 20, X, X}},
1584 {{3, 14, 18, 22, X, X}},
1585 {{11, 4, 0, 5, X, X}},
1586 {{15, 4, 1, 6, X, X}},
1587 {{19, 2, 5, 7, X, X}},
1588 {{23, 6, 7, 3, X, X}},
1589 }},
1590 // new line is (4, 9)
1591 {{
1592 {{8, 13, 16, 0, X, X}},
1593 {{9, 12, 1, 21, X, X}},
1594 {{10, 2, 17, 20, X, X}},
1595 {{3, 14, 18, 22, X, X}},
1596 {{11, 4, 5, 2, X, X}},
1597 {{15, 6, 7, 3, X, X}},
1598 {{19, 5, 0, 6, X, X}},
1599 {{23, 1, 4, 7, X, X}},
1600 }}}};
1601 return cell_quads_tet[refienement_choice];
1602 }
1603
1605 {
1606 constexpr ::ndarray<unsigned int, 8, 6> cell_quads_hex = {{
1607 {{12, 0, 20, 4, 28, 8}}, // bottom children
1608 {{0, 16, 22, 6, 29, 9}}, //
1609 {{13, 1, 4, 24, 30, 10}}, //
1610 {{1, 17, 6, 26, 31, 11}}, //
1611 {{14, 2, 21, 5, 8, 32}}, // top children
1612 {{2, 18, 23, 7, 9, 33}}, //
1613 {{15, 3, 5, 25, 10, 34}}, //
1614 {{3, 19, 7, 27, 11, 35}} //
1615 }};
1616 return cell_quads_hex;
1617 }
1618 default:
1620 }
1621
1622 return {};
1623}
1624
1625
1626
1627constexpr ::ndarray<unsigned int, 8, 4>
1629 const unsigned int refienement_choice) const
1630{
1631 // AssertIndexRange(refienement_choice, n_isotropic_refinement_choices());
1632
1633 switch (this->kind)
1634 {
1636 {
1637 constexpr ::ndarray<unsigned int, 3, 8, 4> cell_vertices_tet = {
1638 {// new line is (6,8)
1639 {{{{0, 4, 6, 7}},
1640 {{4, 1, 5, 8}},
1641 {{6, 5, 2, 9}},
1642 {{7, 8, 9, 3}},
1643 {{4, 5, 6, 8}},
1644 {{4, 7, 8, 6}},
1645 {{6, 9, 7, 8}},
1646 {{5, 8, 9, 6}}}},
1647 // new line is (5,7)
1648 {{{{0, 4, 6, 7}},
1649 {{4, 1, 5, 8}},
1650 {{6, 5, 2, 9}},
1651 {{7, 8, 9, 3}},
1652 {{4, 5, 6, 7}},
1653 {{4, 7, 8, 5}},
1654 {{6, 9, 7, 5}},
1655 {{5, 8, 9, 7}}}},
1656 // new line is (4,9)
1657 {{{{0, 4, 6, 7}},
1658 {{4, 1, 5, 8}},
1659 {{6, 5, 2, 9}},
1660 {{7, 8, 9, 3}},
1661 {{4, 5, 6, 9}},
1662 {{4, 7, 8, 9}},
1663 {{6, 9, 7, 4}},
1664 {{5, 8, 9, 4}}}}}};
1665 return cell_vertices_tet[refienement_choice];
1666 }
1667 default:
1669 }
1670
1671 return {};
1672}
1673
1674
1675
1676inline bool
1678{
1679 return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
1682}
1683
1684
1685
1686inline bool
1688{
1689 return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
1690 *this == ReferenceCells::Triangle ||
1692}
1693
1694
1695
1696inline unsigned int
1698{
1699 switch (this->kind)
1700 {
1702 return 0;
1704 return 1;
1707 return 2;
1712 return 3;
1713 default:
1715 }
1716
1718}
1719
1720
1721
1722template <int dim>
1725{
1726 Assert(dim == get_dimension(),
1727 ExcMessage("You can only call this function with arguments for "
1728 "which 'dim' equals the dimension of the space in which "
1729 "the current reference cell lives. You have dim=" +
1730 std::to_string(dim) + " but the reference cell is " +
1731 std::to_string(get_dimension()) + " dimensional."));
1732
1733 return Quadrature<dim>(std::vector<Point<dim>>({barycenter<dim>()}),
1734 std::vector<double>({volume()}));
1735}
1736
1737
1738
1739inline unsigned int
1741{
1742 switch (this->kind)
1743 {
1745 return 1;
1747 return 2;
1749 return 3;
1751 return 4;
1753 return 4;
1755 return 5;
1757 return 6;
1759 return 8;
1760 default:
1762 }
1763
1765}
1766
1767
1768
1769inline unsigned int
1771{
1772 switch (this->kind)
1773 {
1775 return 0;
1777 return 1;
1779 return 3;
1781 return 4;
1783 return 6;
1785 return 8;
1787 return 9;
1789 return 12;
1790 default:
1792 }
1793
1795}
1796
1797
1798
1799template <int dim>
1801ReferenceCell::vertex(const unsigned int v) const
1802{
1803 Assert(dim == get_dimension(),
1804 ExcMessage("You can only call this function with arguments for "
1805 "which 'dim' equals the dimension of the space in which "
1806 "the current reference cell lives. You have dim=" +
1807 std::to_string(dim) + " but the reference cell is " +
1808 std::to_string(get_dimension()) + " dimensional."));
1810
1811 switch (dim)
1812 {
1813 case 0:
1814 {
1815 if (*this == ReferenceCells::Vertex)
1816 return Point<dim>();
1817 break;
1818 }
1819 case 1:
1820 {
1821 static const Point<dim> vertices[2] = {
1822 Point<dim>(), // the origin
1823 Point<dim>::unit_vector(0) // unit point along x-axis
1824 };
1825 if (*this == ReferenceCells::Line)
1826 return vertices[v];
1827 break;
1828 }
1829 case 2:
1830 {
1831 switch (this->kind)
1832 {
1834 {
1835 static const Point<dim> vertices[3] = {
1836 Point<dim>(), // the origin
1837 Point<dim>::unit_vector(0), // unit point along x-axis
1838 Point<dim>::unit_vector(1) // unit point along y-axis
1839 };
1840 return vertices[v];
1841 }
1843 {
1844 static const Point<dim> vertices[4] = {
1845 // First the two points on the x-axis
1846 Point<dim>(),
1848 // Then these two points shifted in the y-direction
1851 return vertices[v];
1852 }
1853 }
1854 break;
1855 }
1856 case 3:
1857 {
1858 switch (this->kind)
1859 {
1861 {
1862 static const Point<dim> vertices[4] = {
1863 Point<dim>(), // the origin
1864 Point<dim>::unit_vector(0), // unit point along x-axis
1865 Point<dim>::unit_vector(1), // unit point along y-axis
1866 Point<dim>::unit_vector(2) // unit point along z-axis
1867 };
1868 return vertices[v];
1869 }
1871 {
1872 static const Point<dim> vertices[5] = {
1873 Point<dim>{-1.0, -1.0, 0.0},
1874 Point<dim>{+1.0, -1.0, 0.0},
1875 Point<dim>{-1.0, +1.0, 0.0},
1876 Point<dim>{+1.0, +1.0, 0.0},
1877 Point<dim>{+0.0, +0.0, 1.0}};
1878 return vertices[v];
1879 }
1881 {
1882 static const Point<dim> vertices[6] = {
1883 // First the three points on the triangular base of the
1884 // wedge:
1885 Point<dim>(),
1888 // And now everything shifted in the z-direction again
1892 return vertices[v];
1893 }
1895 {
1896 static const Point<dim> vertices[8] = {
1897 // First the two points on the x-axis
1898 Point<dim>(),
1900 // Then these two points shifted in the y-direction
1903 // And now all four points shifted in the z-direction
1910 return vertices[v];
1911 }
1912 }
1913 break;
1914 }
1915 default:
1917 }
1918
1920 return Point<dim>();
1921}
1922
1923
1924inline unsigned int
1926{
1927 switch (this->kind)
1928 {
1930 return 0;
1932 return 2;
1934 return 3;
1936 return 4;
1938 return 4;
1940 return 5;
1942 return 5;
1944 return 6;
1945 default:
1947 }
1948
1950}
1951
1952
1953
1960
1961
1962
1963inline unsigned int
1965{
1966 switch (this->kind)
1967 {
1969 return 1;
1971 return 1;
1973 return 1;
1975 return 1;
1977 return 3;
1979 return 1;
1981 return 1;
1983 return 1;
1984 default:
1986 }
1987
1989}
1990
1991
1992
1995 const unsigned int ref_choice) const
1996{
1998 switch (this->kind)
1999 {
2001 {
2003 isotropic_ref_choices = {{IsotropicRefinementChoice::cut_tet_68,
2006 return isotropic_ref_choices[ref_choice];
2007 }
2008 default:
2009 {
2011 }
2012 }
2013
2015}
2016
2017
2018
2019inline unsigned int
2021{
2022 switch (this->kind)
2023 {
2025 return 0;
2027 return 2;
2029 return 4;
2031 return 4;
2033 return 8;
2035 // We haven't yet decided how to refine pyramids. Update this when we
2036 // have
2037 return 0;
2039 return 8;
2041 return 8;
2042 default:
2044 }
2045
2047}
2048
2049
2050
2051template <int dim>
2052unsigned int
2054{
2055 Assert(dim == get_dimension(),
2056 ExcMessage("You can only call this function with arguments for "
2057 "which 'dim' equals the dimension of the space in which "
2058 "the current reference cell lives. You have dim=" +
2059 std::to_string(dim) + " but the reference cell is " +
2060 std::to_string(get_dimension()) + " dimensional."));
2061
2062 // Use GeometryInfo here to keep it the single source of truth
2063 if (this->is_hyper_cube())
2064 return GeometryInfo<dim>::n_children(ref_case);
2065 else
2066 return this->n_isotropic_children();
2067}
2068
2069
2070
2077
2078
2079
2086
2087
2088
2095
2096
2097
2098inline ReferenceCell
2099ReferenceCell::face_reference_cell(const unsigned int face_no) const
2100{
2101 AssertIndexRange(face_no, n_faces());
2102
2103 switch (this->kind)
2104 {
2111 return ReferenceCells::Line;
2115 if (face_no == 0)
2117 else
2120 if (face_no > 1)
2122 else
2126 default:
2128 }
2129
2131}
2132
2133
2134
2135inline constexpr types::geometric_orientation
2140
2141
2142
2143inline constexpr types::geometric_orientation
2145{
2146 // For a reversed line 'orientation' is false and neither flip nor rotate are
2147 // defined.
2149}
2150
2151
2152
2153inline unsigned int
2155 const unsigned int subface) const
2156{
2157 return child_cell_on_face(face,
2158 subface,
2160}
2161
2162
2163
2164inline unsigned int
2166 const unsigned int face,
2167 const unsigned int subface,
2168 const types::geometric_orientation combined_face_orientation) const
2169{
2170 AssertIndexRange(face, n_faces());
2172
2173 switch (this->kind)
2174 {
2177 {
2179 break;
2180 }
2182 {
2183 static constexpr ndarray<unsigned int, 3, 2> subcells = {
2184 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
2185
2186 Assert(combined_face_orientation ==
2188 combined_face_orientation ==
2191 return subcells[face][combined_face_orientation ==
2193 subface :
2194 1 - subface];
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 }
2212 {
2214 break;
2215 }
2217 {
2218 const auto [face_orientation, face_rotation, face_flip] =
2219 internal::split_face_orientation(combined_face_orientation);
2220
2223 face,
2224 subface,
2225 face_orientation,
2226 face_flip,
2227 face_rotation);
2228 }
2229 default:
2231 }
2232
2234}
2235
2236
2237
2238inline std::array<unsigned int, 2>
2240 const unsigned int vertex) const
2241{
2243 // Work around a GCC warning at higher optimization levels by making all of
2244 // these tables the same size
2245 constexpr unsigned int X = numbers::invalid_unsigned_int;
2246
2247 switch (this->kind)
2248 {
2252 break;
2254 {
2255 static constexpr ndarray<unsigned int, 6, 2> table = {
2256 {{{0, 0}}, {{0, 1}}, {{1, 1}}, {{X, X}}, {{X, X}}, {{X, X}}}};
2257
2258 return table[vertex];
2259 }
2261 {
2263 vertex);
2264 }
2266 {
2267 static constexpr ndarray<unsigned int, 6, 2> table = {
2268 {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 2}}, {{X, X}}, {{X, X}}}};
2269
2270 return table[vertex];
2271 }
2273 {
2274 static constexpr ndarray<unsigned int, 6, 2> table = {
2275 {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{0, 3}}, {{1, 2}}, {{X, X}}}};
2276
2277 return table[vertex];
2278 }
2280 {
2281 static constexpr ndarray<unsigned int, 6, 2> table = {
2282 {{{0, 1}}, {{0, 0}}, {{0, 2}}, {{1, 0}}, {{1, 1}}, {{1, 2}}}};
2283
2284 return table[vertex];
2285 }
2287 {
2289 vertex);
2290 }
2291 default:
2293 }
2294
2295 return {};
2296}
2297
2298
2299
2300inline std::array<unsigned int, 2>
2302 const unsigned int line) const
2303{
2304 AssertIndexRange(line, n_lines());
2305
2306 switch (this->kind)
2307 {
2312 {
2314 break;
2315 }
2317 {
2318 static const std::array<unsigned int, 2> table[6] = {
2319 {{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}};
2320
2321 return table[line];
2322 }
2324 {
2325 static const std::array<unsigned int, 2> table[8] = {{{0, 0}},
2326 {{0, 1}},
2327 {{0, 2}},
2328 {{0, 3}},
2329 {{1, 2}},
2330 {{2, 1}},
2331 {{1, 1}},
2332 {{2, 2}}};
2333
2334 return table[line];
2335 }
2337 {
2338 static const std::array<unsigned int, 2> table[9] = {{{0, 0}},
2339 {{0, 2}},
2340 {{0, 1}},
2341 {{1, 0}},
2342 {{1, 1}},
2343 {{1, 2}},
2344 {{2, 0}},
2345 {{2, 1}},
2346 {{3, 1}}};
2347
2348 return table[line];
2349 }
2351 {
2353 }
2354 default:
2356 }
2357
2358 return {};
2359}
2360
2361
2362inline unsigned int
2364 const unsigned int vertex) const
2365{
2367 AssertIndexRange(line, n_lines());
2368
2369 switch (this->kind)
2370 {
2373 return vertex;
2375 {
2376 static constexpr ndarray<unsigned int, 3, 2> table = {
2377 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
2378 return table[line][vertex];
2379 }
2381 {
2382 static constexpr ndarray<unsigned int, 4, 2> table = {
2383 {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
2384 return table[line][vertex];
2385 }
2387 {
2388 static constexpr ndarray<unsigned int, 6, 2> table = {
2389 {{{0, 1}}, {{1, 2}}, {{2, 0}}, {{0, 3}}, {{1, 3}}, {{2, 3}}}};
2390 return table[line][vertex];
2391 }
2393 {
2394 static constexpr ndarray<unsigned int, 8, 2> table = {{{{0, 2}},
2395 {{1, 3}},
2396 {{0, 1}},
2397 {{2, 3}},
2398 {{4, 0}},
2399 {{1, 4}},
2400 {{2, 4}},
2401 {{4, 3}}}};
2402 return table[line][vertex];
2403 }
2405 {
2406 static constexpr ndarray<unsigned int, 9, 2> table = {{{{1, 0}},
2407 {{2, 1}},
2408 {{0, 2}},
2409 {{3, 4}},
2410 {{4, 5}},
2411 {{5, 3}},
2412 {{0, 3}},
2413 {{1, 4}},
2414 {{2, 5}}}};
2415 return table[line][vertex];
2416 }
2418 {
2419 // first four lines comprise the bottom face, next four are the top,
2420 // and the last four are 'bottom to top'
2421 static constexpr ndarray<unsigned int, 12, 2> table = {{{{0, 2}},
2422 {{1, 3}},
2423 {{0, 1}},
2424 {{2, 3}},
2425 {{4, 6}},
2426 {{5, 7}},
2427 {{4, 5}},
2428 {{6, 7}},
2429 {{0, 4}},
2430 {{1, 5}},
2431 {{2, 6}},
2432 {{3, 7}}}};
2433 return table[line][vertex];
2434 }
2435
2436 default:
2438 }
2439
2441}
2442
2443
2444inline unsigned int
2446 const unsigned int face,
2447 const unsigned int line,
2448 const types::geometric_orientation combined_face_orientation) const
2449{
2450 AssertIndexRange(face, n_faces());
2452
2453 static constexpr unsigned int X = numbers::invalid_unsigned_int;
2454
2455 switch (this->kind)
2456 {
2458 {
2460 break;
2461 }
2463 {
2464 const auto [face_orientation, face_rotation, face_flip] =
2465 internal::split_face_orientation(combined_face_orientation);
2466
2468 face, line, face_orientation, face_flip, face_rotation);
2469 }
2471 {
2472 return face;
2473 }
2475 {
2476 const auto [face_orientation, face_rotation, face_flip] =
2477 internal::split_face_orientation(combined_face_orientation);
2478
2480 face, line, face_orientation, face_flip, face_rotation);
2481 }
2483 {
2484 static constexpr ndarray<unsigned int, 4, 3> table = {
2485 {{{0, 1, 2}}, {{0, 3, 4}}, {{2, 5, 3}}, {{1, 4, 5}}}};
2486
2487 return table[face][standard_to_real_face_line(
2488 line, face, combined_face_orientation)];
2489 }
2491 {
2492 static constexpr ndarray<unsigned int, 5, 4> table = {
2493 {{{0, 1, 2, 3}},
2494 {{0, 6, 4, X}},
2495 {{1, 5, 7, X}},
2496 {{2, 4, 5, X}},
2497 {{3, 7, 6, X}}}};
2498
2499 return table[face][standard_to_real_face_line(
2500 line, face, combined_face_orientation)];
2501 }
2503 {
2504 static constexpr ndarray<unsigned int, 5, 4> table = {
2505 {{{0, 2, 1, X}},
2506 {{3, 4, 5, X}},
2507 {{6, 7, 0, 3}},
2508 {{7, 8, 1, 4}},
2509 {{8, 6, 5, 2}}}};
2510
2511 return table[face][standard_to_real_face_line(
2512 line, face, combined_face_orientation)];
2513 }
2515 {
2516 const auto [face_orientation, face_rotation, face_flip] =
2517 internal::split_face_orientation(combined_face_orientation);
2518
2520 face, line, face_orientation, face_flip, face_rotation);
2521 }
2522 default:
2524 }
2525
2527}
2528
2529
2530
2531inline unsigned int
2533 const unsigned int face,
2534 const unsigned int vertex,
2535 const types::geometric_orientation combined_face_orientation) const
2536{
2537 AssertIndexRange(face, n_faces());
2539 AssertIndexRange(combined_face_orientation, n_face_orientations(face));
2540
2541 switch (this->kind)
2542 {
2544 {
2546 break;
2547 }
2549 {
2550 const auto [face_orientation, face_rotation, face_flip] =
2551 internal::split_face_orientation(combined_face_orientation);
2552
2554 face, vertex, face_orientation, face_flip, face_rotation);
2555 }
2557 {
2558 static constexpr ndarray<unsigned int, 3, 2> table = {
2559 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
2560
2561 return table[face][combined_face_orientation ==
2563 vertex :
2564 (1 - vertex)];
2565 }
2567 {
2568 const auto [face_orientation, face_rotation, face_flip] =
2569 internal::split_face_orientation(combined_face_orientation);
2570
2572 face, vertex, face_orientation, face_flip, face_rotation);
2573 }
2575 {
2576 static constexpr ndarray<unsigned int, 4, 3> table = {
2577 {{{0, 1, 2}}, {{1, 0, 3}}, {{0, 2, 3}}, {{2, 1, 3}}}};
2578
2579 return table[face][standard_to_real_face_vertex(
2580 vertex, face, combined_face_orientation)];
2581 }
2583 {
2584 constexpr auto X = numbers::invalid_unsigned_int;
2585 static constexpr ndarray<unsigned int, 5, 4> table = {
2586 {{{0, 1, 2, 3}},
2587 {{0, 2, 4, X}},
2588 {{3, 1, 4, X}},
2589 {{1, 0, 4, X}},
2590 {{2, 3, 4, X}}}};
2591
2592 return table[face][standard_to_real_face_vertex(
2593 vertex, face, combined_face_orientation)];
2594 }
2596 {
2597 constexpr auto X = numbers::invalid_unsigned_int;
2598 static constexpr ndarray<unsigned int, 6, 4> table = {
2599 {{{1, 0, 2, X}},
2600 {{3, 4, 5, X}},
2601 {{0, 1, 3, 4}},
2602 {{1, 2, 4, 5}},
2603 {{2, 0, 5, 3}}}};
2604
2605 return table[face][standard_to_real_face_vertex(
2606 vertex, face, combined_face_orientation)];
2607 }
2609 {
2610 const auto [face_orientation, face_rotation, face_flip] =
2611 internal::split_face_orientation(combined_face_orientation);
2612
2614 face, vertex, face_orientation, face_flip, face_rotation);
2615 }
2616 default:
2618 }
2619
2621}
2622
2623
2624
2625template <int dim>
2628 const unsigned int vertex) const
2629{
2630 Assert(dim == get_dimension(),
2631 ExcMessage("You can only call this function with arguments for "
2632 "which 'dim' equals the dimension of the space in which "
2633 "the current reference cell lives. You have dim=" +
2634 std::to_string(dim) + " but the reference cell is " +
2635 std::to_string(get_dimension()) + " dimensional."));
2636 return this->template vertex<dim>(face_to_cell_vertices(
2638}
2639
2640
2641
2642inline unsigned int
2644 const unsigned int vertex,
2645 const unsigned int face,
2646 const types::geometric_orientation face_orientation) const
2647{
2648 AssertIndexRange(face, n_faces());
2650
2651 switch (this->kind)
2652 {
2655 break;
2658 ExcMessage(
2659 "In 1D, all faces must have the default orientation."));
2660 return vertex;
2663 return line_vertex_permutations[face_orientation][vertex];
2665 return triangle_vertex_permutations[face_orientation][vertex];
2667 // face 0 is a quadrilateral
2668 if (face == 0)
2669 return quadrilateral_vertex_permutations[face_orientation][vertex];
2670 else
2671 return triangle_vertex_permutations[face_orientation][vertex];
2673 // faces 0 and 1 are triangles
2674 if (face > 1)
2675 return quadrilateral_vertex_permutations[face_orientation][vertex];
2676 else
2677 return triangle_vertex_permutations[face_orientation][vertex];
2679 return quadrilateral_vertex_permutations[face_orientation][vertex];
2680 default:
2682 }
2683
2686}
2687
2688
2689
2690inline unsigned int
2692 const unsigned int line,
2693 const unsigned int face,
2694 const types::geometric_orientation combined_face_orientation) const
2695{
2696 AssertIndexRange(face, n_faces());
2698
2699 static constexpr ndarray<unsigned int, 6, 3> triangle_table = {{{{0, 1, 2}},
2700 {{2, 1, 0}},
2701 {{2, 0, 1}},
2702 {{1, 0, 2}},
2703 {{1, 2, 0}},
2704 {{0, 2, 1}}}};
2705
2706
2707 switch (this->kind)
2708 {
2714 break;
2716 return triangle_table[combined_face_orientation][line];
2718 if (face == 0) // The quadrilateral face
2719 {
2720 const auto [face_orientation, face_rotation, face_flip] =
2721 internal::split_face_orientation(combined_face_orientation);
2722
2724 face_orientation,
2725 face_flip,
2726 face_rotation);
2727 }
2728 else // One of the triangular faces
2729 {
2730 return triangle_table[combined_face_orientation][line];
2731 }
2733 if (face > 1) // One of the quadrilateral faces
2734 {
2735 const auto [face_orientation, face_rotation, face_flip] =
2736 internal::split_face_orientation(combined_face_orientation);
2737
2739 face_orientation,
2740 face_flip,
2741 face_rotation);
2742 }
2743 else // One of the triangular faces
2744 return triangle_table[combined_face_orientation][line];
2746 {
2747 static constexpr ndarray<unsigned int, 8, 4> table = {{
2748 {{0, 1, 2, 3}},
2749 {{2, 3, 0, 1}},
2750 {{3, 2, 0, 1}},
2751 {{0, 1, 3, 2}},
2752 {{1, 0, 3, 2}},
2753 {{3, 2, 1, 0}},
2754 {{2, 3, 1, 0}},
2755 {{1, 0, 2, 3}},
2756 }};
2757 return table[combined_face_orientation][line];
2758 }
2759 default:
2761 }
2762
2764}
2765
2766
2767
2768namespace ReferenceCells
2769{
2770 template <int dim>
2771 inline constexpr const ReferenceCell &
2773 {
2774 switch (dim)
2775 {
2776 case 0:
2778 case 1:
2779 return ReferenceCells::Line;
2780 case 2:
2782 case 3:
2784 default:
2786 }
2788 }
2789
2790
2791
2792 template <int dim>
2793 inline constexpr const ReferenceCell &
2795 {
2796 switch (dim)
2797 {
2798 case 0:
2800 case 1:
2801 return ReferenceCells::Line;
2802 case 2:
2804 case 3:
2806 default:
2808 }
2810 }
2811
2812
2813
2814 template <int structdim>
2815 inline constexpr unsigned int
2820
2821
2822
2823 template <int structdim>
2824 inline constexpr unsigned int
2829
2830
2831
2832 template <int structdim>
2833 inline constexpr unsigned int
2838} // namespace ReferenceCells
2839
2840
2841inline ReferenceCell
2842ReferenceCell::n_vertices_to_type(const int dim, const unsigned int n_vertices)
2843{
2844 AssertIndexRange(dim, 4);
2846
2847 const auto X = ReferenceCells::Invalid;
2848 static const ndarray<ReferenceCell, 4, 9> table = {
2849 {// dim 0
2850 {{X, ReferenceCells::Vertex, X, X, X, X, X, X, X}},
2851 // dim 1
2852 {{X, X, ReferenceCells::Line, X, X, X, X, X, X}},
2853 // dim 2
2854 {{X,
2855 X,
2856 X,
2859 X,
2860 X,
2861 X,
2862 X}},
2863 // dim 3
2864 {{X,
2865 X,
2866 X,
2867 X,
2871 X,
2874 ExcMessage("The combination of dim = " + std::to_string(dim) +
2875 " and n_vertices = " + std::to_string(n_vertices) +
2876 " does not correspond to a known reference cell type."));
2877 return table[dim][n_vertices];
2878}
2879
2880
2881
2882template <int dim>
2883inline double
2885 const unsigned int i) const
2886{
2887 Assert(dim == get_dimension(),
2888 ExcMessage("You can only call this function with arguments for "
2889 "which 'dim' equals the dimension of the space in which "
2890 "the current reference cell lives. You have dim=" +
2891 std::to_string(dim) + " but the reference cell is " +
2892 std::to_string(get_dimension()) + " dimensional."));
2894 switch (this->kind)
2895 {
2897 return 1.0;
2901 if constexpr (dim > 0)
2904 // see also BarycentricPolynomials<2>::compute_value
2906 {
2907 switch (i)
2908 {
2909 case 0:
2910 return 1.0 - xi[std::min(0, dim - 1)] -
2911 xi[std::min(1, dim - 1)];
2912 case 1:
2913 return xi[std::min(0, dim - 1)];
2914 case 2:
2915 return xi[std::min(1, dim - 1)];
2916 default:
2918 }
2919 }
2920 // see also BarycentricPolynomials<3>::compute_value
2922 {
2923 switch (i)
2924 {
2925 case 0:
2926 return 1.0 - xi[std::min(0, dim - 1)] -
2927 xi[std::min(1, dim - 1)] - xi[std::min(2, dim - 1)];
2928 case 1:
2929 return xi[std::min(0, dim - 1)];
2930 case 2:
2931 return xi[std::min(1, dim - 1)];
2932 case 3:
2933 return xi[std::min(2, dim - 1)];
2934 default:
2936 }
2937 }
2938 // see also ScalarLagrangePolynomialPyramid::compute_value()
2940 {
2941 const double Q14 = 0.25;
2942
2943 const double r = xi[std::min(0, dim - 1)];
2944 const double s = xi[std::min(1, dim - 1)];
2945 const double t = xi[std::min(2, dim - 1)];
2946
2947 const double ratio =
2948 (std::fabs(t - 1.0) > 1.0e-14 ? (r * s * t) / (1.0 - t) : 0.0);
2949
2950 if (i == 0)
2951 return Q14 * ((1.0 - r) * (1.0 - s) - t + ratio);
2952 if (i == 1)
2953 return Q14 * ((1.0 + r) * (1.0 - s) - t - ratio);
2954 if (i == 2)
2955 return Q14 * ((1.0 - r) * (1.0 + s) - t - ratio);
2956 if (i == 3)
2957 return Q14 * ((1.0 + r) * (1.0 + s) - t + ratio);
2958 else
2959 return t;
2960 }
2961 // see also ScalarLagrangePolynomialWedge::compute_value()
2964 .d_linear_shape_function<2>(Point<2>(xi[std::min(0, dim - 1)],
2965 xi[std::min(1, dim - 1)]),
2966 i % 3) *
2968 .d_linear_shape_function<1>(Point<1>(xi[std::min(2, dim - 1)]),
2969 i / 3);
2970 default:
2972 }
2973
2975 return 0.0;
2976}
2977
2978
2979
2980template <int dim>
2981inline Tensor<1, dim>
2983 const unsigned int i) const
2984{
2985 Assert(dim == get_dimension(),
2986 ExcMessage("You can only call this function with arguments for "
2987 "which 'dim' equals the dimension of the space in which "
2988 "the current reference cell lives. You have dim=" +
2989 std::to_string(dim) + " but the reference cell is " +
2990 std::to_string(get_dimension()) + " dimensional."));
2991 switch (this->kind)
2992 {
2998 // see also BarycentricPolynomials<2>::compute_grad()
3000 switch (i)
3001 {
3002 case 0:
3003 return Point<dim>(-1.0, -1.0);
3004 case 1:
3005 return Point<dim>(+1.0, +0.0);
3006 case 2:
3007 return Point<dim>(+0.0, +1.0);
3008 default:
3010 }
3011 default:
3013 }
3014
3015 return Point<dim>(+0.0, +0.0, +0.0);
3016}
3017
3018
3019
3020inline double
3022{
3023 switch (this->kind)
3024 {
3026 return 1;
3028 return 1;
3030 return 1. / 2.;
3032 return 1;
3034 return 1. / 6.;
3036 return 4. / 3.;
3038 return 1. / 2.;
3040 return 1;
3041 default:
3043 }
3044
3045 return 0.0;
3046}
3047
3048
3049
3050inline double
3051ReferenceCell::face_measure(const unsigned int face_no) const
3052{
3053 AssertIndexRange(face_no, n_faces());
3054 constexpr double X = std::numeric_limits<double>::signaling_NaN();
3055
3056 static const std::array<double, 5> tri_faces{
3057 {1.0, std::sqrt(2.0), 1.0, X, X}};
3058 static const std::array<double, 5> tet_faces{
3059 {0.5, 0.5, 0.5, std::sqrt(3.0) / 2.0, X}};
3060 static const std::array<double, 5> pyramid_faces{
3061 {4, std::sqrt(2.0), std::sqrt(2.0), std::sqrt(2.0), std::sqrt(2.0)}};
3062 static const std::array<double, 5> wedge_faces{
3063 {0.5, 0.5, 1.0, std::sqrt(2.0), 1.0}};
3064
3065 switch (this->kind)
3066 {
3071 return 1.0;
3073 return tri_faces[face_no];
3075 return tet_faces[face_no];
3077 return pyramid_faces[face_no];
3079 return wedge_faces[face_no];
3080 }
3082 return X;
3083}
3084
3085
3086
3087template <int dim>
3088inline Point<dim>
3090{
3091 Assert(dim == get_dimension(),
3092 ExcMessage("You can only call this function with arguments for "
3093 "which 'dim' equals the dimension of the space in which "
3094 "the current reference cell lives. You have dim=" +
3095 std::to_string(dim) + " but the reference cell is " +
3096 std::to_string(get_dimension()) + " dimensional."));
3097
3098 switch (this->kind)
3099 {
3101 return Point<dim>();
3103 return Point<dim>(1. / 2.);
3105 return Point<dim>(1. / 3., 1. / 3.);
3107 return Point<dim>(1. / 2., 1. / 2.);
3109 return Point<dim>(1. / 4., 1. / 4., 1. / 4.);
3111 return Point<dim>(0, 0, 1. / 4.);
3113 return Point<dim>(1. / 3, 1. / 3, 1. / 2.);
3115 return Point<dim>(1. / 2., 1. / 2., 1. / 2.);
3116 default:
3118 }
3119
3120 return Point<dim>();
3121}
3122
3123
3124
3125template <int dim>
3126inline bool
3127ReferenceCell::contains_point(const Point<dim> &p, const double tolerance) const
3128{
3129 Assert(dim == get_dimension(),
3130 ExcMessage("You can only call this function with arguments for "
3131 "which 'dim' equals the dimension of the space in which "
3132 "the current reference cell lives. You have dim=" +
3133 std::to_string(dim) + " but the reference cell is " +
3134 std::to_string(get_dimension()) + " dimensional."));
3135
3136 // Introduce abbreviations to silence compiler warnings about invalid
3137 // array accesses (that can't happen, of course, but the compiler
3138 // doesn't know that).
3139 constexpr unsigned int x_coordinate = 0;
3140 constexpr unsigned int y_coordinate = (dim >= 2 ? 1 : 0);
3141 constexpr unsigned int z_coordinate = (dim >= 3 ? 2 : 0);
3142
3143 switch (this->kind)
3144 {
3146 {
3147 // Vertices are special cases in that they do not actually
3148 // have coordinates. Error out if this function is called
3149 // with a vertex:
3150 Assert(false,
3151 ExcMessage("Vertices are zero-dimensional objects and "
3152 "as a consequence have no coordinates. You "
3153 "cannot meaningfully ask whether a point is "
3154 "inside a vertex (within a certain tolerance) "
3155 "without coordinate values."));
3156 return false;
3157 }
3161 {
3162 for (unsigned int d = 0; d < dim; ++d)
3163 if ((p[d] < -tolerance) || (p[d] > 1 + tolerance))
3164 return false;
3165 return true;
3166 }
3169 {
3170 // First make sure that we are in the first quadrant or octant
3171 for (unsigned int d = 0; d < dim; ++d)
3172 if (p[d] < -tolerance)
3173 return false;
3174
3175 // Now we also need to make sure that we are below the diagonal line
3176 // or plane that delineates the simplex. This diagonal is given by
3177 // sum(p[d])<=1, and a diagonal a distance eps away is given by
3178 // sum(p[d])<=1+eps*sqrt(d). (For example, the point at (1,1) is a
3179 // distance of 1/sqrt(2) away from the diagonal. That is, its
3180 // sum satisfies
3181 // sum(p[d]) = 2 <= 1 + (1/sqrt(2)) * sqrt(2)
3182 // in other words, it satisfies the predicate with eps=1/sqrt(2).)
3183 double sum = 0;
3184 for (unsigned int d = 0; d < dim; ++d)
3185 sum += p[d];
3186 return (sum <= 1 + tolerance * std::sqrt(1. * dim));
3187 }
3189 {
3190 // A pyramid only lives in the upper half-space:
3191 if (p[z_coordinate] < -tolerance)
3192 return false;
3193
3194 // It also only lives in the space below z=1:
3195 if (p[z_coordinate] > 1 + tolerance)
3196 return false;
3197
3198 // Within what's left of the space, a pyramid is a cone that tapers
3199 // towards the top. First compute the distance of the point to the
3200 // axis in the max norm (this is the right norm because the vertices
3201 // of the pyramid are at points +/-1, +/-1):
3202 const double distance_from_axis =
3203 std::max(std::fabs(p[x_coordinate]), std::fabs(p[y_coordinate]));
3204
3205 // We are inside the pyramid if the distance from the axis is less
3206 // than (1-z)
3207 return (distance_from_axis <= 1 + tolerance - p[z_coordinate]);
3208 }
3210 {
3211 // The wedge we use is a triangle extruded into the third
3212 // dimension by one unit. So we can use the same logic as for
3213 // triangles above (i.e., for the simplex above, using dim==2)
3214 // and then check the third dimension separately.
3215
3216 if ((p[x_coordinate] < -tolerance) || (p[y_coordinate] < -tolerance))
3217 return false;
3218
3219 const double sum = p[x_coordinate] + p[y_coordinate];
3220 if (sum > 1 + tolerance * std::sqrt(2.0))
3221 return false;
3222
3223 if (p[z_coordinate] < -tolerance)
3224 return false;
3225 if (p[z_coordinate] > 1 + tolerance)
3226 return false;
3227
3228 return true;
3229 }
3230 default:
3232 }
3233
3234 return false;
3235}
3236
3237
3238
3239template <int dim>
3240inline Tensor<1, dim>
3242 const unsigned int i) const
3243{
3244 return face_tangent_vector<dim>(face_no, i);
3245}
3246
3247
3248
3249template <int dim>
3250inline Tensor<1, dim>
3251ReferenceCell::face_tangent_vector(const unsigned int face_no,
3252 const unsigned int i) const
3253{
3254 Assert(dim == get_dimension(),
3255 ExcMessage("You can only call this function with arguments for "
3256 "which 'dim' equals the dimension of the space in which "
3257 "the current reference cell lives. You have dim=" +
3258 std::to_string(dim) + " but the reference cell is " +
3259 std::to_string(get_dimension()) + " dimensional."));
3260 AssertIndexRange(i, dim - 1);
3261
3262 switch (this->kind)
3263 {
3271 {
3272 AssertIndexRange(face_no, 3);
3273 static const std::array<Tensor<1, dim>, 3> table = {
3274 {Point<dim>(1, 0),
3275 Point<dim>(-std::sqrt(0.5), +std::sqrt(0.5)),
3276 Point<dim>(0, -1)}};
3277
3278 return table[face_no];
3279 }
3281 {
3282 AssertIndexRange(face_no, 4);
3283 static const ndarray<Tensor<1, dim>, 4, 2> table = {
3284 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
3285 {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
3286 {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
3287 {{Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
3288 +std::pow(1.0 / 3.0, 1.0 / 4.0),
3289 0),
3290 Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
3291 0,
3292 +std::pow(1.0 / 3.0, 1.0 / 4.0))}}}};
3293
3294 return table[face_no][i];
3295 }
3297 {
3298 AssertIndexRange(face_no, 5);
3299 static const ndarray<Tensor<1, dim>, 5, 2> table = {
3300 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
3301 {{Point<dim>(+1.0 / std::sqrt(2.0), 0, +1.0 / std::sqrt(2.0)),
3302 Point<dim>(0, 1, 0)}},
3303 {{Point<dim>(+1.0 / std::sqrt(2.0), 0, -1.0 / std::sqrt(2.0)),
3304 Point<dim>(0, 1, 0)}},
3305 {{Point<dim>(1, 0, 0),
3306 Point<dim>(0, +1.0 / std::sqrt(2.0), +1.0 / std::sqrt(2.0))}},
3307 {{Point<dim>(1, 0, 0),
3308 Point<dim>(0, +1.0 / std::sqrt(2.0), -1.0 / std::sqrt(2.0))}}}};
3309
3310 return table[face_no][i];
3311 }
3313 {
3314 AssertIndexRange(face_no, 5);
3315 static const ndarray<Tensor<1, dim>, 5, 2> table = {
3316 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
3317 {{Point<dim>(1, 0, 0), Point<dim>(0, 1, 0)}},
3318 {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
3319 {{Point<dim>(-1 / std::sqrt(2.0), +1 / std::sqrt(2.0), 0),
3320 Point<dim>(0, 0, 1)}},
3321 {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}}}};
3322
3323 return table[face_no][i];
3324 }
3325 default:
3327 }
3328
3329
3330 return {};
3331}
3332
3333
3334
3335template <int dim>
3336inline Tensor<1, dim>
3337ReferenceCell::unit_normal_vectors(const unsigned int face_no) const
3338{
3339 return face_normal_vector<dim>(face_no);
3340}
3341
3342
3343
3344template <int dim>
3345inline Tensor<1, dim>
3346ReferenceCell::face_normal_vector(const unsigned int face_no) const
3347{
3348 Assert(dim == get_dimension(),
3349 ExcMessage("You can only call this function with arguments for "
3350 "which 'dim' equals the dimension of the space in which "
3351 "the current reference cell lives. You have dim=" +
3352 std::to_string(dim) + " but the reference cell is " +
3353 std::to_string(get_dimension()) + " dimensional."));
3354
3355 if (is_hyper_cube())
3356 {
3359 }
3360 else if (dim == 2)
3361 {
3363
3364 // Return the rotated vector
3365 return cross_product_2d(face_tangent_vector<dim>(face_no, 0));
3366 }
3367 else if (dim == 3)
3368 {
3369 return cross_product_3d(face_tangent_vector<dim>(face_no, 0),
3370 face_tangent_vector<dim>(face_no, 1));
3371 }
3372
3374
3375 return {};
3376}
3377
3378
3379
3380inline unsigned int
3381ReferenceCell::n_face_orientations(const unsigned int face_no) const
3382{
3383 AssertIndexRange(face_no, n_faces());
3384 if (get_dimension() == 1)
3385 return 1;
3386 if (get_dimension() == 2)
3387 return 2;
3389 return 8;
3390 else if (face_reference_cell(face_no) == ReferenceCells::Triangle)
3391 return 6;
3392
3395}
3396
3397
3398
3401 const unsigned int line,
3402 const unsigned int face,
3403 const types::geometric_orientation combined_face_orientation,
3404 const types::geometric_orientation line_orientation) const
3405{
3407 face,
3408 combined_face_orientation,
3409 line_orientation);
3410}
3411
3412
3415 const unsigned int face_line_no,
3416 const unsigned int face_no,
3417 const types::geometric_orientation combined_face_orientation,
3418 const types::geometric_orientation line_orientation) const
3419{
3420 constexpr auto D = numbers::default_geometric_orientation;
3421 constexpr auto R = numbers::reverse_line_orientation;
3422 if (*this == ReferenceCells::Hexahedron)
3423 {
3424 static constexpr ::ndarray<types::geometric_orientation, 2, 8>
3425 table{{{{D, D, D, R, R, R, R, D}}, {{D, D, R, D, R, R, D, R}}}};
3426 // We use face_line_no / 2 here since lines i and i + 1 are parallel and,
3427 // on a given face, have the same relative orientations.
3428 const bool match =
3429 line_orientation == table[face_line_no / 2][combined_face_orientation];
3430
3433 }
3434 else if (*this == ReferenceCells::Tetrahedron)
3435 {
3436 static constexpr unsigned int X = numbers::invalid_unsigned_int;
3437 static constexpr ::ndarray<unsigned int, 4, 3> combined_lines{
3438 {{{0, 0, 0}}, {{X, 0, 1}}, {{X, 0, X}}, {{X, X, X}}}};
3439
3440 const auto combined_line = combined_lines[face_no][face_line_no];
3441
3442 Assert(combined_line != X,
3443 ExcMessage(
3444 "This function can only be called for following face-line "
3445 "combinations: (0,0), (0,1), (0,2), (1,1), (1,2), (2,1),"));
3446
3447 static constexpr ::ndarray<types::geometric_orientation, 2, 6>
3448 table{{{{D, R, D, R, D, R}}, {{R, D, R, D, R, D}}}};
3449
3450 const bool match =
3451 line_orientation == table[combined_line][combined_face_orientation];
3452
3455 }
3456 else
3457 // TODO: This might actually be wrong for some of the other
3458 // kinds of objects. We should check this
3460}
3461
3462
3463
3464namespace internal
3465{
3466 template <typename T>
3468 {
3469 public:
3483
3487 virtual ~NoPermutation() noexcept override = default;
3488
3492 virtual void
3493 print_info(std::ostream &out) const override
3494 {
3495 out << '[';
3496
3497 const unsigned int n_vertices = entity_type.n_vertices();
3498
3499 for (unsigned int i = 0; i < n_vertices; ++i)
3500 {
3501 out << vertices_0[i];
3502 if (i + 1 != n_vertices)
3503 out << ',';
3504 }
3505
3506 out << "] is not a valid permutation of [";
3507
3508 for (unsigned int i = 0; i < n_vertices; ++i)
3509 {
3510 out << vertices_1[i];
3511 if (i + 1 != n_vertices)
3512 out << ',';
3513 }
3514
3515 out << "]." << std::endl;
3516 }
3517
3522
3527
3532 };
3533
3545 << "The reference-cell type used on this cell (" << arg1.to_string()
3546 << ") does not match the reference-cell type of the finite element "
3547 << "associated with this cell (" << arg2.to_string() << "). "
3548 << "Did you accidentally use simplex elements on hypercube meshes "
3549 << "(or the other way around), or are you using a mixed mesh and "
3550 << "assigned a simplex element to a hypercube cell (or the other "
3551 << "way around) via the active_fe_index?");
3552} // namespace internal
3553
3554
3555
3556template <typename T, std::size_t N>
3558ReferenceCell::compute_orientation(const std::array<T, N> &vertices_0,
3559 const std::array<T, N> &vertices_1) const
3560{
3561 Assert(N >= n_vertices(),
3562 ExcMessage("The number of array elements must be equal to or "
3563 "greater than the number of vertices of the cell "
3564 "referenced by this object."));
3565
3566 // Call the non-deprecated function, taking care of calling it only with
3567 // those array elements that we actually care about (see the note
3568 // in the documentation about the arguments potentially being
3569 // larger arrays than necessary).
3571 make_array_view(vertices_0.begin(), vertices_0.begin() + n_vertices()),
3572 make_array_view(vertices_1.begin(), vertices_1.begin() + n_vertices()));
3573}
3574
3575
3576
3577template <typename T>
3580 const ArrayView<const T> &vertices_0,
3581 const ArrayView<const T> &vertices_1) const
3582{
3583 Assert(vertices_0.size() == n_vertices(),
3584 ExcMessage("The number of array elements must be equal to "
3585 "the number of vertices of the cell "
3586 "referenced by this object."));
3587 Assert(vertices_1.size() == n_vertices(),
3588 ExcMessage("The number of array elements must be equal to "
3589 "the number of vertices of the cell "
3590 "referenced by this object."));
3591
3592 auto compute_orientation = [&](const auto &table) {
3593 for (types::geometric_orientation o = 0; o < table.size(); ++o)
3594 {
3595 bool match = true;
3596 for (unsigned int j = 0; j < table[o].size(); ++j)
3597 match = (match && vertices_0[j] == vertices_1[table[o][j]]);
3598
3599 if (match)
3600 return o;
3601 }
3602
3603 // Do not use `this` in Assert because nvcc when using C++20 assumes that
3604 // `this` is an integer and we get the following error: invalid type
3605 // argument of unary '*' (have 'int')
3606 [[maybe_unused]] const auto &ref_cell = *this;
3607 Assert(false,
3608 (internal::NoPermutation<T>(ref_cell, vertices_0, vertices_1)));
3609 return std::numeric_limits<types::geometric_orientation>::max();
3610 };
3611
3612 switch (this->kind)
3613 {
3622 default:
3624 }
3625
3626 Assert(false, (internal::NoPermutation<T>(*this, vertices_0, vertices_1)));
3627 return std::numeric_limits<types::geometric_orientation>::max();
3628}
3629
3630
3631
3632template <typename T, std::size_t N>
3633inline std::array<T, N>
3635 const std::array<T, N> &vertices,
3636 const unsigned int orientation) const
3637{
3638 Assert(N >= n_vertices(),
3639 ExcMessage("The number of array elements must be equal to or "
3640 "greater than the number of vertices of the cell "
3641 "referenced by this object."));
3642
3643 // Call the non-deprecated function, taking care of calling it only with
3644 // those array elements that we actually care about (see the note
3645 // in the documentation about the arguments potentially being
3646 // larger arrays than necessary).
3647 const auto permutation = permute_by_combined_orientation(
3648 make_array_view(vertices.begin(), vertices.begin() + n_vertices()),
3649 orientation);
3650
3651 std::array<T, N> temp;
3652 std::copy(permutation.begin(), permutation.end(), temp.begin());
3653
3654 return temp;
3655}
3656
3657
3658
3659template <typename T>
3660boost::container::small_vector<T, 8>
3662 const ArrayView<const T> &vertices,
3663 const types::geometric_orientation orientation) const
3664{
3665 AssertDimension(vertices.size(), n_vertices());
3666 boost::container::small_vector<T, 8> result(n_vertices());
3667
3668 auto permute = [&](const auto &table) {
3669 AssertIndexRange(orientation, table.size());
3670 for (unsigned int j = 0; j < table[orientation].size(); ++j)
3671 result[j] = vertices[table[orientation][j]];
3672 };
3673
3674 switch (this->kind)
3675 {
3678 break;
3680 permute(line_vertex_permutations);
3681 break;
3684 break;
3687 break;
3688 default:
3690 }
3691
3692 return result;
3693}
3694
3695
3696
3699 const types::geometric_orientation orientation) const
3700{
3701 switch (this->kind)
3702 {
3704 // Things are always default-oriented in 1D
3705 return orientation;
3706
3708 // the 1d orientations are the identity and a flip: i.e., the identity
3709 // and an involutory mapping
3710 return orientation;
3711
3713 {
3714 AssertIndexRange(orientation, 6);
3715 constexpr std::array<types::geometric_orientation, 6> inverses{
3716 {0, 1, 4, 3, 2, 5}};
3717 return inverses[orientation];
3718 }
3720 {
3721 AssertIndexRange(orientation, 8);
3722 constexpr std::array<types::geometric_orientation, 8> inverses{
3723 {0, 1, 6, 3, 4, 5, 2, 7}};
3724 return inverses[orientation];
3725 }
3726 default:
3728 }
3729
3730 return std::numeric_limits<types::geometric_orientation>::max();
3731}
3732
3733
3734
3735template <>
3736unsigned int
3737ReferenceCell::vtk_lexicographic_to_node_index<0>(
3738 const std::array<unsigned, 0> &node_indices,
3739 const std::array<unsigned, 0> &nodes_per_direction,
3740 const bool legacy_format) const;
3741
3742template <>
3743unsigned int
3744ReferenceCell::vtk_lexicographic_to_node_index<1>(
3745 const std::array<unsigned, 1> &node_indices,
3746 const std::array<unsigned, 1> &nodes_per_direction,
3747 const bool legacy_format) const;
3748
3749template <>
3750unsigned int
3751ReferenceCell::vtk_lexicographic_to_node_index<2>(
3752 const std::array<unsigned, 2> &node_indices,
3753 const std::array<unsigned, 2> &nodes_per_direction,
3754 const bool legacy_format) const;
3755
3756template <>
3757unsigned int
3758ReferenceCell::vtk_lexicographic_to_node_index<3>(
3759 const std::array<unsigned, 3> &node_indices,
3760 const std::array<unsigned, 3> &nodes_per_direction,
3761 const bool legacy_format) const;
3762
3764
3765#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:954
std::size_t size() const
Definition array_view.h:689
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
double face_measure(const unsigned int face_no) 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:281
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#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 ReferenceCell Triangle
constexpr ReferenceCell Hexahedron
constexpr ReferenceCell Pyramid
constexpr ReferenceCell Wedge
constexpr ReferenceCell Vertex
constexpr unsigned int max_n_lines()
constexpr const ReferenceCell & get_hypercube()
constexpr ReferenceCell Invalid
constexpr unsigned int max_n_vertices()
constexpr unsigned int max_n_faces()
constexpr ReferenceCell Quadrilateral
constexpr ReferenceCell Tetrahedron
constexpr const ReferenceCell & get_simplex()
constexpr ReferenceCell Line
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:238
constexpr types::geometric_orientation reverse_line_orientation
Definition types.h:365
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:352
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 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)