deal.II version GIT relicensing-6809-ge913b9bb34 2026-09-25 17:20:01+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
geometry_info.h
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 1998 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13#ifndef dealii_geometry_info_h
14#define dealii_geometry_info_h
15
16
17#include <deal.II/base/config.h>
18
22#include <deal.II/base/point.h>
24#include <deal.II/base/tensor.h>
25#include <deal.II/base/types.h>
26
27#include <Kokkos_Macros.hpp>
28
29#include <algorithm>
30#include <array>
31#include <cstddef>
32#include <cstdint>
33#include <ostream>
34#include <string>
35
36
37
39
40#ifndef DOXYGEN
41namespace internal
42{
43 namespace GeometryInfoHelper
44 {
45 // A struct that holds the values for all the arrays we want to initialize
46 // in GeometryInfo
47 template <int dim>
48 struct Initializers;
49
50 template <>
51 struct Initializers<1>
52 {
53 static constexpr std::array<unsigned int, 2>
54 ucd_to_deal()
55 {
56 return {{0, 1}};
57 }
58
59 static constexpr std::array<unsigned int, 2>
60 unit_normal_direction()
61 {
62 return {{0, 0}};
63 }
64
65 static constexpr std::array<int, 2>
66 unit_normal_orientation()
67 {
68 return {{-1, 1}};
69 }
70
71 static constexpr std::array<Tensor<1, 1>, 2>
72 unit_normal_vector()
73 {
74 return {{Tensor<1, 1>{{-1}}, Tensor<1, 1>{{1}}}};
75 }
76
77 static constexpr ::ndarray<Tensor<1, 1>, 2, 0>
78 unit_tangential_vectors()
79 {
80 return {{{{}}, {{}}}};
81 }
83 static constexpr std::array<unsigned int, 2>
84 opposite_face()
85 {
86 return {{1, 0}};
87 }
88
89 static constexpr std::array<unsigned int, 2>
90 dx_to_deal()
91 {
92 return {{0, 1}};
93 }
94
95 static constexpr ::ndarray<unsigned int, 2, 1>
96 vertex_to_face()
97 {
98 return {{{{0}}, {{1}}}};
99 }
100 };
101
102 template <>
103 struct Initializers<2>
104 {
105 static constexpr std::array<unsigned int, 4>
106 ucd_to_deal()
107 {
108 return {{0, 1, 3, 2}};
109 }
110
111 static constexpr std::array<unsigned int, 4>
112 unit_normal_direction()
113 {
114 return {{0, 0, 1, 1}};
115 }
116
117 static constexpr std::array<int, 4>
118 unit_normal_orientation()
119 {
120 return {{-1, 1, -1, 1}};
121 }
122
123 static constexpr std::array<Tensor<1, 2>, 4>
124 unit_normal_vector()
125 {
126 return {{Tensor<1, 2>{{-1., 0.}},
127 Tensor<1, 2>{{1., 0.}},
128 Tensor<1, 2>{{0., -1.}},
129 Tensor<1, 2>{{0., 1.}}}};
130 }
131
132 static constexpr ::ndarray<Tensor<1, 2>, 4, 1>
133 unit_tangential_vectors()
134 {
135 return {{{{Tensor<1, 2>{{0, -1}}}},
136 {{Tensor<1, 2>{{0, 1}}}},
137 {{Tensor<1, 2>{{1, 0}}}},
138 {{Tensor<1, 2>{{-1, 0}}}}}};
139 }
140
141 static constexpr std::array<unsigned int, 4>
142 opposite_face()
143 {
144 return {{1, 0, 3, 2}};
145 }
146
147 static constexpr std::array<unsigned int, 4>
148 dx_to_deal()
149 {
150 return {{0, 2, 1, 3}};
151 }
152
153 static constexpr ::ndarray<unsigned int, 4, 2>
154 vertex_to_face()
155 {
156 return {{{{0, 2}}, {{1, 2}}, {{0, 3}}, {{1, 3}}}};
157 }
158 };
159
160 template <>
161 struct Initializers<3>
162 {
163 static constexpr std::array<unsigned int, 8>
164 ucd_to_deal()
165 {
166 return {{0, 4, 5, 1, 2, 6, 7, 3}};
167 }
168
169 static constexpr std::array<unsigned int, 6>
170 unit_normal_direction()
171 {
172 return {{0, 0, 1, 1, 2, 2}};
173 }
174
175 static constexpr std::array<int, 6>
176 unit_normal_orientation()
177 {
178 return {{-1, 1, -1, 1, -1, 1}};
179 }
180
181 static constexpr std::array<Tensor<1, 3>, 6>
182 unit_normal_vector()
183 {
184 return {{Tensor<1, 3>{{-1, 0, 0}},
185 Tensor<1, 3>{{1, 0, 0}},
186 Tensor<1, 3>{{0, -1, 0}},
187 Tensor<1, 3>{{0, 1, 0}},
188 Tensor<1, 3>{{0, 0, -1}},
189 Tensor<1, 3>{{0, 0, 1}}}};
190 }
191
192 static constexpr ::ndarray<Tensor<1, 3>, 6, 2>
193 unit_tangential_vectors()
194 {
195 return {{{{Tensor<1, 3>{{0, -1, 0}}, Tensor<1, 3>{{0, 0, 1}}}},
196 {{Tensor<1, 3>{{0, 1, 0}}, Tensor<1, 3>{{0, 0, 1}}}},
197 {{Tensor<1, 3>{{0, 0, -1}}, Tensor<1, 3>{{1, 0, 0}}}},
198 {{Tensor<1, 3>{{0, 0, 1}}, Tensor<1, 3>{{1, 0, 0}}}},
199 {{Tensor<1, 3>{{-1, 0, 0}}, Tensor<1, 3>{{0, 1, 0}}}},
200 {{Tensor<1, 3>{{1, 0, 0}}, Tensor<1, 3>{{0, 1, 0}}}}}};
201 }
202
203 static constexpr std::array<unsigned int, 6>
204 opposite_face()
205 {
206 return {{1, 0, 3, 2, 5, 4}};
207 }
208
209 static constexpr std::array<unsigned int, 8>
210 dx_to_deal()
211 {
212 return {{0, 4, 2, 6, 1, 5, 3, 7}};
213 }
214
215 static constexpr ::ndarray<unsigned int, 8, 3>
216 vertex_to_face()
217 {
218 return {{{{0, 2, 4}},
219 {{1, 2, 4}},
220 {{0, 3, 4}},
221 {{1, 3, 4}},
222 {{0, 2, 5}},
223 {{1, 2, 5}},
224 {{0, 3, 5}},
225 {{1, 3, 5}}}};
226 }
227 };
228
229 template <>
230 struct Initializers<4>
231 {
232 static constexpr std::array<unsigned int, 16>
233 ucd_to_deal()
234 {
251 }
252
253 static constexpr std::array<unsigned int, 8>
254 unit_normal_direction()
255 {
256 return {{0, 0, 1, 1, 2, 2, 3, 3}};
257 }
258
259 static constexpr std::array<int, 8>
260 unit_normal_orientation()
261 {
262 return {{-1, 1, -1, 1, -1, 1, -1, 1}};
263 }
264
265 static constexpr std::array<Tensor<1, 4>, 8>
266 unit_normal_vector()
267 {
268 return {{Tensor<1, 4>{{-1, 0, 0, 0}},
269 Tensor<1, 4>{{1, 0, 0, 0}},
270 Tensor<1, 4>{{0, -1, 0, 0}},
271 Tensor<1, 4>{{0, 1, 0, 0}},
272 Tensor<1, 4>{{0, 0, -1, 0}},
273 Tensor<1, 4>{{0, 0, 1, 0}},
274 Tensor<1, 4>{{0, 0, 0, -1}},
275 Tensor<1, 4>{{0, 0, 0, 1}}}};
276 }
277
278 static constexpr ::ndarray<Tensor<1, 4>, 8, 3>
279 unit_tangential_vectors()
280 {
281 return {{{{Tensor<1, 4>{{0, -1, 0, 0}},
282 Tensor<1, 4>{{0, 0, 1, 0}},
283 Tensor<1, 4>{{0, 0, 0, 1}}}},
284 {{Tensor<1, 4>{{0, 1, 0, 0}},
285 Tensor<1, 4>{{0, 0, 1, 0}},
286 Tensor<1, 4>{{0, 0, 0, 1}}}},
287 {{Tensor<1, 4>{{0, 0, -1, 0}},
288 Tensor<1, 4>{{0, 0, 0, 1}},
289 Tensor<1, 4>{{1, 0, 0, 0}}}},
290 {{Tensor<1, 4>{{0, 0, 1, 0}},
291 Tensor<1, 4>{{0, 0, 0, 1}},
292 Tensor<1, 4>{{1, 0, 0, 0}}}},
293 {{Tensor<1, 4>{{0, 0, 0, -1}},
294 Tensor<1, 4>{{1, 0, 0, 0}},
295 Tensor<1, 4>{{0, 1, 0, 0}}}},
296 {{Tensor<1, 4>{{0, 0, 0, 1}},
297 Tensor<1, 4>{{1, 0, 0, 0}},
298 Tensor<1, 4>{{0, 1, 0, 0}}}},
299 {{Tensor<1, 4>{{-1, 0, 0, 0}},
300 Tensor<1, 4>{{0, 1, 0, 0}},
301 Tensor<1, 4>{{0, 0, 1, 0}}}},
302 {{Tensor<1, 4>{{1, 0, 0, 0}},
303 Tensor<1, 4>{{0, 1, 0, 0}},
304 Tensor<1, 4>{{0, 0, 1, 0}}}}}};
305 }
306
307 static constexpr std::array<unsigned int, 8>
308 opposite_face()
309 {
310 return {{1, 0, 3, 2, 5, 4, 7, 6}};
311 }
312
313 static constexpr std::array<unsigned int, 16>
314 dx_to_deal()
315 {
332 }
333
334 static constexpr ::ndarray<unsigned int, 16, 4>
335 vertex_to_face()
336 {
401 }
402 };
403 } // namespace GeometryInfoHelper
404} // namespace internal
405#endif // DOXYGEN
406
407
424{
425public:
433 {
441 line = 1,
445 quad = 2,
449 hex = 3
450 };
451
457
463 GeometryPrimitive(const unsigned int object_dimension);
464
469 operator unsigned int() const;
470
471private:
476};
477
478
479
490template <int dim>
492{
527 enum Possibilities : std::uint8_t
528 {
533
546 };
547};
548
549
550
558template <>
560{
595 enum Possibilities : std::uint8_t
596 {
604 cut_x = 1,
609 };
610};
611
612
613
622template <>
624{
659 enum Possibilities : std::uint8_t
660 {
668 cut_x = 1,
672 cut_y = 2,
676 cut_xy = cut_x | cut_y,
677
681 isotropic_refinement = cut_xy
682 };
683};
684
685
686
695template <>
697{
732 enum Possibilities : std::uint8_t
733 {
741 cut_x = 1,
745 cut_y = 2,
749 cut_xy = cut_x | cut_y,
753 cut_z = 4,
757 cut_xz = cut_x | cut_z,
761 cut_yz = cut_y | cut_z,
765 cut_xyz = cut_x | cut_y | cut_z,
766
770 isotropic_refinement = cut_xyz
771 };
772};
773
774
775
784template <int dim>
786{
787public:
791 static constexpr unsigned int n_refinement_cases = (1 << dim);
792
796 constexpr RefinementCase();
797
802 constexpr RefinementCase(
803 const typename RefinementPossibilities<dim>::Possibilities refinement_case);
804
810 constexpr explicit RefinementCase(const std::uint8_t refinement_case);
811
824 constexpr operator std::uint8_t() const;
825
830 constexpr RefinementCase
831 operator|(const RefinementCase &r) const;
832
837 constexpr RefinementCase
838 operator&(const RefinementCase &r) const;
839
847 constexpr RefinementCase
848 operator~() const;
849
850
856 static constexpr RefinementCase
857 cut_axis(const unsigned int i);
858
867 static constexpr std::array<RefinementCase<dim>, n_refinement_cases>
869
873 static constexpr std::size_t
875
881 template <class Archive>
882 void
883 serialize(Archive &ar, const unsigned int version);
884
890 int,
891 << "The refinement flags given (" << arg1
892 << ") contain set bits that do not "
893 << "make sense for the space dimension of the object to which they are applied.");
894
895private:
900 std::uint8_t value : (dim > 0 ? dim : 1);
901};
902
903
904
905namespace internal
906{
924 template <int dim>
926 {
931 {
936
940 case_isotropic = static_cast<std::uint8_t>(-1)
941 };
942 };
943
944
951 template <>
953 {
971 };
972
973
974
982 template <>
984 {
1002 };
1003
1004
1005
1014 template <>
1016 {
1024 {
1032 case_x = 1,
1036 case_isotropic = case_x
1038 };
1039
1040
1041
1131 template <>
1133 {
1141 {
1143 case_x = 1,
1144 case_x1y = 2,
1145 case_x2y = 3,
1146 case_x1y2y = 4,
1147 case_y = 5,
1148 case_y1x = 6,
1149 case_y2x = 7,
1150 case_y1x2x = 8,
1151 case_xy = 9,
1152
1153 case_isotropic = case_xy
1155 };
1156
1157
1158
1163 template <int dim>
1165 {
1166 public:
1173 subface_possibility);
1174
1186 operator std::uint8_t() const;
1187
1191 static constexpr std::size_t
1193
1199 int,
1200 << "The subface case given (" << arg1 << ") does not make sense "
1201 << "for the space dimension of the object to which they are applied.");
1202
1203 private:
1208 std::uint8_t value : (dim == 3 ? 4 : 1);
1209 };
1210
1211} // namespace internal
1212
1213
1214
1215template <int dim>
1216struct GeometryInfo;
1217
1218
1219
1239template <>
1241{
1249 static constexpr unsigned int max_children_per_cell = 1;
1250
1254 static constexpr unsigned int faces_per_cell = 0;
1255
1272 static std::array<unsigned int, 0>
1274
1282 static constexpr unsigned int max_children_per_face = 0;
1283
1289 static unsigned int
1290 n_children(const RefinementCase<0> &refinement_case);
1291
1295 static constexpr unsigned int vertices_per_cell = 1;
1296
1315 static std::array<unsigned int, vertices_per_cell>
1317
1341 static unsigned int
1342 face_to_cell_vertices(const unsigned int face,
1343 const unsigned int vertex,
1344 const bool face_orientation = true,
1345 const bool face_flip = false,
1346 const bool face_rotation = false);
1347
1362 static unsigned int
1363 face_to_cell_lines(const unsigned int face,
1364 const unsigned int line,
1365 const bool face_orientation = true,
1366 const bool face_flip = false,
1367 const bool face_rotation = false);
1368
1375 static constexpr unsigned int vertices_per_face = 0;
1376
1380 static constexpr unsigned int lines_per_face = 0;
1381
1385 static constexpr unsigned int quads_per_face = 0;
1386
1390 static constexpr unsigned int lines_per_cell = 0;
1391
1395 static constexpr unsigned int quads_per_cell = 0;
1396
1400 static constexpr unsigned int hexes_per_cell = 0;
1401
1419 static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1420
1434 static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1435};
1436
1437
1438
1969template <int dim>
1971{
1979 static constexpr unsigned int max_children_per_cell = 1 << dim;
1980
1984 static constexpr unsigned int faces_per_cell = 2 * dim;
1985
2004
2012 static constexpr unsigned int max_children_per_face =
2014
2018 static constexpr unsigned int vertices_per_cell = 1 << dim;
2019
2037
2041 static constexpr unsigned int vertices_per_face =
2043
2047 static constexpr unsigned int lines_per_face =
2049
2053 static constexpr unsigned int quads_per_face =
2055
2065 static constexpr unsigned int lines_per_cell =
2066 (2 * GeometryInfo<dim - 1>::lines_per_cell +
2068
2076 static constexpr unsigned int quads_per_cell =
2077 (2 * GeometryInfo<dim - 1>::quads_per_cell +
2079
2083 static constexpr unsigned int hexes_per_cell =
2084 (2 * GeometryInfo<dim - 1>::hexes_per_cell +
2086
2104 static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
2105 internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal();
2106
2120 static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
2121 internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal();
2122
2135 internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face();
2136
2141 static unsigned int
2142 n_children(const RefinementCase<dim> &refinement_case);
2143
2148 static unsigned int
2150
2160 static double
2162 const unsigned int subface_no);
2163
2169 static RefinementCase<dim - 1>
2170 face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2171 const unsigned int face_no,
2172 const bool face_orientation = true,
2173 const bool face_flip = false,
2174 const bool face_rotation = false);
2175
2181 static RefinementCase<dim>
2184 const unsigned int face_no,
2185 const bool face_orientation = true,
2186 const bool face_flip = false,
2187 const bool face_rotation = false);
2188
2193 static RefinementCase<1>
2194 line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2195 const unsigned int line_no);
2196
2201 static RefinementCase<dim>
2203
2250 static unsigned int
2252 const unsigned int face,
2253 const unsigned int subface,
2254 const bool face_orientation = true,
2255 const bool face_flip = false,
2256 const bool face_rotation = false,
2259
2273 static unsigned int
2274 line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
2275
2296 static unsigned int
2297 face_to_cell_vertices(const unsigned int face,
2298 const unsigned int vertex,
2299 const bool face_orientation = true,
2300 const bool face_flip = false,
2301 const bool face_rotation = false);
2302
2314 static unsigned int
2315 face_to_cell_lines(const unsigned int face,
2316 const unsigned int line,
2317 const bool face_orientation = true,
2318 const bool face_flip = false,
2319 const bool face_rotation = false);
2320
2330 static unsigned int
2331 standard_to_real_face_vertex(const unsigned int vertex,
2332 const bool face_orientation = true,
2333 const bool face_flip = false,
2334 const bool face_rotation = false);
2335
2345 static unsigned int
2346 real_to_standard_face_vertex(const unsigned int vertex,
2347 const bool face_orientation = true,
2348 const bool face_flip = false,
2349 const bool face_rotation = false);
2350
2360 static unsigned int
2361 standard_to_real_face_line(const unsigned int line,
2362 const bool face_orientation = true,
2363 const bool face_flip = false,
2364 const bool face_rotation = false);
2365
2371 static unsigned int
2372 standard_to_real_line_vertex(const unsigned int vertex,
2373 const bool line_orientation = true);
2374
2382 static std::array<unsigned int, 2>
2384
2392 static std::array<unsigned int, 2>
2394
2402 static std::array<unsigned int, 2>
2404
2414 static unsigned int
2415 real_to_standard_face_line(const unsigned int line,
2416 const bool face_orientation = true,
2417 const bool face_flip = false,
2418 const bool face_rotation = false);
2419
2425 static Point<dim>
2426 unit_cell_vertex(const unsigned int vertex);
2427
2437 static unsigned int
2439
2447 static Point<dim>
2449 const unsigned int child_index,
2450 const RefinementCase<dim> refine_case =
2452
2458 static Point<dim>
2460 const unsigned int child_index,
2461 const RefinementCase<dim> refine_case =
2463
2468 static bool
2470
2482 static bool
2483 is_inside_unit_cell(const Point<dim> &p, const double eps);
2484
2489 template <typename Number = double>
2490 static Point<dim, Number>
2492
2498 static double
2500
2505 static double
2506 d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2507
2512 static Tensor<1, dim>
2513 d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2514
2566 template <int spacedim>
2567 static void
2569#ifndef DEAL_II_CXX14_CONSTEXPR_BUG
2570 (const Point<spacedim> (&vertices)[vertices_per_cell],
2571 Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2572#else
2573 (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2574#endif
2575
2585 static constexpr std::array<unsigned int, faces_per_cell>
2587 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction();
2588
2605 static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2606 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation();
2607
2618 static constexpr std::array<Tensor<1, dim>, faces_per_cell>
2620 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_vector();
2621
2635 static constexpr ndarray<Tensor<1, dim>, faces_per_cell, dim - 1>
2636 unit_tangential_vectors = internal::GeometryInfoHelper::Initializers<
2638
2644 static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2645 internal::GeometryInfoHelper::Initializers<dim>::opposite_face();
2646
2647
2652 double,
2653 << "The coordinates must satisfy 0 <= x_i <= 1, "
2654 << "but here we have x_i=" << arg1);
2655
2660 int,
2661 int,
2662 int,
2663 << "RefinementCase<dim> " << arg1 << ": face " << arg2
2664 << " has no subface " << arg3);
2665};
2666
2667
2668
2669#ifndef DOXYGEN
2670
2671
2672/* -------------- declaration of explicit specializations ------------- */
2673
2674
2675template <>
2678 const unsigned int i);
2679template <>
2682 const unsigned int i);
2683template <>
2686 const unsigned int i);
2687
2688
2689
2690/* -------------- inline functions ------------- */
2691
2692
2693inline GeometryPrimitive::GeometryPrimitive(const Object object)
2694 : object(object)
2695{}
2696
2697
2698
2699inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2700 : object(static_cast<Object>(object_dimension))
2701{}
2702
2703
2704inline GeometryPrimitive::operator unsigned int() const
2705{
2706 return static_cast<unsigned int>(object);
2707}
2708
2709
2710
2711namespace internal
2712{
2713 template <int dim>
2715 const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2716 : value(subface_possibility)
2717 {}
2718
2719
2720 template <int dim>
2721 inline SubfaceCase<dim>::operator std::uint8_t() const
2722 {
2723 return value;
2724 }
2725
2726
2727} // namespace internal
2728
2729
2730template <int dim>
2731inline constexpr RefinementCase<dim>
2732RefinementCase<dim>::cut_axis(const unsigned int)
2733{
2735 return static_cast<std::uint8_t>(-1);
2736}
2737
2738
2739template <>
2740inline constexpr RefinementCase<1>
2741RefinementCase<1>::cut_axis(const unsigned int i)
2742{
2743 AssertIndexRange(i, 1);
2744
2746 return options[i];
2747}
2748
2749
2750
2751template <>
2752inline constexpr RefinementCase<2>
2753RefinementCase<2>::cut_axis(const unsigned int i)
2754{
2755 AssertIndexRange(i, 2);
2756
2759 return options[i];
2760}
2761
2762
2763
2764template <>
2765inline constexpr RefinementCase<3>
2766RefinementCase<3>::cut_axis(const unsigned int i)
2767{
2768 AssertIndexRange(i, 3);
2769
2773 return options[i];
2774}
2775
2776
2777
2778template <>
2779inline constexpr std::array<RefinementCase<1>, 2>
2781{
2784}
2785
2786
2787
2788template <>
2789inline constexpr std::array<RefinementCase<2>, 4>
2791{
2796}
2797
2798
2799
2800template <>
2801inline constexpr std::array<RefinementCase<3>, 8>
2803{
2812}
2813
2814
2815
2816template <int dim>
2817inline constexpr RefinementCase<dim>::RefinementCase()
2818 : value(RefinementPossibilities<dim>::no_refinement)
2819{}
2820
2821
2822
2823template <int dim>
2825 const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2826 : value(refinement_case)
2827{
2828 // check that only those bits of
2829 // the given argument are set that
2830 // make sense for a given space
2831 // dimension
2832 Assert((refinement_case &
2834 refinement_case,
2835 ExcInvalidRefinementCase(refinement_case));
2836}
2837
2838
2839
2840template <int dim>
2842 const std::uint8_t refinement_case)
2843 : value(refinement_case)
2844{
2845 // check that only those bits of
2846 // the given argument are set that
2847 // make sense for a given space
2848 // dimension
2849 Assert((refinement_case &
2851 refinement_case,
2852 ExcInvalidRefinementCase(refinement_case));
2853}
2854
2855
2856
2857template <int dim>
2858inline constexpr DEAL_II_HOST_DEVICE
2859 RefinementCase<dim>::operator std::uint8_t() const
2860{
2861 return value;
2862}
2863
2864
2865
2866template <int dim>
2867inline constexpr RefinementCase<dim>
2869{
2870 return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2871}
2872
2873
2874
2875template <int dim>
2876inline constexpr RefinementCase<dim>
2878{
2879 return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2880}
2881
2882
2883
2884template <int dim>
2885inline constexpr RefinementCase<dim>
2887{
2888 return RefinementCase<dim>(static_cast<std::uint8_t>(
2890}
2891
2892
2893
2894template <int dim>
2895inline constexpr std::size_t
2897{
2898 return sizeof(RefinementCase<dim>);
2899}
2900
2901
2902
2903template <int dim>
2904template <class Archive>
2905inline void
2906RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2907{
2908 // serialization can't deal with bitfields, so copy from/to a full sized
2909 // std::uint8_t
2910 std::uint8_t uchar_value = value;
2911 ar &uchar_value;
2912 value = uchar_value;
2913}
2914
2915
2916
2917template <>
2918inline Point<1>
2919GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2920{
2921 AssertIndexRange(vertex, vertices_per_cell);
2922
2923 return Point<1>(static_cast<double>(vertex));
2924}
2925
2926
2927
2928template <>
2929inline Point<2>
2930GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2931{
2932 AssertIndexRange(vertex, vertices_per_cell);
2933
2934 return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2935}
2936
2937
2938
2939template <>
2940inline Point<3>
2941GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2942{
2943 AssertIndexRange(vertex, vertices_per_cell);
2944
2945 return {static_cast<double>(vertex % 2),
2946 static_cast<double>(vertex / 2 % 2),
2947 static_cast<double>(vertex / 4)};
2948}
2949
2950
2951
2952inline std::array<unsigned int, 0>
2954{
2955 return {{}};
2956}
2957
2958
2959
2960inline std::array<unsigned int, 1>
2962{
2963 return {{0}};
2964}
2965
2966
2967
2968template <int dim>
2971{
2973 0U, faces_per_cell);
2974}
2975
2976
2977
2978template <int dim>
2981{
2983 0U, vertices_per_cell);
2984}
2985
2986
2987
2988template <int dim>
2989inline Point<dim>
2990GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2991{
2993
2994 return {};
2995}
2996
2997
2998
2999template <>
3000inline unsigned int
3002{
3003 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
3004
3005 return (p[0] <= 0.5 ? 0 : 1);
3006}
3007
3008
3009
3010template <>
3011inline unsigned int
3013{
3014 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
3015 Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
3016
3017 return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
3018}
3019
3020
3021
3022template <>
3023inline unsigned int
3025{
3026 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
3027 Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
3028 Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
3029
3030 return (p[0] <= 0.5 ?
3031 (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
3032 (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
3033}
3034
3035
3036template <int dim>
3037inline unsigned int
3039{
3041
3042 return 0;
3043}
3044
3045
3046
3047template <>
3048inline Point<1>
3050 const unsigned int child_index,
3051 const RefinementCase<1> refine_case)
3052
3053{
3054 AssertIndexRange(child_index, 2);
3056 (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3057
3058 return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
3059}
3060
3061
3062
3063template <>
3064inline Point<2>
3066 const unsigned int child_index,
3067 const RefinementCase<2> refine_case)
3068
3069{
3070 AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3071
3072 Point<2> point = p;
3073 switch (refine_case)
3074 {
3076 point[0] *= 2.0;
3077 if (child_index == 1)
3078 point[0] -= 1.0;
3079 break;
3081 point[1] *= 2.0;
3082 if (child_index == 1)
3083 point[1] -= 1.0;
3084 break;
3086 point *= 2.0;
3087 point -= unit_cell_vertex(child_index);
3088 break;
3089 default:
3091 }
3092
3093 return point;
3094}
3095
3096
3097
3098template <>
3099inline Point<3>
3101 const unsigned int child_index,
3102 const RefinementCase<3> refine_case)
3103
3104{
3105 AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3106
3107 Point<3> point = p;
3108 // there might be a cleverer way to do
3109 // this, but since this function is called
3110 // in very few places for initialization
3111 // purposes only, I don't care at the
3112 // moment
3113 switch (refine_case)
3114 {
3116 point[0] *= 2.0;
3117 if (child_index == 1)
3118 point[0] -= 1.0;
3119 break;
3121 point[1] *= 2.0;
3122 if (child_index == 1)
3123 point[1] -= 1.0;
3124 break;
3126 point[2] *= 2.0;
3127 if (child_index == 1)
3128 point[2] -= 1.0;
3129 break;
3131 point[0] *= 2.0;
3132 point[1] *= 2.0;
3133 if (child_index % 2 == 1)
3134 point[0] -= 1.0;
3135 if (child_index / 2 == 1)
3136 point[1] -= 1.0;
3137 break;
3139 // careful, this is slightly
3140 // different from xy and yz due to
3141 // different internal numbering of
3142 // children!
3143 point[0] *= 2.0;
3144 point[2] *= 2.0;
3145 if (child_index / 2 == 1)
3146 point[0] -= 1.0;
3147 if (child_index % 2 == 1)
3148 point[2] -= 1.0;
3149 break;
3151 point[1] *= 2.0;
3152 point[2] *= 2.0;
3153 if (child_index % 2 == 1)
3154 point[1] -= 1.0;
3155 if (child_index / 2 == 1)
3156 point[2] -= 1.0;
3157 break;
3159 point *= 2.0;
3160 point -= unit_cell_vertex(child_index);
3161 break;
3162 default:
3164 }
3165
3166 return point;
3167}
3168
3169
3170
3171template <int dim>
3172inline Point<dim>
3174 const Point<dim> & /*p*/,
3175 const unsigned int /*child_index*/,
3176 const RefinementCase<dim> /*refine_case*/)
3177
3178{
3180 return {};
3181}
3182
3183
3184
3185template <>
3186inline Point<1>
3188 const unsigned int child_index,
3189 const RefinementCase<1> refine_case)
3190
3191{
3192 AssertIndexRange(child_index, 2);
3194 (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3195
3196 return (p + unit_cell_vertex(child_index)) * 0.5;
3197}
3198
3199
3200
3201template <>
3202inline Point<3>
3204 const unsigned int child_index,
3205 const RefinementCase<3> refine_case)
3206
3207{
3208 AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3209
3210 Point<3> point = p;
3211 // there might be a cleverer way to do
3212 // this, but since this function is called
3213 // in very few places for initialization
3214 // purposes only, I don't care at the
3215 // moment
3216 switch (refine_case)
3217 {
3219 if (child_index == 1)
3220 point[0] += 1.0;
3221 point[0] *= 0.5;
3222 break;
3224 if (child_index == 1)
3225 point[1] += 1.0;
3226 point[1] *= 0.5;
3227 break;
3229 if (child_index == 1)
3230 point[2] += 1.0;
3231 point[2] *= 0.5;
3232 break;
3234 if (child_index % 2 == 1)
3235 point[0] += 1.0;
3236 if (child_index / 2 == 1)
3237 point[1] += 1.0;
3238 point[0] *= 0.5;
3239 point[1] *= 0.5;
3240 break;
3242 // careful, this is slightly
3243 // different from xy and yz due to
3244 // different internal numbering of
3245 // children!
3246 if (child_index / 2 == 1)
3247 point[0] += 1.0;
3248 if (child_index % 2 == 1)
3249 point[2] += 1.0;
3250 point[0] *= 0.5;
3251 point[2] *= 0.5;
3252 break;
3254 if (child_index % 2 == 1)
3255 point[1] += 1.0;
3256 if (child_index / 2 == 1)
3257 point[2] += 1.0;
3258 point[1] *= 0.5;
3259 point[2] *= 0.5;
3260 break;
3262 point += unit_cell_vertex(child_index);
3263 point *= 0.5;
3264 break;
3265 default:
3267 }
3268
3269 return point;
3270}
3271
3272
3273
3274template <>
3275inline Point<2>
3277 const unsigned int child_index,
3278 const RefinementCase<2> refine_case)
3279{
3280 AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3281
3282 Point<2> point = p;
3283 switch (refine_case)
3284 {
3286 if (child_index == 1)
3287 point[0] += 1.0;
3288 point[0] *= 0.5;
3289 break;
3291 if (child_index == 1)
3292 point[1] += 1.0;
3293 point[1] *= 0.5;
3294 break;
3296 point += unit_cell_vertex(child_index);
3297 point *= 0.5;
3298 break;
3299 default:
3301 }
3302
3303 return point;
3304}
3305
3306
3307
3308template <int dim>
3309inline Point<dim>
3311 const Point<dim> & /*p*/,
3312 const unsigned int /*child_index*/,
3313 const RefinementCase<dim> /*refine_case*/)
3314{
3316 return {};
3317}
3318
3319
3320
3321template <int dim>
3322inline bool
3324{
3326 return false;
3327}
3328
3329template <>
3330inline bool
3332{
3333 return (p[0] >= 0.) && (p[0] <= 1.);
3334}
3335
3336
3337
3338template <>
3339inline bool
3341{
3342 return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
3343}
3344
3345
3346
3347template <>
3348inline bool
3350{
3351 return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
3352 (p[2] >= 0.) && (p[2] <= 1.);
3353}
3354
3355
3356
3357template <int dim>
3358inline bool
3360{
3362 return false;
3363}
3364
3365template <>
3366inline bool
3367GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
3368{
3369 return (p[0] >= -eps) && (p[0] <= 1. + eps);
3370}
3371
3372
3373
3374template <>
3375inline bool
3376GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
3377{
3378 const double l = -eps, u = 1 + eps;
3379 return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
3380}
3381
3382
3383
3384template <>
3385inline bool
3386GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
3387{
3388 const double l = -eps, u = 1.0 + eps;
3389 return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
3390 (p[2] >= l) && (p[2] <= u);
3391}
3392
3393
3394
3395template <>
3396inline unsigned int
3397GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
3398 const unsigned int vertex)
3399{
3400 (void)line;
3401 AssertIndexRange(line, lines_per_cell);
3402 AssertIndexRange(vertex, 2);
3403
3404 return vertex;
3405}
3406
3407
3408template <>
3409inline unsigned int
3410GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
3411 const unsigned int vertex)
3412{
3413 constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
3414 return cell_vertices[line][vertex];
3415}
3416
3417
3418
3419template <>
3420inline unsigned int
3421GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
3422 const unsigned int vertex)
3423{
3424 AssertIndexRange(line, lines_per_cell);
3425 AssertIndexRange(vertex, 2);
3426
3427 constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3428 {1, 3},
3429 {0, 1},
3430 {2, 3},
3431 {4, 6}, // top face
3432 {5, 7},
3433 {4, 5},
3434 {6, 7},
3435 {0,
3436 4}, // connects of bottom
3437 {1, 5}, // top face
3438 {2, 6},
3439 {3, 7}};
3440 return vertices[line][vertex];
3441}
3442
3443
3444template <>
3445inline unsigned int
3446GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3447{
3450}
3451
3452template <>
3453inline unsigned int
3454GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3455 const bool face_orientation,
3456 const bool face_flip,
3457 const bool face_rotation)
3458{
3460
3461 // set up a table to make sure that
3462 // we handle non-standard faces correctly
3463 //
3464 // so set up a table that for each vertex (of
3465 // a quad in standard position) describes
3466 // which vertex to take
3467 //
3468 // first index: four vertices 0...3
3469 //
3470 // second index: face_orientation; 0:
3471 // opposite normal, 1: standard
3472 //
3473 // third index: face_flip; 0: standard, 1:
3474 // face rotated by 180 degrees
3475 //
3476 // forth index: face_rotation: 0: standard,
3477 // 1: face rotated by 90 degrees
3478
3479 constexpr unsigned int vertex_translation[4][2][2][2] = {
3480 {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3481 // face_rotation=false and true
3482 {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3483 // face_rotation=false and true
3484 {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3485 // face_rotation=false and true
3486 {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3487 // face_rotation=false and true
3488
3489 {{{2, 3}, // vertex 1 ...
3490 {1, 0}},
3491 {{1, 0}, {2, 3}}},
3492
3493 {{{1, 0}, // vertex 2 ...
3494 {2, 3}},
3495 {{2, 3}, {1, 0}}},
3496
3497 {{{3, 1}, // vertex 3 ...
3498 {0, 2}},
3499 {{3, 1}, {0, 2}}}};
3500
3501 return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3502}
3503
3504
3505
3506template <int dim>
3507inline unsigned int
3508GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3509 const bool,
3510 const bool,
3511 const bool)
3512{
3513 Assert(dim > 1, ExcImpossibleInDim(dim));
3515 return vertex;
3516}
3517
3518template <int dim>
3519inline unsigned int
3521{
3522 constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3523 0, 2, 2, 4, 2, 4, 4, 8};
3524
3525 return n_children[ref_case];
3526}
3527
3528
3529
3530template <int dim>
3531inline unsigned int
3533{
3535 return 0;
3536}
3537
3538template <>
3539inline unsigned int
3541{
3542 Assert(false, ExcImpossibleInDim(1));
3543 return 0;
3544}
3545
3546template <>
3547inline unsigned int
3549{
3550 return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3551}
3552
3553
3554
3555template <>
3556inline unsigned int
3558{
3559 const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3560 0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3561 return nsubs[subface_case];
3562}
3563
3564
3565
3566template <int dim>
3567inline double
3569 const unsigned int)
3570{
3572 return 0.;
3573}
3574
3575template <>
3576inline double
3578 const unsigned int)
3579{
3580 return 1;
3581}
3582
3583
3584template <>
3585inline double
3587 const unsigned int)
3588{
3589 double ratio = 1;
3590 switch (subface_case)
3591 {
3593 // Here, an
3594 // Assert(false,ExcInternalError())
3595 // would be the right
3596 // choice, but
3597 // unfortunately the
3598 // current function is
3599 // also called for faces
3600 // without children (see
3601 // tests/fe/mapping.cc).
3602 // Assert(false, ExcMessage("Face has no subfaces."));
3603 // Furthermore, assign
3604 // following value as
3605 // otherwise the
3606 // bits/volume_x tests
3607 // break
3609 break;
3611 ratio = 0.5;
3612 break;
3613 default:
3614 // there should be no
3615 // cases left
3617 break;
3618 }
3619
3620 return ratio;
3621}
3622
3623
3624template <>
3625inline double
3627 const unsigned int subface_no)
3628{
3629 double ratio = 1;
3630 switch (subface_case)
3631 {
3633 // Here, an
3634 // Assert(false,ExcInternalError())
3635 // would be the right
3636 // choice, but
3637 // unfortunately the
3638 // current function is
3639 // also called for faces
3640 // without children (see
3641 // tests/bits/mesh_3d_16.cc). Add
3642 // following switch to
3643 // avoid diffs in
3644 // tests/bits/mesh_3d_16
3646 break;
3649 ratio = 0.5;
3650 break;
3654 ratio = 0.25;
3655 break;
3658 if (subface_no < 2)
3659 ratio = 0.25;
3660 else
3661 ratio = 0.5;
3662 break;
3665 if (subface_no == 0)
3666 ratio = 0.5;
3667 else
3668 ratio = 0.25;
3669 break;
3670 default:
3671 // there should be no
3672 // cases left
3674 break;
3675 }
3676
3677 return ratio;
3678}
3679
3680
3681
3682template <int dim>
3684 const RefinementCase<dim> &,
3685 const unsigned int,
3686 const bool,
3687 const bool,
3688 const bool)
3689{
3691 return RefinementCase<dim - 1>::no_refinement;
3692}
3693
3694template <>
3696 const RefinementCase<1> &,
3697 const unsigned int,
3698 const bool,
3699 const bool,
3700 const bool)
3701{
3702 Assert(false, ExcImpossibleInDim(1));
3703
3705}
3706
3707
3708template <>
3709inline RefinementCase<1>
3711 const RefinementCase<2> &cell_refinement_case,
3712 const unsigned int face_no,
3713 const bool,
3714 const bool,
3715 const bool)
3716{
3717 const unsigned int dim = 2;
3718 AssertIndexRange(cell_refinement_case,
3721
3722 // simple special case
3723 if (cell_refinement_case == RefinementCase<dim>::cut_xy)
3725
3726 const RefinementCase<dim - 1>
3729 {RefinementCase<dim - 1>::no_refinement, // no_refinement
3730 RefinementCase<dim - 1>::no_refinement},
3731
3732 {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3733
3734 {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3735
3736 {RefinementCase<dim - 1>::cut_x, // cut_xy
3737 RefinementCase<dim - 1>::cut_x}};
3738
3739 return ref_cases[cell_refinement_case][face_no / 2];
3740}
3741
3742
3743template <>
3744inline RefinementCase<2>
3746 const RefinementCase<3> &cell_refinement_case,
3747 const unsigned int face_no,
3748 const bool face_orientation,
3749 const bool /*face_flip*/,
3750 const bool face_rotation)
3751{
3752 const unsigned int dim = 3;
3753 AssertIndexRange(cell_refinement_case,
3756
3757 // simple special case
3758 if (cell_refinement_case == RefinementCase<dim>::cut_xyz)
3759 return RefinementCase<dim - 1>::cut_xy;
3760
3761 const RefinementCase<dim - 1>
3764 {RefinementCase<dim - 1>::no_refinement, // no_refinement
3765 RefinementCase<dim - 1>::no_refinement,
3766 RefinementCase<dim - 1>::no_refinement},
3767
3768 {RefinementCase<dim - 1>::no_refinement, // cut_x
3769 RefinementCase<dim - 1>::cut_y,
3770 RefinementCase<dim - 1>::cut_x},
3771
3772 {RefinementCase<dim - 1>::cut_x, // cut_y
3773 RefinementCase<dim - 1>::no_refinement,
3774 RefinementCase<dim - 1>::cut_y},
3775
3776 {RefinementCase<dim - 1>::cut_x, // cut_xy
3777 RefinementCase<dim - 1>::cut_y,
3778 RefinementCase<dim - 1>::cut_xy},
3779
3780 {RefinementCase<dim - 1>::cut_y, // cut_z
3781 RefinementCase<dim - 1>::cut_x,
3782 RefinementCase<dim - 1>::no_refinement},
3783
3784 {RefinementCase<dim - 1>::cut_y, // cut_xz
3785 RefinementCase<dim - 1>::cut_xy,
3786 RefinementCase<dim - 1>::cut_x},
3787
3788 {RefinementCase<dim - 1>::cut_xy, // cut_yz
3789 RefinementCase<dim - 1>::cut_x,
3790 RefinementCase<dim - 1>::cut_y},
3791
3792 {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3793 RefinementCase<dim - 1>::cut_xy,
3794 RefinementCase<dim - 1>::cut_xy},
3795 };
3796
3797 const RefinementCase<dim - 1> ref_case =
3798 ref_cases[cell_refinement_case][face_no / 2];
3799
3800 const RefinementCase<dim - 1> flip[4] = {
3801 RefinementCase<dim - 1>::no_refinement,
3802 RefinementCase<dim - 1>::cut_y,
3803 RefinementCase<dim - 1>::cut_x,
3804 RefinementCase<dim - 1>::cut_xy};
3805
3806 // correct the ref_case for face_orientation
3807 // and face_rotation. for face_orientation,
3808 // 'true' is the default value whereas for
3809 // face_rotation, 'false' is standard. If
3810 // <tt>face_rotation==face_orientation</tt>,
3811 // then one of them is non-standard and we
3812 // have to swap cut_x and cut_y, otherwise no
3813 // change is necessary. face_flip has no
3814 // influence. however, in order to keep the
3815 // interface consistent with other functions,
3816 // we still include it as an argument to this
3817 // function
3818 return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3819}
3820
3821
3822
3823template <int dim>
3824inline RefinementCase<1>
3826 const unsigned int)
3827{
3830}
3831
3832template <>
3833inline RefinementCase<1>
3835 const RefinementCase<1> &cell_refinement_case,
3836 const unsigned int line_no)
3837{
3838 (void)line_no;
3839 const unsigned int dim = 1;
3840 (void)dim;
3841 AssertIndexRange(cell_refinement_case,
3844
3845 return cell_refinement_case;
3846}
3847
3848
3849template <>
3850inline RefinementCase<1>
3852 const RefinementCase<2> &cell_refinement_case,
3853 const unsigned int line_no)
3854{
3855 // Assertions are in face_refinement_case()
3856 return face_refinement_case(cell_refinement_case, line_no);
3857}
3858
3859
3860template <>
3861inline RefinementCase<1>
3863 const RefinementCase<3> &cell_refinement_case,
3864 const unsigned int line_no)
3865{
3866 const unsigned int dim = 3;
3867 AssertIndexRange(cell_refinement_case,
3870
3871 // simple special case
3872 if (cell_refinement_case == RefinementCase<dim>::cut_xyz)
3874
3875 // array indicating, which simple refine
3876 // case cuts a line in direction x, y or
3877 // z. For example, cut_y and everything
3878 // containing cut_y (cut_xy, cut_yz,
3879 // cut_xyz) cuts lines, which are in y
3880 // direction.
3881 const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3884
3885 // order the direction of lines
3886 // 0->x, 1->y, 2->z
3887 const unsigned int direction[lines_per_cell] = {
3888 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3889
3890 return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3892 RefinementCase<1>::no_refinement);
3893}
3894
3895
3896
3897template <int dim>
3901 const unsigned int,
3902 const bool,
3903 const bool,
3904 const bool)
3905{
3907
3909}
3910
3911template <>
3912inline RefinementCase<1>
3914 const RefinementCase<0> &,
3915 const unsigned int,
3916 const bool,
3917 const bool,
3918 const bool)
3919{
3920 const unsigned int dim = 1;
3921 Assert(false, ExcImpossibleInDim(dim));
3922
3924}
3925
3926
3927template <>
3928inline RefinementCase<2>
3930 const RefinementCase<1> &face_refinement_case,
3931 const unsigned int face_no,
3932 const bool,
3933 const bool,
3934 const bool)
3935{
3936 const unsigned int dim = 2;
3937 AssertIndexRange(face_refinement_case,
3940
3941 if (face_refinement_case == RefinementCase<dim>::cut_x)
3942 return (face_no / 2) != 0u ? RefinementCase<dim>::cut_x :
3944 else
3946}
3947
3948
3949template <>
3950inline RefinementCase<3>
3952 const RefinementCase<2> &face_refinement_case,
3953 const unsigned int face_no,
3954 const bool face_orientation,
3955 const bool /*face_flip*/,
3956 const bool face_rotation)
3957{
3958 const unsigned int dim = 3;
3959 AssertIndexRange(face_refinement_case,
3962
3967
3968 // correct the face_refinement_case for
3969 // face_orientation and face_rotation. for
3970 // face_orientation, 'true' is the default
3971 // value whereas for face_rotation, 'false'
3972 // is standard. If
3973 // <tt>face_rotation==face_orientation</tt>,
3974 // then one of them is non-standard and we
3975 // have to swap cut_x and cut_y, otherwise no
3976 // change is necessary. face_flip has no
3977 // influence. however, in order to keep the
3978 // interface consistent with other functions,
3979 // we still include it as an argument to this
3980 // function
3981 const RefinementCase<dim - 1> std_face_ref =
3982 (face_orientation == face_rotation) ? flip[face_refinement_case] :
3983 face_refinement_case;
3984
3985 const RefinementCase<dim> face_to_cell[3][4] = {
3986 {RefinementCase<dim>::no_refinement, // faces 0 and 1
3987 RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3990
3991 {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3992 // "exchanged on faces 2 and 3")
3996
3997 {RefinementCase<dim>::no_refinement, // faces 4 and 5
4001
4002 return face_to_cell[face_no / 2][std_face_ref];
4003}
4004
4005
4006
4007template <int dim>
4010 const unsigned int)
4011{
4013
4015}
4016
4017template <>
4018inline RefinementCase<1>
4020 const unsigned int line_no)
4021{
4022 (void)line_no;
4023 AssertIndexRange(line_no, 1);
4024
4026}
4027
4028
4029template <>
4030inline RefinementCase<2>
4032 const unsigned int line_no)
4033{
4034 const unsigned int dim = 2;
4035 (void)dim;
4037
4038 return (line_no / 2) != 0u ? RefinementCase<2>::cut_x :
4040}
4041
4042
4043template <>
4044inline RefinementCase<3>
4046 const unsigned int line_no)
4047{
4048 const unsigned int dim = 3;
4050
4051 const RefinementCase<dim> ref_cases[6] = {
4052 RefinementCase<dim>::cut_y, // lines 0 and 1
4053 RefinementCase<dim>::cut_x, // lines 2 and 3
4054 RefinementCase<dim>::cut_y, // lines 4 and 5
4055 RefinementCase<dim>::cut_x, // lines 6 and 7
4056 RefinementCase<dim>::cut_z, // lines 8 and 9
4057 RefinementCase<dim>::cut_z}; // lines 10 and 11
4058
4059 return ref_cases[line_no / 2];
4060}
4061
4062
4063
4064template <>
4065inline unsigned int
4066GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
4067 const bool face_orientation,
4068 const bool face_flip,
4069 const bool face_rotation)
4070{
4072
4073 // set up a table to make sure that
4074 // we handle non-standard faces correctly
4075 //
4076 // so set up a table that for each vertex (of
4077 // a quad in standard position) describes
4078 // which vertex to take
4079 //
4080 // first index: four vertices 0...3
4081 //
4082 // second index: face_orientation; 0:
4083 // opposite normal, 1: standard
4084 //
4085 // third index: face_flip; 0: standard, 1:
4086 // face rotated by 180 degrees
4087 //
4088 // forth index: face_rotation: 0: standard,
4089 // 1: face rotated by 90 degrees
4090
4091 const unsigned int vertex_translation[4][2][2][2] = {
4092 {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
4093 // face_rotation=false and true
4094 {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
4095 // face_rotation=false and true
4096 {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
4097 // face_rotation=false and true
4098 {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
4099 // face_rotation=false and true
4100
4101 {{{2, 3}, // vertex 1 ...
4102 {1, 0}},
4103 {{1, 3}, {2, 0}}},
4104
4105 {{{1, 0}, // vertex 2 ...
4106 {2, 3}},
4107 {{2, 0}, {1, 3}}},
4108
4109 {{{3, 1}, // vertex 3 ...
4110 {0, 2}},
4111 {{3, 2}, {0, 1}}}};
4112
4113 return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
4114}
4115
4116
4117
4118template <int dim>
4119inline unsigned int
4120GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
4121 const bool,
4122 const bool,
4123 const bool)
4124{
4125 Assert(dim > 1, ExcImpossibleInDim(dim));
4127 return vertex;
4128}
4129
4130
4131
4132template <>
4133inline unsigned int
4134GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
4135 const bool face_orientation,
4136 const bool face_flip,
4137 const bool face_rotation)
4138{
4140
4141
4142 // make sure we handle
4143 // non-standard faces correctly
4144 //
4145 // so set up a table that for each line (of a
4146 // quad) describes which line to take
4147 //
4148 // first index: four lines 0...3
4149 //
4150 // second index: face_orientation; 0:
4151 // opposite normal, 1: standard
4152 //
4153 // third index: face_flip; 0: standard, 1:
4154 // face rotated by 180 degrees
4155 //
4156 // forth index: face_rotation: 0: standard,
4157 // 1: face rotated by 90 degrees
4158
4159 const unsigned int line_translation[4][2][2][2] = {
4160 {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4161 // face_rotation=false and true
4162 {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4163 // face_rotation=false and true
4164 {{0, 3}, // line 0, face_orientation=true, face_flip=false,
4165 // face_rotation=false and true
4166 {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
4167 // face_rotation=false and true
4168
4169 {{{3, 1}, // line 1 ...
4170 {2, 0}},
4171 {{1, 2}, {0, 3}}},
4172
4173 {{{0, 3}, // line 2 ...
4174 {1, 2}},
4175 {{2, 0}, {3, 1}}},
4176
4177 {{{1, 2}, // line 3 ...
4178 {0, 3}},
4179 {{3, 1}, {2, 0}}}};
4180
4181 return line_translation[line][face_orientation][face_flip][face_rotation];
4182}
4183
4184
4185
4186template <int dim>
4187inline unsigned int
4188GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
4189 const bool,
4190 const bool,
4191 const bool)
4192{
4194 return line;
4195}
4196
4197
4198
4199template <>
4200inline unsigned int
4201GeometryInfo<2>::standard_to_real_line_vertex(const unsigned int vertex,
4202 const bool line_orientation)
4203{
4204 return line_orientation ? vertex : (1 - vertex);
4205}
4206
4207
4208
4209template <int dim>
4210inline unsigned int
4211GeometryInfo<dim>::standard_to_real_line_vertex(const unsigned int vertex,
4212 const bool)
4213{
4215 return vertex;
4216}
4217
4218
4219
4220template <>
4221inline std::array<unsigned int, 2>
4223 const unsigned int vertex)
4224{
4225 return {{vertex % 2, vertex / 2}};
4226}
4227
4228
4229
4230template <int dim>
4231inline std::array<unsigned int, 2>
4233 const unsigned int vertex)
4234{
4236 (void)vertex;
4237 return {{0, 0}};
4238}
4239
4240
4241
4242template <>
4243inline std::array<unsigned int, 2>
4245{
4246 // set up a table that for each
4247 // line describes a) from which
4248 // quad to take it, b) which line
4249 // therein it is if the face is
4250 // oriented correctly
4251 static const unsigned int lookup_table[GeometryInfo<3>::lines_per_cell][2] = {
4252 {4, 0}, // take first four lines from bottom face
4253 {4, 1},
4254 {4, 2},
4255 {4, 3},
4256
4257 {5, 0}, // second four lines from top face
4258 {5, 1},
4259 {5, 2},
4260 {5, 3},
4261
4262 {0, 0}, // the rest randomly
4263 {1, 0},
4264 {0, 1},
4265 {1, 1}};
4266
4267 return {{lookup_table[i][0], lookup_table[i][1]}};
4268}
4269
4270
4271
4272template <int dim>
4273inline std::array<unsigned int, 2>
4275{
4277 (void)line;
4278 return {{0, 0}};
4279}
4280
4281
4282
4283template <>
4284inline std::array<unsigned int, 2>
4286 const unsigned int vertex)
4287{
4288 // get the corner indices by asking either the bottom or the top face for its
4289 // vertices. handle non-standard faces by calling the vertex reordering
4290 // function from GeometryInfo
4291
4292 // bottom face (4) for first four vertices, top face (5) for the rest
4293 return {{4 + vertex / 4, vertex % 4}};
4294}
4295
4296
4297
4298template <int dim>
4299inline std::array<unsigned int, 2>
4301 const unsigned int vertex)
4302{
4304 (void)vertex;
4305 return {{0, 0}};
4306}
4307
4308
4309
4310template <>
4311inline unsigned int
4312GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
4313 const bool face_orientation,
4314 const bool face_flip,
4315 const bool face_rotation)
4316{
4318
4319
4320 // make sure we handle
4321 // non-standard faces correctly
4322 //
4323 // so set up a table that for each line (of a
4324 // quad) describes which line to take
4325 //
4326 // first index: four lines 0...3
4327 //
4328 // second index: face_orientation; 0:
4329 // opposite normal, 1: standard
4330 //
4331 // third index: face_flip; 0: standard, 1:
4332 // face rotated by 180 degrees
4333 //
4334 // forth index: face_rotation: 0: standard,
4335 // 1: face rotated by 90 degrees
4336
4337 const unsigned int line_translation[4][2][2][2] = {
4338 {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4339 // face_rotation=false and true
4340 {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4341 // face_rotation=false and true
4342 {{0, 2}, // line 0, face_orientation=true, face_flip=false,
4343 // face_rotation=false and true
4344 {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
4345 // face_rotation=false and true
4346
4347 {{{3, 1}, // line 1 ...
4348 {2, 0}},
4349 {{1, 3}, {0, 2}}},
4350
4351 {{{0, 3}, // line 2 ...
4352 {1, 2}},
4353 {{2, 1}, {3, 0}}},
4354
4355 {{{1, 2}, // line 3 ...
4356 {0, 3}},
4357 {{3, 0}, {2, 1}}}};
4358
4359 return line_translation[line][face_orientation][face_flip][face_rotation];
4360}
4361
4362
4363
4364template <int dim>
4365inline unsigned int
4366GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
4367 const bool,
4368 const bool,
4369 const bool)
4370{
4372 return line;
4373}
4374
4375
4376
4377template <>
4378inline unsigned int
4380 const unsigned int face,
4381 const unsigned int subface,
4382 const bool,
4383 const bool,
4384 const bool,
4385 const RefinementCase<0> &)
4386{
4387 (void)subface;
4388 AssertIndexRange(face, faces_per_cell);
4389 AssertIndexRange(subface, max_children_per_face);
4390
4391 return face;
4392}
4393
4394
4395
4396template <>
4397inline unsigned int
4399 const unsigned int face,
4400 const unsigned int subface,
4401 const bool face_orientation,
4402 const bool /*face_flip*/,
4403 const bool /*face_rotation*/,
4404 const RefinementCase<1> &)
4405{
4406 AssertIndexRange(face, faces_per_cell);
4407 AssertIndexRange(subface, max_children_per_face);
4408
4409 // always return the child adjacent to the specified
4410 // subface. if the face of a cell is not refined, don't
4411 // throw an assertion but deliver the child adjacent to
4412 // the face nevertheless, i.e. deliver the child of
4413 // this cell adjacent to the subface of a possibly
4414 // refined neighbor. this simplifies setting neighbor
4415 // information in execute_refinement.
4416 constexpr unsigned int
4417 subcells[/* possible face orientation */ 2]
4418 [/* number of different ways to refine a cell */ 4]
4419 [/* faces_per_cell */ 4][/* max_children_per_face */ 2] = {
4420 {
4421 // Flipped orientation (face_orientation = false)
4422 {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
4423 {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
4424 {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
4425 },
4426 {
4427 // Normal orientation (face_orientation = true)
4428 {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
4429 {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
4430 {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
4431 }};
4432
4433 return subcells[face_orientation][ref_case - 1][face][subface];
4434}
4435
4436
4437
4438template <>
4439inline unsigned int
4441 const unsigned int face,
4442 const unsigned int subface,
4443 const bool face_orientation,
4444 const bool face_flip,
4445 const bool face_rotation,
4446 const RefinementCase<2> &face_ref_case)
4447{
4448 const unsigned int dim = 3;
4449
4451 ExcMessage("Cell has no children."));
4452 AssertIndexRange(face, faces_per_cell);
4453 if (!(subface == 0 &&
4455 {
4456 AssertIndexRange(subface,
4457 GeometryInfo<dim - 1>::n_children(face_ref_case));
4458 }
4459
4460 // invalid number used for invalid cases,
4461 // e.g. when the children are more refined at
4462 // a given face than the face itself
4463 const unsigned int e = numbers::invalid_unsigned_int;
4464
4465 // the whole process of finding a child cell
4466 // at a given subface considering the
4467 // possibly anisotropic refinement cases of
4468 // the cell and the face as well as
4469 // orientation, flip and rotation of the face
4470 // is quite complicated. thus, we break it
4471 // down into several steps.
4472
4473 // first step: convert the given face refine
4474 // case to a face refine case concerning the
4475 // face in standard orientation (, flip and
4476 // rotation). This only affects cut_x and
4477 // cut_y
4478 const RefinementCase<dim - 1> flip[4] = {
4479 RefinementCase<dim - 1>::no_refinement,
4480 RefinementCase<dim - 1>::cut_y,
4481 RefinementCase<dim - 1>::cut_x,
4482 RefinementCase<dim - 1>::cut_xy};
4483 // for face_orientation, 'true' is the
4484 // default value whereas for face_rotation,
4485 // 'false' is standard. If
4486 // <tt>face_rotation==face_orientation</tt>,
4487 // then one of them is non-standard and we
4488 // have to swap cut_x and cut_y, otherwise no
4489 // change is necessary.
4490 const RefinementCase<dim - 1> std_face_ref =
4491 (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4492
4493 // second step: convert the given subface
4494 // index to the one for a standard face
4495 // respecting face_orientation, face_flip and
4496 // face_rotation
4497
4498 // first index: face_ref_case
4499 // second index: face_orientation
4500 // third index: face_flip
4501 // forth index: face_rotation
4502 // fifth index: subface index
4503 const unsigned int subface_exchange[4][2][2][2][4] = {
4504 // no_refinement (subface 0 stays 0,
4505 // all others are invalid)
4506 {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4507 {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4508 // cut_x (here, if the face is only
4509 // rotated OR only falsely oriented,
4510 // then subface 0 of the non-standard
4511 // face does NOT correspond to one of
4512 // the subfaces of a standard
4513 // face. Thus we indicate the subface
4514 // which is located at the lower left
4515 // corner (the origin of the face's
4516 // local coordinate system) with
4517 // '0'. The rest of this issue is
4518 // taken care of using the above
4519 // conversion to a 'standard face
4520 // refine case')
4521 {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4522 {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4523 // cut_y (the same applies as for
4524 // cut_x)
4525 {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4526 {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4527 // cut_xyz: this information is
4528 // identical to the information
4529 // returned by
4530 // GeometryInfo<3>::real_to_standard_face_vertex()
4531 {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4532 // face_rotation=false, subfaces 0,1,2,3
4533 {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4534 // face_rotation=true, subfaces 0,1,2,3
4535 {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4536 // face_rotation=false, subfaces 0,1,2,3
4537 {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4538 // face_rotation=true, subfaces 0,1,2,3
4539 {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4540 // face_rotation=false, subfaces 0,1,2,3
4541 {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4542 // face_rotation=true, subfaces 0,1,2,3
4543 {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4544 // face_rotation=false, subfaces 0,1,2,3
4545 {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4546 // face_rotation=true, subfaces 0,1,2,3
4547
4548 const unsigned int std_subface =
4549 subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4550 [subface];
4551 Assert(std_subface != e, ExcInternalError());
4552
4553 // third step: these are the children, which
4554 // can be found at the given subfaces of an
4555 // isotropically refined (standard) face
4556 //
4557 // first index: (refinement_case-1)
4558 // second index: face_index
4559 // third index: subface_index (isotropic refinement)
4560 const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4561 [max_children_per_face] = {
4562 // cut_x
4563 {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4564 {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4565 {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4566 {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4567 {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4568 {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4569 // cut_y
4570 {{0, 1, 0, 1},
4571 {0, 1, 0, 1},
4572 {0, 0, 0, 0},
4573 {1, 1, 1, 1},
4574 {0, 0, 1, 1},
4575 {0, 0, 1, 1}},
4576 // cut_xy
4577 {{0, 2, 0, 2},
4578 {1, 3, 1, 3},
4579 {0, 0, 1, 1},
4580 {2, 2, 3, 3},
4581 {0, 1, 2, 3},
4582 {0, 1, 2, 3}},
4583 // cut_z
4584 {{0, 0, 1, 1},
4585 {0, 0, 1, 1},
4586 {0, 1, 0, 1},
4587 {0, 1, 0, 1},
4588 {0, 0, 0, 0},
4589 {1, 1, 1, 1}},
4590 // cut_xz
4591 {{0, 0, 1, 1},
4592 {2, 2, 3, 3},
4593 {0, 1, 2, 3},
4594 {0, 1, 2, 3},
4595 {0, 2, 0, 2},
4596 {1, 3, 1, 3}},
4597 // cut_yz
4598 {{0, 1, 2, 3},
4599 {0, 1, 2, 3},
4600 {0, 2, 0, 2},
4601 {1, 3, 1, 3},
4602 {0, 0, 1, 1},
4603 {2, 2, 3, 3}},
4604 // cut_xyz
4605 {{0, 2, 4, 6},
4606 {1, 3, 5, 7},
4607 {0, 4, 1, 5},
4608 {2, 6, 3, 7},
4609 {0, 1, 2, 3},
4610 {4, 5, 6, 7}}};
4611
4612 // forth step: check, whether the given face
4613 // refine case is valid for the given cell
4614 // refine case. this is the case, if the
4615 // given face refine case is at least as
4616 // refined as the face is for the given cell
4617 // refine case
4618
4619 // note, that we are considering standard
4620 // face refinement cases here and thus must
4621 // not pass the given orientation, flip and
4622 // rotation flags
4623 if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4624 face_refinement_case(ref_case, face))
4625 {
4626 // all is fine. for anisotropic face
4627 // refine cases, select one of the
4628 // isotropic subfaces which neighbors the
4629 // same child
4630
4631 // first index: (standard) face refine case
4632 // second index: subface index
4633 const unsigned int equivalent_iso_subface[4][4] = {
4634 {0, e, e, e}, // no_refinement
4635 {0, 3, e, e}, // cut_x
4636 {0, 3, e, e}, // cut_y
4637 {0, 1, 2, 3}}; // cut_xy
4638
4639 const unsigned int equ_std_subface =
4640 equivalent_iso_subface[std_face_ref][std_subface];
4641 Assert(equ_std_subface != e, ExcInternalError());
4642
4643 return iso_children[ref_case - 1][face][equ_std_subface];
4644 }
4645 else
4646 {
4647 // the face_ref_case was too coarse,
4648 // throw an error
4649 Assert(false,
4650 ExcMessage("The face RefineCase is too coarse "
4651 "for the given cell RefineCase."));
4652 }
4653 // we only get here in case of an error
4654 return e;
4655}
4656
4657
4658
4659template <>
4660inline unsigned int
4662 const unsigned int,
4663 const unsigned int,
4664 const bool,
4665 const bool,
4666 const bool,
4667 const RefinementCase<3> &)
4668{
4671}
4672
4673
4674
4675template <>
4676inline unsigned int
4677GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4678 const unsigned int line,
4679 const bool,
4680 const bool,
4681 const bool)
4682{
4683 (void)line;
4684 AssertIndexRange(face, faces_per_cell);
4685 AssertIndexRange(line, lines_per_face);
4686
4687 // The face is a line itself.
4688 return face;
4689}
4690
4691
4692
4693template <>
4694inline unsigned int
4695GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4696 const unsigned int line,
4697 const bool face_orientation,
4698 const bool face_flip,
4699 const bool face_rotation)
4700{
4701 AssertIndexRange(face, faces_per_cell);
4702 AssertIndexRange(line, lines_per_face);
4703
4704 const unsigned lines[faces_per_cell][lines_per_face] = {
4705 {8, 10, 0, 4}, // left face
4706 {9, 11, 1, 5}, // right face
4707 {2, 6, 8, 9}, // front face
4708 {3, 7, 10, 11}, // back face
4709 {0, 1, 2, 3}, // bottom face
4710 {4, 5, 6, 7}}; // top face
4711 return lines[face][real_to_standard_face_line(
4712 line, face_orientation, face_flip, face_rotation)];
4713}
4714
4715
4716
4717inline unsigned int
4718GeometryInfo<0>::face_to_cell_lines(const unsigned int,
4719 const unsigned int,
4720 const bool,
4721 const bool,
4722 const bool)
4723{
4726}
4727
4728
4729
4730template <int dim>
4731inline unsigned int
4732GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4733 const unsigned int,
4734 const bool,
4735 const bool,
4736 const bool)
4737{
4740}
4741
4742
4743
4744template <int dim>
4745inline unsigned int
4746GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4747 const unsigned int vertex,
4748 const bool face_orientation,
4749 const bool face_flip,
4750 const bool face_rotation)
4751{
4752 return child_cell_on_face(RefinementCase<dim>::isotropic_refinement,
4753 face,
4754 vertex,
4755 face_orientation,
4756 face_flip,
4757 face_rotation);
4758}
4759
4760
4761
4762inline unsigned int
4764 const unsigned int,
4765 const bool,
4766 const bool,
4767 const bool)
4768{
4771}
4772
4773
4774
4775template <int dim>
4776template <typename Number>
4777inline Point<dim, Number>
4779{
4781 for (unsigned int i = 0; i < dim; ++i)
4782 p[i] = std::min(std::max(q[i], Number(0.)), Number(1.));
4783
4784 return p;
4785}
4786
4787
4788
4789template <int dim>
4790inline double
4792{
4793 double result = 0.0;
4794
4795 for (unsigned int i = 0; i < dim; ++i)
4796 {
4797 result = std::max(result, -p[i]);
4798 result = std::max(result, p[i] - 1.);
4799 }
4800
4801 return result;
4802}
4803
4804
4805
4806template <int dim>
4807inline double
4809 const unsigned int i)
4810{
4812
4813 switch (dim)
4814 {
4815 case 1:
4816 {
4817 const double x = xi[0];
4818 switch (i)
4819 {
4820 case 0:
4821 return 1 - x;
4822 case 1:
4823 return x;
4824 }
4825 break;
4826 }
4827
4828 case 2:
4829 {
4830 const double x = xi[0];
4831 const double y = xi[1];
4832 switch (i)
4833 {
4834 case 0:
4835 return (1 - x) * (1 - y);
4836 case 1:
4837 return x * (1 - y);
4838 case 2:
4839 return (1 - x) * y;
4840 case 3:
4841 return x * y;
4842 }
4843 break;
4844 }
4845
4846 case 3:
4847 {
4848 const double x = xi[0];
4849 const double y = xi[1];
4850 const double z = xi[2];
4851 switch (i)
4852 {
4853 case 0:
4854 return (1 - x) * (1 - y) * (1 - z);
4855 case 1:
4856 return x * (1 - y) * (1 - z);
4857 case 2:
4858 return (1 - x) * y * (1 - z);
4859 case 3:
4860 return x * y * (1 - z);
4861 case 4:
4862 return (1 - x) * (1 - y) * z;
4863 case 5:
4864 return x * (1 - y) * z;
4865 case 6:
4866 return (1 - x) * y * z;
4867 case 7:
4868 return x * y * z;
4869 }
4870 break;
4871 }
4872
4873 default:
4875 }
4876 return -1e9;
4877}
4878
4879
4880
4881template <>
4883 const Point<1> &,
4884 const unsigned int i)
4885{
4887
4888 switch (i)
4889 {
4890 case 0:
4891 return Point<1>(-1.);
4892 case 1:
4893 return Point<1>(1.);
4894 }
4895
4896 return Point<1>(-1e9);
4897}
4898
4899
4900
4901template <>
4903 const Point<2> &xi,
4904 const unsigned int i)
4905{
4907
4908 const double x = xi[0];
4909 const double y = xi[1];
4910 switch (i)
4911 {
4912 case 0:
4913 return Point<2>(-(1 - y), -(1 - x));
4914 case 1:
4915 return Point<2>(1 - y, -x);
4916 case 2:
4917 return Point<2>(-y, 1 - x);
4918 case 3:
4919 return Point<2>(y, x);
4920 }
4921 return Point<2>(-1e9, -1e9);
4922}
4923
4924
4925
4926template <>
4928 const Point<3> &xi,
4929 const unsigned int i)
4930{
4932
4933 const double x = xi[0];
4934 const double y = xi[1];
4935 const double z = xi[2];
4936 switch (i)
4937 {
4938 case 0:
4939 return Point<3>(-(1 - y) * (1 - z),
4940 -(1 - x) * (1 - z),
4941 -(1 - x) * (1 - y));
4942 case 1:
4943 return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4944 case 2:
4945 return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4946 case 3:
4947 return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4948 case 4:
4949 return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4950 case 5:
4951 return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4952 case 6:
4953 return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4954 case 7:
4955 return Point<3>(y * z, x * z, x * y);
4956 }
4957
4958 return Point<3>(-1e9, -1e9, -1e9);
4959}
4960
4961
4962
4963template <int dim>
4964inline Tensor<1, dim>
4966 const unsigned int)
4967{
4969 return Tensor<1, dim>();
4970}
4971
4972
4973unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4974{
4975 return 0;
4976}
4977
4978
4979namespace internal
4980{
4981 namespace GeometryInfoHelper
4982 {
4983 // wedge product of a single
4984 // vector in 2d: we just have to
4985 // rotate it by 90 degrees to the
4986 // right
4987 inline Tensor<1, 2>
4988 wedge_product(const Tensor<1, 2> (&derivative)[1])
4989 {
4990 Tensor<1, 2> result;
4991 result[0] = derivative[0][1];
4992 result[1] = -derivative[0][0];
4993
4994 return result;
4995 }
4996
4997
4998 // wedge product of 2 vectors in
4999 // 3d is the cross product
5000 inline Tensor<1, 3>
5001 wedge_product(const Tensor<1, 3> (&derivative)[2])
5002 {
5003 return cross_product_3d(derivative[0], derivative[1]);
5004 }
5005
5006
5007 // wedge product of dim vectors
5008 // in dim-d: that's the
5009 // determinant of the matrix
5010 template <int dim>
5011 inline Tensor<0, dim>
5012 wedge_product(const Tensor<1, dim> (&derivative)[dim])
5013 {
5014 Tensor<2, dim> jacobian;
5015 for (unsigned int i = 0; i < dim; ++i)
5016 jacobian[i] = derivative[i];
5017
5018 return determinant(jacobian);
5019 }
5020 } // namespace GeometryInfoHelper
5021} // namespace internal
5022
5023
5024template <int dim>
5025template <int spacedim>
5026inline void
5028# ifndef DEAL_II_CXX14_CONSTEXPR_BUG
5029 (const Point<spacedim> (&vertices)[vertices_per_cell],
5030 Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell])
5031# else
5032 (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms)
5033# endif
5034{
5035 // for each of the vertices,
5036 // compute the alternating form
5037 // of the mapped unit
5038 // vectors. consider for
5039 // example the case of a quad
5040 // in spacedim==3: to do so, we
5041 // need to see how the
5042 // infinitesimal vectors
5043 // (d\xi_1,0) and (0,d\xi_2)
5044 // are transformed into
5045 // spacedim-dimensional space
5046 // and then form their cross
5047 // product (i.e. the wedge product
5048 // of two vectors). to this end, note
5049 // that
5050 // \vec x = sum_i \vec v_i phi_i(\vec xi)
5051 // so the transformed vectors are
5052 // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
5053 // and
5054 // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
5055 // which boils down to the columns
5056 // of the 3x2 matrix \grad_\xi x(\xi)
5057 //
5058 // a similar reasoning would
5059 // hold for all dim,spacedim
5060 // pairs -- we only have to
5061 // compute the wedge product of
5062 // the columns of the
5063 // derivatives
5064 for (unsigned int i = 0; i < vertices_per_cell; ++i)
5065 {
5066 Tensor<1, spacedim> derivatives[dim];
5067
5068 for (unsigned int j = 0; j < vertices_per_cell; ++j)
5069 {
5070 const Tensor<1, dim> grad_phi_j =
5071 d_linear_shape_function_gradient(unit_cell_vertex(i), j);
5072 for (unsigned int l = 0; l < dim; ++l)
5073 derivatives[l] += vertices[j] * grad_phi_j[l];
5074 }
5075
5076 forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
5077 }
5078}
5079
5080#endif // DOXYGEN
5082
5083#endif
GeometryPrimitive(const unsigned int object_dimension)
GeometryPrimitive(const Object object)
Definition point.h:111
constexpr RefinementCase()
constexpr RefinementCase(const typename RefinementPossibilities< dim >::Possibilities refinement_case)
constexpr RefinementCase operator&(const RefinementCase &r) const
constexpr RefinementCase operator~() const
static constexpr unsigned int n_refinement_cases
std::uint8_t value
static constexpr RefinementCase cut_axis(const unsigned int i)
void serialize(Archive &ar, const unsigned int version)
constexpr RefinementCase operator|(const RefinementCase &r) const
constexpr RefinementCase(const std::uint8_t refinement_case)
static constexpr std::array< RefinementCase< dim >, n_refinement_cases > all_refinement_cases()
static constexpr std::size_t memory_consumption()
static constexpr std::size_t memory_consumption()
SubfaceCase(const typename SubfacePossibilities< dim >::Possibilities subface_possibility)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
#define DEAL_II_HOST_DEVICE
Definition config.h:171
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static ::ExceptionBase & ExcInvalidCoordinate(double arg1)
static ::ExceptionBase & ExcInvalidSubface(int arg1, int arg2, int arg3)
static ::ExceptionBase & ExcInvalidRefinementCase(int arg1)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
#define DeclException3(Exception3, type1, type2, type3, outsequence)
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcInvalidSubfaceCase(int arg1)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:210
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
*  *  *  RotationFunction< dim, Number >::RotationFunction Number(dim)
constexpr unsigned int invalid_unsigned_int
Definition types.h:228
boost::integer_range< IncrementableType > iota_view
Definition iota_view.h:43
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:105
static unsigned int n_children(const RefinementCase< 0 > &refinement_case)
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, 0 > face_indices()
static std::array< unsigned int, vertices_per_cell > vertex_indices()
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 unsigned int child_cell_from_point(const Point< dim > &p)
static unsigned int standard_to_real_line_vertex(const unsigned int vertex, const bool line_orientation=true)
static constexpr std::array< unsigned int, vertices_per_cell > ucd_to_deal
static constexpr ndarray< Tensor< 1, dim >, faces_per_cell, dim - 1 > unit_tangential_vectors
static bool is_inside_unit_cell(const Point< dim > &p, const double eps)
static double distance_to_unit_cell(const Point< dim > &p)
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 constexpr unsigned int quads_per_face
static RefinementCase< dim - 1 > face_refinement_case(const RefinementCase< dim > &cell_refinement_case, const unsigned int face_no, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static unsigned int standard_to_real_face_vertex(const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static RefinementCase< dim > min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
static constexpr std::array< unsigned int, faces_per_cell > opposite_face
static constexpr unsigned int max_children_per_face
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 constexpr unsigned int vertices_per_cell
static constexpr std::array< unsigned int, vertices_per_cell > dx_to_deal
static constexpr unsigned int lines_per_cell
static constexpr std::array< unsigned int, faces_per_cell > unit_normal_direction
static std::array< unsigned int, 2 > standard_hex_line_to_quad_line_index(const unsigned int line)
static constexpr unsigned int faces_per_cell
static unsigned int standard_to_real_face_line(const unsigned int line, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
static constexpr unsigned int hexes_per_cell
static Point< dim, Number > project_to_unit_cell(const Point< dim, Number > &p)
static RefinementCase< 1 > line_refinement_case(const RefinementCase< dim > &cell_refinement_case, const unsigned int line_no)
static std::array< unsigned int, 2 > standard_hex_vertex_to_quad_vertex_index(const unsigned int vertex)
static std::array< unsigned int, 2 > standard_quad_vertex_to_line_vertex_index(const unsigned int vertex)
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
static constexpr std::array< Tensor< 1, dim >, faces_per_cell > unit_normal_vector
static constexpr unsigned int vertices_per_face
static unsigned int real_to_standard_face_line(const unsigned int line, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
static bool is_inside_unit_cell(const Point< dim > &p)
static RefinementCase< dim > min_cell_refinement_case_for_face_refinement(const RefinementCase< dim - 1 > &face_refinement_case, const unsigned int face_no, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
static unsigned int real_to_standard_face_vertex(const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static unsigned int face_to_cell_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 Point< dim > cell_to_child_coordinates(const Point< dim > &p, const unsigned int child_index, const RefinementCase< dim > refine_case=RefinementCase< dim >::isotropic_refinement)
static constexpr ndarray< unsigned int, vertices_per_cell, dim > vertex_to_face
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
static void alternating_form_at_vertices(const Point< spacedim >(&vertices)[vertices_per_cell], Tensor< spacedim - dim, spacedim >(&forms)[vertices_per_cell])
static Point< dim > unit_cell_vertex(const unsigned int vertex)
static constexpr std::array< int, faces_per_cell > unit_normal_orientation
static constexpr unsigned int quads_per_cell
static Point< dim > child_to_cell_coordinates(const Point< dim > &p, const unsigned int child_index, const RefinementCase< dim > refine_case=RefinementCase< dim >::isotropic_refinement)
static constexpr unsigned int max_children_per_cell
static unsigned int n_subfaces(const internal::SubfaceCase< dim > &subface_case)
static constexpr unsigned int lines_per_face
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)