deal.II version GIT relicensing-3540-g7552a02177 2025-06-20 13:50: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
418 template <int dim>
419 std::pair<unsigned int, RefinementCase<dim - 1>>
421 const types::geometric_orientation combined_face_orientation,
422 const internal::SubfaceCase<dim> subface_case,
423 const unsigned int subface_no) const;
424
438 face_reference_cell(const unsigned int face_no) const;
439
455 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
456 "Use numbers::default_geometric_orientation instead.")
457 static constexpr types::geometric_orientation
459
468 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
469 "Use numbers::reverse_line_orientation instead.")
470 static constexpr types::geometric_orientation
472
480 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
481 "Use the version of this function which takes the orientation as the third "
482 "parameter instead.")
483 unsigned int
484 child_cell_on_face(const unsigned int face, const unsigned int subface) const;
485
526 unsigned int
528 const unsigned int face,
529 const unsigned int subface,
530 const types::geometric_orientation combined_orientation) const;
531
541 std::array<unsigned int, 2>
542 standard_vertex_to_face_and_vertex_index(const unsigned int vertex) const;
543
553 std::array<unsigned int, 2>
554 standard_line_to_face_and_line_index(const unsigned int line) const;
555
564 unsigned int
565 line_to_cell_vertices(const unsigned int line,
566 const unsigned int vertex) const;
567
571 unsigned int
572 face_to_cell_lines(const unsigned int face,
573 const unsigned int line,
574 const types::geometric_orientation face_orientation) const;
575
579 unsigned int
581 const unsigned int face,
582 const unsigned int vertex,
583 const types::geometric_orientation face_orientation) const;
584
608 template <int dim>
610 face_vertex_location(const unsigned int face,
611 const unsigned int vertex) const;
612
616 unsigned int
618 const unsigned int vertex,
619 const unsigned int face,
620 const types::geometric_orientation face_orientation) const;
621
625 unsigned int
627 const unsigned int line,
628 const unsigned int face,
629 const types::geometric_orientation face_orientation) const;
630
638 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
639 "Use face_to_cell_line_orientation() instead.")
642 const unsigned int line,
643 const unsigned int face,
644 const types::geometric_orientation face_orientation,
645 const types::geometric_orientation line_orientation) const;
646
671 const unsigned int face_line_no,
672 const unsigned int face_no,
673 const types::geometric_orientation face_orientation,
674 const types::geometric_orientation line_orientation) const;
675
699 double
700 volume() const;
701
717 double
718 face_measure(const unsigned int face_no) const;
719
736 template <int dim>
738 barycenter() const;
739
763 template <int dim>
764 bool
765 contains_point(const Point<dim> &p, const double tolerance = 0) const;
766
775 template <int dim>
777 closest_point(const Point<dim> &p) const;
778
790 template <int dim>
792 face_tangent_vector(const unsigned int face_no, const unsigned int i) const;
793
807 template <int dim>
808 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT("Use face_tangent_vector() instead.")
809 Tensor<1, dim> unit_tangential_vectors(const unsigned int face_no,
810 const unsigned int i) const;
811
819 template <int dim>
821 face_normal_vector(const unsigned int face_no) const;
822
832 template <int dim>
833 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT("Use face_normal_vector() instead.")
834 Tensor<1, dim> unit_normal_vectors(const unsigned int face_no) const;
835
841 unsigned int
842 n_face_orientations(const unsigned int face_no) const;
843
860 template <typename T, std::size_t N>
862 compute_orientation(const std::array<T, N> &vertices_0,
863 const std::array<T, N> &vertices_1) const;
864
889 template <typename T>
892 const ArrayView<const T> &vertices_1) const;
893
894
908 template <typename T, std::size_t N>
909 DEAL_II_DEPRECATED std::array<T, N>
910 permute_according_orientation(const std::array<T, N> &vertices,
911 const unsigned int orientation) const;
912
924 template <typename T>
925 boost::container::small_vector<T, 8>
927 const ArrayView<const T> &vertices,
928 const types::geometric_orientation orientation) const;
929
937 const types::geometric_orientation orientation) const;
938
943 faces_for_given_vertex(const unsigned int vertex_index) const;
944
949 constexpr ::ndarray<unsigned int, 12, 4>
950 new_isotropic_child_face_lines(const unsigned int refinement_choice) const;
951
956 constexpr ::ndarray<unsigned int, 12, 4, 2>
958 const unsigned int refinement_choice) const;
959
964 constexpr ::ndarray<unsigned int, 8, 6>
965 new_isotropic_child_cell_faces(const unsigned int refinement_choice) const;
966
971 constexpr ::ndarray<unsigned int, 8, 4>
972 new_isotropic_child_cell_vertices(const unsigned int refinement_choice) const;
973
986 unsigned int
987 exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const;
988
992 unsigned int
993 exodusii_face_to_deal_face(const unsigned int face_n) const;
994
998 unsigned int
999 unv_vertex_to_deal_vertex(const unsigned int vertex_n) const;
1000
1004 unsigned int
1005 vtk_linear_type() const;
1006
1011 unsigned int
1012 vtk_quadratic_type() const;
1013
1018 unsigned int
1019 vtk_lagrange_type() const;
1020
1037 template <int dim>
1038 unsigned int
1040 const std::array<unsigned, dim> &node_indices,
1041 const std::array<unsigned, dim> &nodes_per_direction,
1042 const bool legacy_format) const;
1043
1047 unsigned int
1048 vtk_vertex_to_deal_vertex(const unsigned int vertex_index) const;
1049
1053 unsigned int
1054 gmsh_element_type() const;
1055
1069 std::string
1070 to_string() const;
1071
1075 constexpr operator std::uint8_t() const;
1076
1080 constexpr bool
1081 operator==(const ReferenceCell &type) const;
1082
1086 constexpr bool
1087 operator!=(const ReferenceCell &type) const;
1088
1094 template <class Archive>
1095 void
1096 serialize(Archive &archive, const unsigned int /*version*/);
1097
1101 static constexpr std::size_t
1103
1108private:
1112 std::uint8_t kind;
1113
1120 constexpr ReferenceCell(const std::uint8_t kind);
1121
1128 {{{0}}}};
1129
1134 {{{0, 1}}, {{1, 0}}}};
1135
1140 {{{0, 1, 2}},
1141 {{0, 2, 1}},
1142 {{2, 0, 1}},
1143 {{2, 1, 0}},
1144 {{1, 2, 0}},
1145 {{1, 0, 2}}}};
1146
1150 static constexpr ndarray<unsigned int, 8, 4>
1152 {{0, 1, 2, 3}},
1153 {{0, 2, 1, 3}},
1154 {{2, 0, 3, 1}},
1155 {{2, 3, 0, 1}},
1156 {{3, 2, 1, 0}},
1157 {{3, 1, 2, 0}},
1158 {{1, 3, 0, 2}},
1159 {{1, 0, 3, 2}},
1160 }};
1161
1166 friend constexpr ReferenceCell
1168
1169 friend std::ostream &
1170 operator<<(std::ostream &out, const ReferenceCell &reference_cell);
1171
1172 friend std::istream &
1173 operator>>(std::istream &in, ReferenceCell &reference_cell);
1174};
1175
1176
1184std::ostream &
1185operator<<(std::ostream &out, const ReferenceCell &reference_cell);
1186
1194std::istream &
1195operator>>(std::istream &in, ReferenceCell &reference_cell);
1196
1197
1198inline constexpr ReferenceCell::ReferenceCell(const std::uint8_t kind)
1199 : kind(kind)
1200{}
1201
1202
1203
1204inline constexpr ReferenceCell::operator std::uint8_t() const
1205{
1206 return kind;
1207}
1208
1209
1210
1211inline constexpr bool
1213{
1214 return kind == type.kind;
1215}
1216
1217
1218
1219inline constexpr bool
1221{
1222 return kind != type.kind;
1223}
1224
1225
1226
1227namespace internal
1228{
1229 inline constexpr ReferenceCell
1230 make_reference_cell_from_int(const std::uint8_t kind)
1231 {
1232#ifndef DEAL_II_CXX14_CONSTEXPR_BUG
1233 // Make sure these are the only indices from which objects can be
1234 // created.
1235 Assert((kind == std::numeric_limits<std::uint8_t>::max()) || (kind < 8),
1237#endif
1238
1239 // Call the private constructor, which we can from here because this
1240 // function is a 'friend'.
1241 return {kind};
1242 }
1243} // namespace internal
1244
1245
1246
1255{
1268 std::numeric_limits<std::uint8_t>::max());
1269
1275 template <int dim>
1276 constexpr const ReferenceCell &
1277 get_simplex();
1278
1284 template <int dim>
1285 constexpr const ReferenceCell &
1286 get_hypercube();
1287
1295 template <int structdim>
1296 constexpr unsigned int
1298
1306 template <int structdim>
1307 constexpr unsigned int
1308 max_n_lines();
1309
1325 template <int structdim>
1326 constexpr unsigned int
1327 max_n_faces();
1328} // namespace ReferenceCells
1329
1330
1331
1333 : ReferenceCell(ReferenceCells::Invalid)
1334{}
1335
1336
1337
1338template <class Archive>
1339inline void
1340ReferenceCell::serialize(Archive &archive, const unsigned int /*version*/)
1341{
1342 archive &kind;
1343}
1344
1345
1346
1347inline constexpr std::size_t
1352
1353
1354
1356ReferenceCell::faces_for_given_vertex(const unsigned int vertex) const
1357{
1359 switch (this->kind)
1360 {
1362 return {&GeometryInfo<2>::vertex_to_face[vertex][0], 1};
1364 return {&GeometryInfo<2>::vertex_to_face[vertex][0], 2};
1366 {
1367 static constexpr ndarray<unsigned int, 3, 2> table = {
1368 {{{0, 2}}, {{0, 1}}, {{1, 2}}}};
1369 return table[vertex];
1370 }
1372 {
1373 static constexpr ndarray<unsigned int, 4, 3> table = {
1374 {{{0, 1, 2}}, {{0, 1, 3}}, {{0, 2, 3}}, {{1, 2, 3}}}};
1375
1376 return table[vertex];
1377 }
1379 {
1380 static constexpr unsigned int X = numbers::invalid_unsigned_int;
1381 static constexpr ndarray<unsigned int, 5, 4> table = {
1382 {{{0, 1, 3, X}},
1383 {{0, 2, 3, X}},
1384 {{0, 1, 4, X}},
1385 {{0, 2, 4, X}},
1386 {{1, 2, 3, 4}}}};
1387
1388 return {&table[vertex][0], vertex == 4 ? 4u : 3u};
1389 }
1391 {
1393 static constexpr ndarray<unsigned int, 6, 3> table = {{{{0, 2, 4}},
1394 {{0, 2, 3}},
1395 {{0, 3, 4}},
1396 {{1, 2, 4}},
1397 {{1, 2, 3}},
1398 {{1, 3, 4}}}};
1399
1400 return table[vertex];
1401 }
1403 return {&GeometryInfo<3>::vertex_to_face[vertex][0], 3};
1404 default:
1406 }
1407
1408 return {};
1409}
1410
1411
1412constexpr ::ndarray<unsigned int, 12, 4>
1414 const unsigned int refienement_choice) const
1415{
1416 // AssertIndexRange(refienement_choice, n_isotropic_refinement_choices());
1417 const unsigned int X = numbers::invalid_unsigned_int;
1418 switch (this->kind)
1419 {
1421 {
1422 constexpr ::ndarray<unsigned int, 3, 12, 4> new_quad_lines_tet =
1423 {{// new line is (6,8)
1424 {{{{2, 3, 8, X}},
1425 {{0, 9, 5, X}},
1426 {{1, 6, 11, X}},
1427 {{4, 10, 7, X}},
1428 {{2, 12, 5, X}},
1429 {{1, 9, 12, X}},
1430 {{4, 8, 12, X}},
1431 {{6, 12, 10, X}},
1432 {{X, X, X, X}},
1433 {{X, X, X, X}},
1434 {{X, X, X, X}},
1435 {{X, X, X, X}}}},
1436 // new line is (5,7)
1437 {{{{2, 3, 8, X}},
1438 {{0, 9, 5, X}},
1439 {{1, 6, 11, X}},
1440 {{4, 10, 7, X}},
1441 {{0, 3, 12, X}},
1442 {{1, 12, 8, X}},
1443 {{4, 12, 9, X}},
1444 {{7, 11, 12, X}},
1445 {{X, X, X, X}},
1446 {{X, X, X, X}},
1447 {{X, X, X, X}},
1448 {{X, X, X, X}}}},
1449 // new line is (4,9)
1450 {{{{2, 3, 8, X}},
1451 {{0, 9, 5, X}},
1452 {{1, 6, 11, X}},
1453 {{4, 10, 7, X}},
1454 {{0, 12, 11, X}},
1455 {{2, 6, 12, X}},
1456 {{3, 12, 7, X}},
1457 {{5, 10, 12, X}},
1458 {{X, X, X, X}},
1459 {{X, X, X, X}},
1460 {{X, X, X, X}},
1461 {{X, X, X, X}}}}}};
1462 return new_quad_lines_tet[refienement_choice];
1463 }
1464
1466 {
1467 constexpr ::ndarray<unsigned int, 12, 4> new_quad_lines_hex = {
1468 {{{10, 28, 16, 24}},
1469 {{28, 14, 17, 25}},
1470 {{11, 29, 24, 20}},
1471 {{29, 15, 25, 21}},
1472 {{18, 26, 0, 28}},
1473 {{26, 22, 1, 29}},
1474 {{19, 27, 28, 4}},
1475 {{27, 23, 29, 5}},
1476 {{2, 24, 8, 26}},
1477 {{24, 6, 9, 27}},
1478 {{3, 25, 26, 12}},
1479 {{25, 7, 27, 13}}}};
1480 return new_quad_lines_hex;
1481 }
1482 default:
1484 }
1485
1486 return {};
1487}
1488
1489
1490
1491constexpr ::ndarray<unsigned int, 12, 4, 2>
1493 const unsigned int refienement_choice) const
1494{
1495 // AssertIndexRange(refienement_choice, n_isotropic_refinement_choices());
1496 const unsigned int X = numbers::invalid_unsigned_int;
1497 switch (this->kind)
1498 {
1500 {
1501 constexpr ::ndarray<unsigned int, 3, 12, 4, 2> table_tet = {
1502 {// new line is (6, 8)
1503 {{{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
1504 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
1505 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
1506 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
1507 {{{{4, 6}}, {{6, 8}}, {{8, 4}}, {{X, X}}}},
1508 {{{{6, 5}}, {{5, 8}}, {{8, 6}}, {{X, X}}}},
1509 {{{{8, 7}}, {{7, 6}}, {{6, 8}}, {{X, X}}}},
1510 {{{{9, 6}}, {{6, 8}}, {{8, 9}}, {{X, X}}}},
1511 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1512 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1513 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1514 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}},
1515 // new line is (5, 7)
1516 {{{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
1517 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
1518 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
1519 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
1520 {{{{5, 4}}, {{4, 7}}, {{7, 5}}, {{X, X}}}},
1521 {{{{6, 5}}, {{5, 7}}, {{7, 6}}, {{X, X}}}},
1522 {{{{8, 7}}, {{7, 5}}, {{5, 8}}, {{X, X}}}},
1523 {{{{7, 9}}, {{9, 5}}, {{5, 7}}, {{X, X}}}},
1524 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1525 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1526 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1527 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}},
1528 // new line is (4, 9)
1529 {{{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}},
1530 {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}},
1531 {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}},
1532 {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}},
1533 {{{{5, 4}}, {{4, 9}}, {{9, 5}}, {{X, X}}}},
1534 {{{{4, 6}}, {{6, 9}}, {{9, 4}}, {{X, X}}}},
1535 {{{{7, 4}}, {{4, 9}}, {{9, 7}}, {{X, X}}}},
1536 {{{{4, 8}}, {{8, 9}}, {{9, 4}}, {{X, X}}}},
1537 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1538 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1539 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}},
1540 {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}}}};
1541 return table_tet[refienement_choice];
1542 }
1543
1545 {
1546 constexpr ::ndarray<unsigned int, 12, 4, 2> table_hex = {
1547 {{{{{10, 22}}, {{24, 26}}, {{10, 24}}, {{22, 26}}}},
1548 {{{{24, 26}}, {{11, 23}}, {{24, 11}}, {{26, 23}}}},
1549 {{{{22, 14}}, {{26, 25}}, {{22, 26}}, {{14, 25}}}},
1550 {{{{26, 25}}, {{23, 15}}, {{26, 23}}, {{25, 15}}}},
1551 {{{{8, 24}}, {{20, 26}}, {{8, 20}}, {{24, 26}}}},
1552 {{{{20, 26}}, {{12, 25}}, {{20, 12}}, {{26, 25}}}},
1553 {{{{24, 9}}, {{26, 21}}, {{24, 26}}, {{9, 21}}}},
1554 {{{{26, 21}}, {{25, 13}}, {{26, 25}}, {{21, 13}}}},
1555 {{{{16, 20}}, {{22, 26}}, {{16, 22}}, {{20, 26}}}},
1556 {{{{22, 26}}, {{17, 21}}, {{22, 17}}, {{26, 21}}}},
1557 {{{{20, 18}}, {{26, 23}}, {{20, 26}}, {{18, 23}}}},
1558 {{{{26, 23}}, {{21, 19}}, {{26, 21}}, {{23, 19}}}}}};
1559 return table_hex;
1560 }
1561 default:
1563 }
1564
1565 return {};
1566}
1567
1568
1569
1570constexpr ::ndarray<unsigned int, 8, 6>
1572 const unsigned int refienement_choice) const
1573{
1574 // AssertIndexRange(refienement_choice, n_isotropic_refinement_choices());
1575 const unsigned int X = numbers::invalid_unsigned_int;
1576 switch (this->kind)
1577 {
1579 {
1580 constexpr ::ndarray<unsigned int, 3, 8, 6> cell_quads_tet = {
1581 {// new line is (6, 8)
1582 {{{{8, 13, 16, 0, X, X}},
1583 {{9, 12, 1, 21, X, X}},
1584 {{10, 2, 17, 20, X, X}},
1585 {{3, 14, 18, 22, X, X}},
1586 {{11, 1, 4, 5, X, X}},
1587 {{15, 0, 4, 6, X, X}},
1588 {{19, 7, 6, 3, X, X}},
1589 {{23, 5, 2, 7, X, X}}}},
1590 // new line is (5, 7)
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, 0, 5, X, X}},
1597 {{15, 4, 1, 6, X, X}},
1598 {{19, 2, 5, 7, X, X}},
1599 {{23, 6, 7, 3, X, X}},
1600 }},
1601 // new line is (4, 9)
1602 {{
1603 {{8, 13, 16, 0, X, X}},
1604 {{9, 12, 1, 21, X, X}},
1605 {{10, 2, 17, 20, X, X}},
1606 {{3, 14, 18, 22, X, X}},
1607 {{11, 4, 5, 2, X, X}},
1608 {{15, 6, 7, 3, X, X}},
1609 {{19, 5, 0, 6, X, X}},
1610 {{23, 1, 4, 7, X, X}},
1611 }}}};
1612 return cell_quads_tet[refienement_choice];
1613 }
1614
1616 {
1617 constexpr ::ndarray<unsigned int, 8, 6> cell_quads_hex = {{
1618 {{12, 0, 20, 4, 28, 8}}, // bottom children
1619 {{0, 16, 22, 6, 29, 9}}, //
1620 {{13, 1, 4, 24, 30, 10}}, //
1621 {{1, 17, 6, 26, 31, 11}}, //
1622 {{14, 2, 21, 5, 8, 32}}, // top children
1623 {{2, 18, 23, 7, 9, 33}}, //
1624 {{15, 3, 5, 25, 10, 34}}, //
1625 {{3, 19, 7, 27, 11, 35}} //
1626 }};
1627 return cell_quads_hex;
1628 }
1629 default:
1631 }
1632
1633 return {};
1634}
1635
1636
1637
1638constexpr ::ndarray<unsigned int, 8, 4>
1640 const unsigned int refienement_choice) const
1641{
1642 // AssertIndexRange(refienement_choice, n_isotropic_refinement_choices());
1643
1644 switch (this->kind)
1645 {
1647 {
1648 constexpr ::ndarray<unsigned int, 3, 8, 4> cell_vertices_tet = {
1649 {// new line is (6,8)
1650 {{{{0, 4, 6, 7}},
1651 {{4, 1, 5, 8}},
1652 {{6, 5, 2, 9}},
1653 {{7, 8, 9, 3}},
1654 {{4, 5, 6, 8}},
1655 {{4, 7, 8, 6}},
1656 {{6, 9, 7, 8}},
1657 {{5, 8, 9, 6}}}},
1658 // new line is (5,7)
1659 {{{{0, 4, 6, 7}},
1660 {{4, 1, 5, 8}},
1661 {{6, 5, 2, 9}},
1662 {{7, 8, 9, 3}},
1663 {{4, 5, 6, 7}},
1664 {{4, 7, 8, 5}},
1665 {{6, 9, 7, 5}},
1666 {{5, 8, 9, 7}}}},
1667 // new line is (4,9)
1668 {{{{0, 4, 6, 7}},
1669 {{4, 1, 5, 8}},
1670 {{6, 5, 2, 9}},
1671 {{7, 8, 9, 3}},
1672 {{4, 5, 6, 9}},
1673 {{4, 7, 8, 9}},
1674 {{6, 9, 7, 4}},
1675 {{5, 8, 9, 4}}}}}};
1676 return cell_vertices_tet[refienement_choice];
1677 }
1678 default:
1680 }
1681
1682 return {};
1683}
1684
1685
1686
1687inline bool
1689{
1690 return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
1693}
1694
1695
1696
1697inline bool
1699{
1700 return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
1701 *this == ReferenceCells::Triangle ||
1703}
1704
1705
1706
1707inline unsigned int
1709{
1710 switch (this->kind)
1711 {
1713 return 0;
1715 return 1;
1718 return 2;
1723 return 3;
1724 default:
1726 }
1727
1729}
1730
1731
1732
1733template <int dim>
1736{
1737 Assert(dim == get_dimension(),
1738 ExcMessage("You can only call this function with arguments for "
1739 "which 'dim' equals the dimension of the space in which "
1740 "the current reference cell lives. You have dim=" +
1741 std::to_string(dim) + " but the reference cell is " +
1742 std::to_string(get_dimension()) + " dimensional."));
1743
1744 return Quadrature<dim>(std::vector<Point<dim>>({barycenter<dim>()}),
1745 std::vector<double>({volume()}));
1746}
1747
1748
1749
1750inline unsigned int
1752{
1753 switch (this->kind)
1754 {
1756 return 1;
1758 return 2;
1760 return 3;
1762 return 4;
1764 return 4;
1766 return 5;
1768 return 6;
1770 return 8;
1771 default:
1773 }
1774
1776}
1777
1778
1779
1780inline unsigned int
1782{
1783 switch (this->kind)
1784 {
1786 return 0;
1788 return 1;
1790 return 3;
1792 return 4;
1794 return 6;
1796 return 8;
1798 return 9;
1800 return 12;
1801 default:
1803 }
1804
1806}
1807
1808
1809
1810template <int dim>
1812ReferenceCell::vertex(const unsigned int v) const
1813{
1814 Assert(dim == get_dimension(),
1815 ExcMessage("You can only call this function with arguments for "
1816 "which 'dim' equals the dimension of the space in which "
1817 "the current reference cell lives. You have dim=" +
1818 std::to_string(dim) + " but the reference cell is " +
1819 std::to_string(get_dimension()) + " dimensional."));
1821
1822 switch (dim)
1823 {
1824 case 0:
1825 {
1826 if (*this == ReferenceCells::Vertex)
1827 return Point<dim>();
1828 break;
1829 }
1830 case 1:
1831 {
1832 static const Point<dim> vertices[2] = {
1833 Point<dim>(), // the origin
1834 Point<dim>::unit_vector(0) // unit point along x-axis
1835 };
1836 if (*this == ReferenceCells::Line)
1837 return vertices[v];
1838 break;
1839 }
1840 case 2:
1841 {
1842 switch (this->kind)
1843 {
1845 {
1846 static const Point<dim> vertices[3] = {
1847 Point<dim>(), // the origin
1848 Point<dim>::unit_vector(0), // unit point along x-axis
1849 Point<dim>::unit_vector(1) // unit point along y-axis
1850 };
1851 return vertices[v];
1852 }
1854 {
1855 static const Point<dim> vertices[4] = {
1856 // First the two points on the x-axis
1857 Point<dim>(),
1859 // Then these two points shifted in the y-direction
1862 return vertices[v];
1863 }
1864 }
1865 break;
1866 }
1867 case 3:
1868 {
1869 switch (this->kind)
1870 {
1872 {
1873 static const Point<dim> vertices[4] = {
1874 Point<dim>(), // the origin
1875 Point<dim>::unit_vector(0), // unit point along x-axis
1876 Point<dim>::unit_vector(1), // unit point along y-axis
1877 Point<dim>::unit_vector(2) // unit point along z-axis
1878 };
1879 return vertices[v];
1880 }
1882 {
1883 static const Point<dim> vertices[5] = {
1884 Point<dim>{-1.0, -1.0, 0.0},
1885 Point<dim>{+1.0, -1.0, 0.0},
1886 Point<dim>{-1.0, +1.0, 0.0},
1887 Point<dim>{+1.0, +1.0, 0.0},
1888 Point<dim>{+0.0, +0.0, 1.0}};
1889 return vertices[v];
1890 }
1892 {
1893 static const Point<dim> vertices[6] = {
1894 // First the three points on the triangular base of the
1895 // wedge:
1896 Point<dim>(),
1899 // And now everything shifted in the z-direction again
1903 return vertices[v];
1904 }
1906 {
1907 static const Point<dim> vertices[8] = {
1908 // First the two points on the x-axis
1909 Point<dim>(),
1911 // Then these two points shifted in the y-direction
1914 // And now all four points shifted in the z-direction
1921 return vertices[v];
1922 }
1923 }
1924 break;
1925 }
1926 default:
1928 }
1929
1931 return Point<dim>();
1932}
1933
1934
1935inline unsigned int
1937{
1938 switch (this->kind)
1939 {
1941 return 0;
1943 return 2;
1945 return 3;
1947 return 4;
1949 return 4;
1951 return 5;
1953 return 5;
1955 return 6;
1956 default:
1958 }
1959
1961}
1962
1963
1964
1971
1972
1973
1974inline unsigned int
1976{
1977 switch (this->kind)
1978 {
1980 return 1;
1982 return 1;
1984 return 1;
1986 return 1;
1988 return 3;
1990 return 1;
1992 return 1;
1994 return 1;
1995 default:
1997 }
1998
2000}
2001
2002
2003
2006 const unsigned int ref_choice) const
2007{
2009 switch (this->kind)
2010 {
2012 {
2014 isotropic_ref_choices = {{IsotropicRefinementChoice::cut_tet_68,
2017 return isotropic_ref_choices[ref_choice];
2018 }
2019 default:
2020 {
2022 }
2023 }
2024
2026}
2027
2028
2029
2030inline unsigned int
2032{
2033 switch (this->kind)
2034 {
2036 return 0;
2038 return 2;
2040 return 4;
2042 return 4;
2044 return 8;
2046 // We haven't yet decided how to refine pyramids. Update this when we
2047 // have
2048 return 0;
2050 return 8;
2052 return 8;
2053 default:
2055 }
2056
2058}
2059
2060
2061
2062template <int dim>
2063unsigned int
2065{
2066 Assert(dim == get_dimension(),
2067 ExcMessage("You can only call this function with arguments for "
2068 "which 'dim' equals the dimension of the space in which "
2069 "the current reference cell lives. You have dim=" +
2070 std::to_string(dim) + " but the reference cell is " +
2071 std::to_string(get_dimension()) + " dimensional."));
2072
2073 // Use GeometryInfo here to keep it the single source of truth
2074 if (this->is_hyper_cube())
2075 return GeometryInfo<dim>::n_children(ref_case);
2076 else
2077 return this->n_isotropic_children();
2078}
2079
2080
2081
2088
2089
2090
2097
2098
2099
2106
2107
2108
2109inline ReferenceCell
2110ReferenceCell::face_reference_cell(const unsigned int face_no) const
2111{
2112 AssertIndexRange(face_no, n_faces());
2113
2114 switch (this->kind)
2115 {
2122 return ReferenceCells::Line;
2126 if (face_no == 0)
2128 else
2131 if (face_no > 1)
2133 else
2137 default:
2139 }
2140
2142}
2143
2144
2145
2146inline constexpr types::geometric_orientation
2151
2152
2153
2154inline constexpr types::geometric_orientation
2156{
2157 // For a reversed line 'orientation' is false and neither flip nor rotate are
2158 // defined.
2160}
2161
2162
2163
2164inline unsigned int
2166 const unsigned int subface) const
2167{
2168 return child_cell_on_face(face,
2169 subface,
2171}
2172
2173
2174
2175inline unsigned int
2177 const unsigned int face,
2178 const unsigned int subface,
2179 const types::geometric_orientation combined_orientation) const
2180{
2181 AssertIndexRange(face, n_faces());
2183
2184 switch (this->kind)
2185 {
2188 {
2190 break;
2191 }
2193 {
2194 static constexpr ndarray<unsigned int, 3, 2> subcells = {
2195 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
2196
2197 Assert(combined_orientation ==
2199 combined_orientation == numbers::reverse_line_orientation,
2201 return subcells[face][combined_orientation ==
2203 subface :
2204 1 - subface];
2205 }
2207 {
2208 const auto [face_orientation, face_rotation, face_flip] =
2209 internal::split_face_orientation(combined_orientation);
2210
2213 face,
2214 subface,
2215 face_orientation,
2216 face_flip,
2217 face_rotation);
2218 }
2222 {
2224 break;
2225 }
2227 {
2228 const auto [face_orientation, face_rotation, face_flip] =
2229 internal::split_face_orientation(combined_orientation);
2230
2233 face,
2234 subface,
2235 face_orientation,
2236 face_flip,
2237 face_rotation);
2238 }
2239 default:
2241 }
2242
2244}
2245
2246
2247
2248inline std::array<unsigned int, 2>
2250 const unsigned int vertex) const
2251{
2253 // Work around a GCC warning at higher optimization levels by making all of
2254 // these tables the same size
2255 constexpr unsigned int X = numbers::invalid_unsigned_int;
2256
2257 switch (this->kind)
2258 {
2262 break;
2264 {
2265 static constexpr ndarray<unsigned int, 6, 2> table = {
2266 {{{0, 0}}, {{0, 1}}, {{1, 1}}, {{X, X}}, {{X, X}}, {{X, X}}}};
2267
2268 return table[vertex];
2269 }
2271 {
2273 vertex);
2274 }
2276 {
2277 static constexpr ndarray<unsigned int, 6, 2> table = {
2278 {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 2}}, {{X, X}}, {{X, X}}}};
2279
2280 return table[vertex];
2281 }
2283 {
2284 static constexpr ndarray<unsigned int, 6, 2> table = {
2285 {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{0, 3}}, {{1, 2}}, {{X, X}}}};
2286
2287 return table[vertex];
2288 }
2290 {
2291 static constexpr ndarray<unsigned int, 6, 2> table = {
2292 {{{0, 1}}, {{0, 0}}, {{0, 2}}, {{1, 0}}, {{1, 1}}, {{1, 2}}}};
2293
2294 return table[vertex];
2295 }
2297 {
2299 vertex);
2300 }
2301 default:
2303 }
2304
2305 return {};
2306}
2307
2308
2309
2310inline std::array<unsigned int, 2>
2312 const unsigned int line) const
2313{
2314 AssertIndexRange(line, n_lines());
2315
2316 switch (this->kind)
2317 {
2322 {
2324 break;
2325 }
2327 {
2328 static const std::array<unsigned int, 2> table[6] = {
2329 {{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}};
2330
2331 return table[line];
2332 }
2334 {
2335 static const std::array<unsigned int, 2> table[8] = {{{0, 0}},
2336 {{0, 1}},
2337 {{0, 2}},
2338 {{0, 3}},
2339 {{1, 2}},
2340 {{2, 1}},
2341 {{1, 1}},
2342 {{2, 2}}};
2343
2344 return table[line];
2345 }
2347 {
2348 static const std::array<unsigned int, 2> table[9] = {{{0, 0}},
2349 {{0, 2}},
2350 {{0, 1}},
2351 {{1, 0}},
2352 {{1, 1}},
2353 {{1, 2}},
2354 {{2, 0}},
2355 {{2, 1}},
2356 {{3, 1}}};
2357
2358 return table[line];
2359 }
2361 {
2363 }
2364 default:
2366 }
2367
2368 return {};
2369}
2370
2371
2372inline unsigned int
2374 const unsigned int vertex) const
2375{
2377 AssertIndexRange(line, n_lines());
2378
2379 switch (this->kind)
2380 {
2383 return vertex;
2385 {
2386 static constexpr ndarray<unsigned int, 3, 2> table = {
2387 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
2388 return table[line][vertex];
2389 }
2391 {
2392 static constexpr ndarray<unsigned int, 4, 2> table = {
2393 {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
2394 return table[line][vertex];
2395 }
2397 {
2398 static constexpr ndarray<unsigned int, 6, 2> table = {
2399 {{{0, 1}}, {{1, 2}}, {{2, 0}}, {{0, 3}}, {{1, 3}}, {{2, 3}}}};
2400 return table[line][vertex];
2401 }
2403 {
2404 static constexpr ndarray<unsigned int, 8, 2> table = {{{{0, 2}},
2405 {{1, 3}},
2406 {{0, 1}},
2407 {{2, 3}},
2408 {{4, 0}},
2409 {{1, 4}},
2410 {{2, 4}},
2411 {{4, 3}}}};
2412 return table[line][vertex];
2413 }
2415 {
2416 static constexpr ndarray<unsigned int, 9, 2> table = {{{{1, 0}},
2417 {{2, 1}},
2418 {{0, 2}},
2419 {{3, 4}},
2420 {{4, 5}},
2421 {{5, 3}},
2422 {{0, 3}},
2423 {{1, 4}},
2424 {{2, 5}}}};
2425 return table[line][vertex];
2426 }
2428 {
2429 // first four lines comprise the bottom face, next four are the top,
2430 // and the last four are 'bottom to top'
2431 static constexpr ndarray<unsigned int, 12, 2> table = {{{{0, 2}},
2432 {{1, 3}},
2433 {{0, 1}},
2434 {{2, 3}},
2435 {{4, 6}},
2436 {{5, 7}},
2437 {{4, 5}},
2438 {{6, 7}},
2439 {{0, 4}},
2440 {{1, 5}},
2441 {{2, 6}},
2442 {{3, 7}}}};
2443 return table[line][vertex];
2444 }
2445
2446 default:
2448 }
2449
2451}
2452
2453
2454inline unsigned int
2456 const unsigned int face,
2457 const unsigned int line,
2458 const types::geometric_orientation combined_face_orientation) const
2459{
2460 AssertIndexRange(face, n_faces());
2462
2463 static constexpr unsigned int X = numbers::invalid_unsigned_int;
2464
2465 switch (this->kind)
2466 {
2468 {
2470 break;
2471 }
2473 {
2474 const auto [face_orientation, face_rotation, face_flip] =
2475 internal::split_face_orientation(combined_face_orientation);
2476
2478 face, line, face_orientation, face_flip, face_rotation);
2479 }
2481 {
2482 return face;
2483 }
2485 {
2486 const auto [face_orientation, face_rotation, face_flip] =
2487 internal::split_face_orientation(combined_face_orientation);
2488
2490 face, line, face_orientation, face_flip, face_rotation);
2491 }
2493 {
2494 static constexpr ndarray<unsigned int, 4, 3> table = {
2495 {{{0, 1, 2}}, {{0, 3, 4}}, {{2, 5, 3}}, {{1, 4, 5}}}};
2496
2497 return table[face][standard_to_real_face_line(
2498 line, face, combined_face_orientation)];
2499 }
2501 {
2502 static constexpr ndarray<unsigned int, 5, 4> table = {
2503 {{{0, 1, 2, 3}},
2504 {{0, 6, 4, X}},
2505 {{1, 5, 7, X}},
2506 {{2, 4, 5, X}},
2507 {{3, 7, 6, X}}}};
2508
2509 return table[face][standard_to_real_face_line(
2510 line, face, combined_face_orientation)];
2511 }
2513 {
2514 static constexpr ndarray<unsigned int, 5, 4> table = {
2515 {{{0, 2, 1, X}},
2516 {{3, 4, 5, X}},
2517 {{6, 7, 0, 3}},
2518 {{7, 8, 1, 4}},
2519 {{8, 6, 5, 2}}}};
2520
2521 return table[face][standard_to_real_face_line(
2522 line, face, combined_face_orientation)];
2523 }
2525 {
2526 const auto [face_orientation, face_rotation, face_flip] =
2527 internal::split_face_orientation(combined_face_orientation);
2528
2530 face, line, face_orientation, face_flip, face_rotation);
2531 }
2532 default:
2534 }
2535
2537}
2538
2539
2540
2541inline unsigned int
2543 const unsigned int face,
2544 const unsigned int vertex,
2545 const types::geometric_orientation combined_face_orientation) const
2546{
2547 AssertIndexRange(face, n_faces());
2549 AssertIndexRange(combined_face_orientation, n_face_orientations(face));
2550
2551 switch (this->kind)
2552 {
2554 {
2556 break;
2557 }
2559 {
2560 const auto [face_orientation, face_rotation, face_flip] =
2561 internal::split_face_orientation(combined_face_orientation);
2562
2564 face, vertex, face_orientation, face_flip, face_rotation);
2565 }
2567 {
2568 static constexpr ndarray<unsigned int, 3, 2> table = {
2569 {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
2570
2571 return table[face][combined_face_orientation ==
2573 vertex :
2574 (1 - vertex)];
2575 }
2577 {
2578 const auto [face_orientation, face_rotation, face_flip] =
2579 internal::split_face_orientation(combined_face_orientation);
2580
2582 face, vertex, face_orientation, face_flip, face_rotation);
2583 }
2585 {
2586 static constexpr ndarray<unsigned int, 4, 3> table = {
2587 {{{0, 1, 2}}, {{1, 0, 3}}, {{0, 2, 3}}, {{2, 1, 3}}}};
2588
2589 return table[face][standard_to_real_face_vertex(
2590 vertex, face, combined_face_orientation)];
2591 }
2593 {
2594 constexpr auto X = numbers::invalid_unsigned_int;
2595 static constexpr ndarray<unsigned int, 5, 4> table = {
2596 {{{0, 1, 2, 3}},
2597 {{0, 2, 4, X}},
2598 {{3, 1, 4, X}},
2599 {{1, 0, 4, X}},
2600 {{2, 3, 4, X}}}};
2601
2602 return table[face][standard_to_real_face_vertex(
2603 vertex, face, combined_face_orientation)];
2604 }
2606 {
2607 constexpr auto X = numbers::invalid_unsigned_int;
2608 static constexpr ndarray<unsigned int, 6, 4> table = {
2609 {{{1, 0, 2, X}},
2610 {{3, 4, 5, X}},
2611 {{0, 1, 3, 4}},
2612 {{1, 2, 4, 5}},
2613 {{2, 0, 5, 3}}}};
2614
2615 return table[face][standard_to_real_face_vertex(
2616 vertex, face, combined_face_orientation)];
2617 }
2619 {
2620 const auto [face_orientation, face_rotation, face_flip] =
2621 internal::split_face_orientation(combined_face_orientation);
2622
2624 face, vertex, face_orientation, face_flip, face_rotation);
2625 }
2626 default:
2628 }
2629
2631}
2632
2633
2634
2635template <int dim>
2638 const unsigned int vertex) const
2639{
2640 Assert(dim == get_dimension(),
2641 ExcMessage("You can only call this function with arguments for "
2642 "which 'dim' equals the dimension of the space in which "
2643 "the current reference cell lives. You have dim=" +
2644 std::to_string(dim) + " but the reference cell is " +
2645 std::to_string(get_dimension()) + " dimensional."));
2646 return this->template vertex<dim>(face_to_cell_vertices(
2648}
2649
2650
2651
2652inline unsigned int
2654 const unsigned int vertex,
2655 const unsigned int face,
2656 const types::geometric_orientation face_orientation) const
2657{
2658 AssertIndexRange(face, n_faces());
2660
2661 switch (this->kind)
2662 {
2665 break;
2668 ExcMessage(
2669 "In 1D, all faces must have the default orientation."));
2670 return vertex;
2673 return line_vertex_permutations[face_orientation][vertex];
2675 return triangle_vertex_permutations[face_orientation][vertex];
2677 // face 0 is a quadrilateral
2678 if (face == 0)
2679 return quadrilateral_vertex_permutations[face_orientation][vertex];
2680 else
2681 return triangle_vertex_permutations[face_orientation][vertex];
2683 // faces 0 and 1 are triangles
2684 if (face > 1)
2685 return quadrilateral_vertex_permutations[face_orientation][vertex];
2686 else
2687 return triangle_vertex_permutations[face_orientation][vertex];
2689 return quadrilateral_vertex_permutations[face_orientation][vertex];
2690 default:
2692 }
2693
2696}
2697
2698
2699
2700inline unsigned int
2702 const unsigned int line,
2703 const unsigned int face,
2704 const types::geometric_orientation combined_face_orientation) const
2705{
2706 AssertIndexRange(face, n_faces());
2708
2709 static constexpr ndarray<unsigned int, 6, 3> triangle_table = {{{{0, 1, 2}},
2710 {{2, 1, 0}},
2711 {{2, 0, 1}},
2712 {{1, 0, 2}},
2713 {{1, 2, 0}},
2714 {{0, 2, 1}}}};
2715
2716
2717 switch (this->kind)
2718 {
2724 break;
2726 return triangle_table[combined_face_orientation][line];
2728 if (face == 0) // The quadrilateral face
2729 {
2730 const auto [face_orientation, face_rotation, face_flip] =
2731 internal::split_face_orientation(combined_face_orientation);
2732
2734 face_orientation,
2735 face_flip,
2736 face_rotation);
2737 }
2738 else // One of the triangular faces
2739 {
2740 return triangle_table[combined_face_orientation][line];
2741 }
2743 if (face > 1) // One of the quadrilateral faces
2744 {
2745 const auto [face_orientation, face_rotation, face_flip] =
2746 internal::split_face_orientation(combined_face_orientation);
2747
2749 face_orientation,
2750 face_flip,
2751 face_rotation);
2752 }
2753 else // One of the triangular faces
2754 return triangle_table[combined_face_orientation][line];
2756 {
2757 static constexpr ndarray<unsigned int, 8, 4> table = {{
2758 {{0, 1, 2, 3}},
2759 {{2, 3, 0, 1}},
2760 {{3, 2, 0, 1}},
2761 {{0, 1, 3, 2}},
2762 {{1, 0, 3, 2}},
2763 {{3, 2, 1, 0}},
2764 {{2, 3, 1, 0}},
2765 {{1, 0, 2, 3}},
2766 }};
2767 return table[combined_face_orientation][line];
2768 }
2769 default:
2771 }
2772
2774}
2775
2776
2777
2778namespace ReferenceCells
2779{
2780 template <int dim>
2781 inline constexpr const ReferenceCell &
2783 {
2784 switch (dim)
2785 {
2786 case 0:
2788 case 1:
2789 return ReferenceCells::Line;
2790 case 2:
2792 case 3:
2794 default:
2796 }
2798 }
2799
2800
2801
2802 template <int dim>
2803 inline constexpr const ReferenceCell &
2805 {
2806 switch (dim)
2807 {
2808 case 0:
2810 case 1:
2811 return ReferenceCells::Line;
2812 case 2:
2814 case 3:
2816 default:
2818 }
2820 }
2821
2822
2823
2824 template <int structdim>
2825 inline constexpr unsigned int
2830
2831
2832
2833 template <int structdim>
2834 inline constexpr unsigned int
2839
2840
2841
2842 template <int structdim>
2843 inline constexpr unsigned int
2848} // namespace ReferenceCells
2849
2850
2851inline ReferenceCell
2852ReferenceCell::n_vertices_to_type(const int dim, const unsigned int n_vertices)
2853{
2854 AssertIndexRange(dim, 4);
2856
2857 const auto X = ReferenceCells::Invalid;
2858 static const ndarray<ReferenceCell, 4, 9> table = {
2859 {// dim 0
2860 {{X, ReferenceCells::Vertex, X, X, X, X, X, X, X}},
2861 // dim 1
2862 {{X, X, ReferenceCells::Line, X, X, X, X, X, X}},
2863 // dim 2
2864 {{X,
2865 X,
2866 X,
2869 X,
2870 X,
2871 X,
2872 X}},
2873 // dim 3
2874 {{X,
2875 X,
2876 X,
2877 X,
2881 X,
2884 ExcMessage("The combination of dim = " + std::to_string(dim) +
2885 " and n_vertices = " + std::to_string(n_vertices) +
2886 " does not correspond to a known reference cell type."));
2887 return table[dim][n_vertices];
2888}
2889
2890
2891
2892template <int dim>
2893inline double
2895 const unsigned int i) const
2896{
2897 Assert(dim == get_dimension(),
2898 ExcMessage("You can only call this function with arguments for "
2899 "which 'dim' equals the dimension of the space in which "
2900 "the current reference cell lives. You have dim=" +
2901 std::to_string(dim) + " but the reference cell is " +
2902 std::to_string(get_dimension()) + " dimensional."));
2904 switch (this->kind)
2905 {
2907 return 1.0;
2911 if constexpr (dim > 0)
2914 // see also BarycentricPolynomials<2>::compute_value
2916 {
2917 switch (i)
2918 {
2919 case 0:
2920 return 1.0 - xi[std::min(0, dim - 1)] -
2921 xi[std::min(1, dim - 1)];
2922 case 1:
2923 return xi[std::min(0, dim - 1)];
2924 case 2:
2925 return xi[std::min(1, dim - 1)];
2926 default:
2928 }
2929 }
2930 // see also BarycentricPolynomials<3>::compute_value
2932 {
2933 switch (i)
2934 {
2935 case 0:
2936 return 1.0 - xi[std::min(0, dim - 1)] -
2937 xi[std::min(1, dim - 1)] - xi[std::min(2, dim - 1)];
2938 case 1:
2939 return xi[std::min(0, dim - 1)];
2940 case 2:
2941 return xi[std::min(1, dim - 1)];
2942 case 3:
2943 return xi[std::min(2, dim - 1)];
2944 default:
2946 }
2947 }
2948 // see also ScalarLagrangePolynomialPyramid::compute_value()
2950 {
2951 const double Q14 = 0.25;
2952
2953 const double r = xi[std::min(0, dim - 1)];
2954 const double s = xi[std::min(1, dim - 1)];
2955 const double t = xi[std::min(2, dim - 1)];
2956
2957 const double ratio =
2958 (std::fabs(t - 1.0) > 1.0e-14 ? (r * s * t) / (1.0 - t) : 0.0);
2959
2960 if (i == 0)
2961 return Q14 * ((1.0 - r) * (1.0 - s) - t + ratio);
2962 if (i == 1)
2963 return Q14 * ((1.0 + r) * (1.0 - s) - t - ratio);
2964 if (i == 2)
2965 return Q14 * ((1.0 - r) * (1.0 + s) - t - ratio);
2966 if (i == 3)
2967 return Q14 * ((1.0 + r) * (1.0 + s) - t + ratio);
2968 else
2969 return t;
2970 }
2971 // see also ScalarLagrangePolynomialWedge::compute_value()
2974 .d_linear_shape_function<2>(Point<2>(xi[std::min(0, dim - 1)],
2975 xi[std::min(1, dim - 1)]),
2976 i % 3) *
2978 .d_linear_shape_function<1>(Point<1>(xi[std::min(2, dim - 1)]),
2979 i / 3);
2980 default:
2982 }
2983
2985 return 0.0;
2986}
2987
2988
2989
2990template <int dim>
2991inline Tensor<1, dim>
2993 const unsigned int i) const
2994{
2995 Assert(dim == get_dimension(),
2996 ExcMessage("You can only call this function with arguments for "
2997 "which 'dim' equals the dimension of the space in which "
2998 "the current reference cell lives. You have dim=" +
2999 std::to_string(dim) + " but the reference cell is " +
3000 std::to_string(get_dimension()) + " dimensional."));
3001 switch (this->kind)
3002 {
3008 // see also BarycentricPolynomials<2>::compute_grad()
3010 switch (i)
3011 {
3012 case 0:
3013 return Point<dim>(-1.0, -1.0);
3014 case 1:
3015 return Point<dim>(+1.0, +0.0);
3016 case 2:
3017 return Point<dim>(+0.0, +1.0);
3018 default:
3020 }
3021 default:
3023 }
3024
3025 return Point<dim>(+0.0, +0.0, +0.0);
3026}
3027
3028
3029
3030inline double
3032{
3033 switch (this->kind)
3034 {
3036 return 1;
3038 return 1;
3040 return 1. / 2.;
3042 return 1;
3044 return 1. / 6.;
3046 return 4. / 3.;
3048 return 1. / 2.;
3050 return 1;
3051 default:
3053 }
3054
3055 return 0.0;
3056}
3057
3058
3059
3060inline double
3061ReferenceCell::face_measure(const unsigned int face_no) const
3062{
3063 AssertIndexRange(face_no, n_faces());
3064 constexpr double X = std::numeric_limits<double>::signaling_NaN();
3065
3066 static const std::array<double, 5> tri_faces{
3067 {1.0, std::sqrt(2.0), 1.0, X, X}};
3068 static const std::array<double, 5> tet_faces{
3069 {0.5, 0.5, 0.5, std::sqrt(3.0) / 2.0, X}};
3070 static const std::array<double, 5> pyramid_faces{
3071 {4, std::sqrt(2.0), std::sqrt(2.0), std::sqrt(2.0), std::sqrt(2.0)}};
3072 static const std::array<double, 5> wedge_faces{
3073 {0.5, 0.5, 1.0, std::sqrt(2.0), 1.0}};
3074
3075 switch (this->kind)
3076 {
3081 return 1.0;
3083 return tri_faces[face_no];
3085 return tet_faces[face_no];
3087 return pyramid_faces[face_no];
3089 return wedge_faces[face_no];
3090 }
3092 return X;
3093}
3094
3095
3096
3097template <int dim>
3098inline Point<dim>
3100{
3101 Assert(dim == get_dimension(),
3102 ExcMessage("You can only call this function with arguments for "
3103 "which 'dim' equals the dimension of the space in which "
3104 "the current reference cell lives. You have dim=" +
3105 std::to_string(dim) + " but the reference cell is " +
3106 std::to_string(get_dimension()) + " dimensional."));
3107
3108 switch (this->kind)
3109 {
3111 return Point<dim>();
3113 return Point<dim>(1. / 2.);
3115 return Point<dim>(1. / 3., 1. / 3.);
3117 return Point<dim>(1. / 2., 1. / 2.);
3119 return Point<dim>(1. / 4., 1. / 4., 1. / 4.);
3121 return Point<dim>(0, 0, 1. / 4.);
3123 return Point<dim>(1. / 3, 1. / 3, 1. / 2.);
3125 return Point<dim>(1. / 2., 1. / 2., 1. / 2.);
3126 default:
3128 }
3129
3130 return Point<dim>();
3131}
3132
3133
3134
3135template <int dim>
3136inline bool
3137ReferenceCell::contains_point(const Point<dim> &p, const double tolerance) const
3138{
3139 Assert(dim == get_dimension(),
3140 ExcMessage("You can only call this function with arguments for "
3141 "which 'dim' equals the dimension of the space in which "
3142 "the current reference cell lives. You have dim=" +
3143 std::to_string(dim) + " but the reference cell is " +
3144 std::to_string(get_dimension()) + " dimensional."));
3145
3146 // Introduce abbreviations to silence compiler warnings about invalid
3147 // array accesses (that can't happen, of course, but the compiler
3148 // doesn't know that).
3149 constexpr unsigned int x_coordinate = 0;
3150 constexpr unsigned int y_coordinate = (dim >= 2 ? 1 : 0);
3151 constexpr unsigned int z_coordinate = (dim >= 3 ? 2 : 0);
3152
3153 switch (this->kind)
3154 {
3156 {
3157 // Vertices are special cases in that they do not actually
3158 // have coordinates. Error out if this function is called
3159 // with a vertex:
3160 Assert(false,
3161 ExcMessage("Vertices are zero-dimensional objects and "
3162 "as a consequence have no coordinates. You "
3163 "cannot meaningfully ask whether a point is "
3164 "inside a vertex (within a certain tolerance) "
3165 "without coordinate values."));
3166 return false;
3167 }
3171 {
3172 for (unsigned int d = 0; d < dim; ++d)
3173 if ((p[d] < -tolerance) || (p[d] > 1 + tolerance))
3174 return false;
3175 return true;
3176 }
3179 {
3180 // First make sure that we are in the first quadrant or octant
3181 for (unsigned int d = 0; d < dim; ++d)
3182 if (p[d] < -tolerance)
3183 return false;
3184
3185 // Now we also need to make sure that we are below the diagonal line
3186 // or plane that delineates the simplex. This diagonal is given by
3187 // sum(p[d])<=1, and a diagonal a distance eps away is given by
3188 // sum(p[d])<=1+eps*sqrt(d). (For example, the point at (1,1) is a
3189 // distance of 1/sqrt(2) away from the diagonal. That is, its
3190 // sum satisfies
3191 // sum(p[d]) = 2 <= 1 + (1/sqrt(2)) * sqrt(2)
3192 // in other words, it satisfies the predicate with eps=1/sqrt(2).)
3193 double sum = 0;
3194 for (unsigned int d = 0; d < dim; ++d)
3195 sum += p[d];
3196 return (sum <= 1 + tolerance * std::sqrt(1. * dim));
3197 }
3199 {
3200 // A pyramid only lives in the upper half-space:
3201 if (p[z_coordinate] < -tolerance)
3202 return false;
3203
3204 // It also only lives in the space below z=1:
3205 if (p[z_coordinate] > 1 + tolerance)
3206 return false;
3207
3208 // Within what's left of the space, a pyramid is a cone that tapers
3209 // towards the top. First compute the distance of the point to the
3210 // axis in the max norm (this is the right norm because the vertices
3211 // of the pyramid are at points +/-1, +/-1):
3212 const double distance_from_axis =
3213 std::max(std::fabs(p[x_coordinate]), std::fabs(p[y_coordinate]));
3214
3215 // We are inside the pyramid if the distance from the axis is less
3216 // than (1-z)
3217 return (distance_from_axis <= 1 + tolerance - p[z_coordinate]);
3218 }
3220 {
3221 // The wedge we use is a triangle extruded into the third
3222 // dimension by one unit. So we can use the same logic as for
3223 // triangles above (i.e., for the simplex above, using dim==2)
3224 // and then check the third dimension separately.
3225
3226 if ((p[x_coordinate] < -tolerance) || (p[y_coordinate] < -tolerance))
3227 return false;
3228
3229 const double sum = p[x_coordinate] + p[y_coordinate];
3230 if (sum > 1 + tolerance * std::sqrt(2.0))
3231 return false;
3232
3233 if (p[z_coordinate] < -tolerance)
3234 return false;
3235 if (p[z_coordinate] > 1 + tolerance)
3236 return false;
3237
3238 return true;
3239 }
3240 default:
3242 }
3243
3244 return false;
3245}
3246
3247
3248
3249template <int dim>
3250inline Tensor<1, dim>
3252 const unsigned int i) const
3253{
3254 return face_tangent_vector<dim>(face_no, i);
3255}
3256
3257
3258
3259template <int dim>
3260inline Tensor<1, dim>
3261ReferenceCell::face_tangent_vector(const unsigned int face_no,
3262 const unsigned int i) const
3263{
3264 Assert(dim == get_dimension(),
3265 ExcMessage("You can only call this function with arguments for "
3266 "which 'dim' equals the dimension of the space in which "
3267 "the current reference cell lives. You have dim=" +
3268 std::to_string(dim) + " but the reference cell is " +
3269 std::to_string(get_dimension()) + " dimensional."));
3270 AssertIndexRange(i, dim - 1);
3271
3272 switch (this->kind)
3273 {
3281 {
3282 AssertIndexRange(face_no, 3);
3283 static const std::array<Tensor<1, dim>, 3> table = {
3284 {Point<dim>(1, 0),
3285 Point<dim>(-std::sqrt(0.5), +std::sqrt(0.5)),
3286 Point<dim>(0, -1)}};
3287
3288 return table[face_no];
3289 }
3291 {
3292 AssertIndexRange(face_no, 4);
3293 static const ndarray<Tensor<1, dim>, 4, 2> table = {
3294 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
3295 {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
3296 {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
3297 {{Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
3298 +std::pow(1.0 / 3.0, 1.0 / 4.0),
3299 0),
3300 Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
3301 0,
3302 +std::pow(1.0 / 3.0, 1.0 / 4.0))}}}};
3303
3304 return table[face_no][i];
3305 }
3307 {
3308 AssertIndexRange(face_no, 5);
3309 static const ndarray<Tensor<1, dim>, 5, 2> table = {
3310 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
3311 {{Point<dim>(+1.0 / std::sqrt(2.0), 0, +1.0 / std::sqrt(2.0)),
3312 Point<dim>(0, 1, 0)}},
3313 {{Point<dim>(+1.0 / std::sqrt(2.0), 0, -1.0 / std::sqrt(2.0)),
3314 Point<dim>(0, 1, 0)}},
3315 {{Point<dim>(1, 0, 0),
3316 Point<dim>(0, +1.0 / std::sqrt(2.0), +1.0 / std::sqrt(2.0))}},
3317 {{Point<dim>(1, 0, 0),
3318 Point<dim>(0, +1.0 / std::sqrt(2.0), -1.0 / std::sqrt(2.0))}}}};
3319
3320 return table[face_no][i];
3321 }
3323 {
3324 AssertIndexRange(face_no, 5);
3325 static const ndarray<Tensor<1, dim>, 5, 2> table = {
3326 {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
3327 {{Point<dim>(1, 0, 0), Point<dim>(0, 1, 0)}},
3328 {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
3329 {{Point<dim>(-1 / std::sqrt(2.0), +1 / std::sqrt(2.0), 0),
3330 Point<dim>(0, 0, 1)}},
3331 {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}}}};
3332
3333 return table[face_no][i];
3334 }
3335 default:
3337 }
3338
3339
3340 return {};
3341}
3342
3343
3344
3345template <int dim>
3346inline Tensor<1, dim>
3347ReferenceCell::unit_normal_vectors(const unsigned int face_no) const
3348{
3349 return face_normal_vector<dim>(face_no);
3350}
3351
3352
3353
3354template <int dim>
3355inline Tensor<1, dim>
3356ReferenceCell::face_normal_vector(const unsigned int face_no) const
3357{
3358 Assert(dim == get_dimension(),
3359 ExcMessage("You can only call this function with arguments for "
3360 "which 'dim' equals the dimension of the space in which "
3361 "the current reference cell lives. You have dim=" +
3362 std::to_string(dim) + " but the reference cell is " +
3363 std::to_string(get_dimension()) + " dimensional."));
3364
3365 if (is_hyper_cube())
3366 {
3369 }
3370 else if (dim == 2)
3371 {
3373
3374 // Return the rotated vector
3375 return cross_product_2d(face_tangent_vector<dim>(face_no, 0));
3376 }
3377 else if (dim == 3)
3378 {
3379 return cross_product_3d(face_tangent_vector<dim>(face_no, 0),
3380 face_tangent_vector<dim>(face_no, 1));
3381 }
3382
3384
3385 return {};
3386}
3387
3388
3389
3390inline unsigned int
3391ReferenceCell::n_face_orientations(const unsigned int face_no) const
3392{
3393 AssertIndexRange(face_no, n_faces());
3394 if (get_dimension() == 1)
3395 return 1;
3396 if (get_dimension() == 2)
3397 return 2;
3399 return 8;
3400 else if (face_reference_cell(face_no) == ReferenceCells::Triangle)
3401 return 6;
3402
3405}
3406
3407
3408
3411 const unsigned int line,
3412 const unsigned int face,
3413 const types::geometric_orientation combined_face_orientation,
3414 const types::geometric_orientation line_orientation) const
3415{
3417 face,
3418 combined_face_orientation,
3419 line_orientation);
3420}
3421
3422
3425 const unsigned int face_line_no,
3426 const unsigned int face_no,
3427 const types::geometric_orientation combined_face_orientation,
3428 const types::geometric_orientation line_orientation) const
3429{
3430 constexpr auto D = numbers::default_geometric_orientation;
3431 constexpr auto R = numbers::reverse_line_orientation;
3432 if (*this == ReferenceCells::Hexahedron)
3433 {
3434 static constexpr ::ndarray<types::geometric_orientation, 2, 8>
3435 table{{{{D, D, D, R, R, R, R, D}}, {{D, D, R, D, R, R, D, R}}}};
3436 // We use face_line_no / 2 here since lines i and i + 1 are parallel and,
3437 // on a given face, have the same relative orientations.
3438 const bool match =
3439 line_orientation == table[face_line_no / 2][combined_face_orientation];
3440
3443 }
3444 else if (*this == ReferenceCells::Tetrahedron)
3445 {
3446 static constexpr unsigned int X = numbers::invalid_unsigned_int;
3447 static constexpr ::ndarray<unsigned int, 4, 3> combined_lines{
3448 {{{0, 0, 0}}, {{X, 0, 1}}, {{X, 0, X}}, {{X, X, X}}}};
3449
3450 const auto combined_line = combined_lines[face_no][face_line_no];
3451
3452 Assert(combined_line != X,
3453 ExcMessage(
3454 "This function can only be called for following face-line "
3455 "combinations: (0,0), (0,1), (0,2), (1,1), (1,2), (2,1),"));
3456
3457 static constexpr ::ndarray<types::geometric_orientation, 2, 6>
3458 table{{{{D, R, D, R, D, R}}, {{R, D, R, D, R, D}}}};
3459
3460 const bool match =
3461 line_orientation == table[combined_line][combined_face_orientation];
3462
3465 }
3466 else
3467 // TODO: This might actually be wrong for some of the other
3468 // kinds of objects. We should check this
3470}
3471
3472
3473
3474namespace internal
3475{
3476 template <typename T>
3478 {
3479 public:
3493
3497 virtual ~NoPermutation() noexcept override = default;
3498
3502 virtual void
3503 print_info(std::ostream &out) const override
3504 {
3505 out << '[';
3506
3507 const unsigned int n_vertices = entity_type.n_vertices();
3508
3509 for (unsigned int i = 0; i < n_vertices; ++i)
3510 {
3511 out << vertices_0[i];
3512 if (i + 1 != n_vertices)
3513 out << ',';
3514 }
3515
3516 out << "] is not a valid permutation of [";
3517
3518 for (unsigned int i = 0; i < n_vertices; ++i)
3519 {
3520 out << vertices_1[i];
3521 if (i + 1 != n_vertices)
3522 out << ',';
3523 }
3524
3525 out << "]." << std::endl;
3526 }
3527
3532
3537
3542 };
3543
3555 << "The reference-cell type used on this cell (" << arg1.to_string()
3556 << ") does not match the reference-cell type of the finite element "
3557 << "associated with this cell (" << arg2.to_string() << "). "
3558 << "Did you accidentally use simplex elements on hypercube meshes "
3559 << "(or the other way around), or are you using a mixed mesh and "
3560 << "assigned a simplex element to a hypercube cell (or the other "
3561 << "way around) via the active_fe_index?");
3562} // namespace internal
3563
3564
3565
3566template <typename T, std::size_t N>
3568ReferenceCell::compute_orientation(const std::array<T, N> &vertices_0,
3569 const std::array<T, N> &vertices_1) const
3570{
3571 Assert(N >= n_vertices(),
3572 ExcMessage("The number of array elements must be equal to or "
3573 "greater than the number of vertices of the cell "
3574 "referenced by this object."));
3575
3576 // Call the non-deprecated function, taking care of calling it only with
3577 // those array elements that we actually care about (see the note
3578 // in the documentation about the arguments potentially being
3579 // larger arrays than necessary).
3581 make_array_view(vertices_0.begin(), vertices_0.begin() + n_vertices()),
3582 make_array_view(vertices_1.begin(), vertices_1.begin() + n_vertices()));
3583}
3584
3585
3586
3587template <typename T>
3590 const ArrayView<const T> &vertices_0,
3591 const ArrayView<const T> &vertices_1) const
3592{
3593 Assert(vertices_0.size() == n_vertices(),
3594 ExcMessage("The number of array elements must be equal to "
3595 "the number of vertices of the cell "
3596 "referenced by this object."));
3597 Assert(vertices_1.size() == n_vertices(),
3598 ExcMessage("The number of array elements must be equal to "
3599 "the number of vertices of the cell "
3600 "referenced by this object."));
3601
3602 auto compute_orientation = [&](const auto &table) {
3603 for (types::geometric_orientation o = 0; o < table.size(); ++o)
3604 {
3605 bool match = true;
3606 for (unsigned int j = 0; j < table[o].size(); ++j)
3607 match = (match && vertices_0[j] == vertices_1[table[o][j]]);
3608
3609 if (match)
3610 return o;
3611 }
3612
3613 // Do not use `this` in Assert because nvcc when using C++20 assumes that
3614 // `this` is an integer and we get the following error: invalid type
3615 // argument of unary '*' (have 'int')
3616 [[maybe_unused]] const auto &ref_cell = *this;
3617 Assert(false,
3618 (internal::NoPermutation<T>(ref_cell, vertices_0, vertices_1)));
3619 return std::numeric_limits<types::geometric_orientation>::max();
3620 };
3621
3622 switch (this->kind)
3623 {
3632 default:
3634 }
3635
3636 Assert(false, (internal::NoPermutation<T>(*this, vertices_0, vertices_1)));
3637 return std::numeric_limits<types::geometric_orientation>::max();
3638}
3639
3640
3641
3642template <typename T, std::size_t N>
3643inline std::array<T, N>
3645 const std::array<T, N> &vertices,
3646 const unsigned int orientation) const
3647{
3648 Assert(N >= n_vertices(),
3649 ExcMessage("The number of array elements must be equal to or "
3650 "greater than the number of vertices of the cell "
3651 "referenced by this object."));
3652
3653 // Call the non-deprecated function, taking care of calling it only with
3654 // those array elements that we actually care about (see the note
3655 // in the documentation about the arguments potentially being
3656 // larger arrays than necessary).
3657 const auto permutation = permute_by_combined_orientation(
3658 make_array_view(vertices.begin(), vertices.begin() + n_vertices()),
3659 orientation);
3660
3661 std::array<T, N> temp;
3662 std::copy(permutation.begin(), permutation.end(), temp.begin());
3663
3664 return temp;
3665}
3666
3667
3668
3669template <typename T>
3670boost::container::small_vector<T, 8>
3672 const ArrayView<const T> &vertices,
3673 const types::geometric_orientation orientation) const
3674{
3675 AssertDimension(vertices.size(), n_vertices());
3676 boost::container::small_vector<T, 8> result(n_vertices());
3677
3678 auto permute = [&](const auto &table) {
3679 AssertIndexRange(orientation, table.size());
3680 for (unsigned int j = 0; j < table[orientation].size(); ++j)
3681 result[j] = vertices[table[orientation][j]];
3682 };
3683
3684 switch (this->kind)
3685 {
3688 break;
3690 permute(line_vertex_permutations);
3691 break;
3694 break;
3697 break;
3698 default:
3700 }
3701
3702 return result;
3703}
3704
3705
3706
3709 const types::geometric_orientation orientation) const
3710{
3711 switch (this->kind)
3712 {
3714 // Things are always default-oriented in 1D
3715 return orientation;
3716
3718 // the 1d orientations are the identity and a flip: i.e., the identity
3719 // and an involutory mapping
3720 return orientation;
3721
3723 {
3724 AssertIndexRange(orientation, 6);
3725 constexpr std::array<types::geometric_orientation, 6> inverses{
3726 {0, 1, 4, 3, 2, 5}};
3727 return inverses[orientation];
3728 }
3730 {
3731 AssertIndexRange(orientation, 8);
3732 constexpr std::array<types::geometric_orientation, 8> inverses{
3733 {0, 1, 6, 3, 4, 5, 2, 7}};
3734 return inverses[orientation];
3735 }
3736 default:
3738 }
3739
3740 return std::numeric_limits<types::geometric_orientation>::max();
3741}
3742
3743
3744
3745template <>
3746unsigned int
3747ReferenceCell::vtk_lexicographic_to_node_index<0>(
3748 const std::array<unsigned, 0> &node_indices,
3749 const std::array<unsigned, 0> &nodes_per_direction,
3750 const bool legacy_format) const;
3751
3752template <>
3753unsigned int
3754ReferenceCell::vtk_lexicographic_to_node_index<1>(
3755 const std::array<unsigned, 1> &node_indices,
3756 const std::array<unsigned, 1> &nodes_per_direction,
3757 const bool legacy_format) const;
3758
3759template <>
3760unsigned int
3761ReferenceCell::vtk_lexicographic_to_node_index<2>(
3762 const std::array<unsigned, 2> &node_indices,
3763 const std::array<unsigned, 2> &nodes_per_direction,
3764 const bool legacy_format) const;
3765
3766template <>
3767unsigned int
3768ReferenceCell::vtk_lexicographic_to_node_index<3>(
3769 const std::array<unsigned, 3> &node_indices,
3770 const std::array<unsigned, 3> &nodes_per_direction,
3771 const bool legacy_format) const;
3772
3774
3775#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
std::pair< unsigned int, RefinementCase< dim - 1 > > equivalent_refinement_case(const types::geometric_orientation combined_face_orientation, const internal::SubfaceCase< dim > subface_case, const unsigned int subface_no) 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)