Reference documentation for deal.II version GIT relicensing-422-gb369f187d8 2024-04-17 18:10:02+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: LGPL-2.1-or-later
4// Copyright (C) 1998 - 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_geometry_info_h
16#define dealii_geometry_info_h
17
18
19#include <deal.II/base/config.h>
20
23#include <deal.II/base/point.h>
25
26#include <array>
27#include <cstdint>
28
29
30
32
33#ifndef DOXYGEN
34namespace internal
35{
36 namespace GeometryInfoHelper
37 {
38 // A struct that holds the values for all the arrays we want to initialize
39 // in GeometryInfo
40 template <int dim>
41 struct Initializers;
42
43 template <>
44 struct Initializers<1>
45 {
46 static constexpr std::array<unsigned int, 2>
47 ucd_to_deal()
48 {
49 return {{0, 1}};
50 }
51
52 static constexpr std::array<unsigned int, 2>
53 unit_normal_direction()
54 {
55 return {{0, 0}};
56 }
57
58 static constexpr std::array<int, 2>
59 unit_normal_orientation()
60 {
61 return {{-1, 1}};
62 }
63
64 static constexpr std::array<Tensor<1, 1>, 2>
65 unit_normal_vector()
66 {
67 return {{Tensor<1, 1>{{-1}}, Tensor<1, 1>{{1}}}};
68 }
69
70 static constexpr ::ndarray<Tensor<1, 1>, 2, 0>
71 unit_tangential_vectors()
72 {
73 return {{{{}}, {{}}}};
74 }
75
76 static constexpr std::array<unsigned int, 2>
77 opposite_face()
78 {
79 return {{1, 0}};
80 }
81
82 static constexpr std::array<unsigned int, 2>
83 dx_to_deal()
84 {
85 return {{0, 1}};
86 }
87
88 static constexpr ::ndarray<unsigned int, 2, 1>
89 vertex_to_face()
90 {
91 return {{{{0}}, {{1}}}};
92 }
93 };
94
95 template <>
96 struct Initializers<2>
97 {
98 static constexpr std::array<unsigned int, 4>
99 ucd_to_deal()
100 {
101 return {{0, 1, 3, 2}};
102 }
103
104 static constexpr std::array<unsigned int, 4>
105 unit_normal_direction()
106 {
107 return {{0, 0, 1, 1}};
108 }
109
110 static constexpr std::array<int, 4>
111 unit_normal_orientation()
112 {
113 return {{-1, 1, -1, 1}};
114 }
115
116 static constexpr std::array<Tensor<1, 2>, 4>
117 unit_normal_vector()
118 {
119 return {{Tensor<1, 2>{{-1., 0.}},
120 Tensor<1, 2>{{1., 0.}},
121 Tensor<1, 2>{{0., -1.}},
122 Tensor<1, 2>{{0., 1.}}}};
123 }
124
125 static constexpr ::ndarray<Tensor<1, 2>, 4, 1>
126 unit_tangential_vectors()
127 {
128 return {{{{Tensor<1, 2>{{0, -1}}}},
129 {{Tensor<1, 2>{{0, 1}}}},
130 {{Tensor<1, 2>{{1, 0}}}},
131 {{Tensor<1, 2>{{-1, 0}}}}}};
132 }
133
134 static constexpr std::array<unsigned int, 4>
135 opposite_face()
136 {
137 return {{1, 0, 3, 2}};
138 }
139
140 static constexpr std::array<unsigned int, 4>
141 dx_to_deal()
142 {
143 return {{0, 2, 1, 3}};
144 }
145
146 static constexpr ::ndarray<unsigned int, 4, 2>
147 vertex_to_face()
148 {
149 return {{{{0, 2}}, {{1, 2}}, {{0, 3}}, {{1, 3}}}};
150 }
151 };
152
153 template <>
154 struct Initializers<3>
155 {
156 static constexpr std::array<unsigned int, 8>
157 ucd_to_deal()
158 {
159 return {{0, 4, 5, 1, 2, 6, 7, 3}};
160 }
161
162 static constexpr std::array<unsigned int, 6>
163 unit_normal_direction()
164 {
165 return {{0, 0, 1, 1, 2, 2}};
166 }
167
168 static constexpr std::array<int, 6>
169 unit_normal_orientation()
170 {
171 return {{-1, 1, -1, 1, -1, 1}};
172 }
173
174 static constexpr std::array<Tensor<1, 3>, 6>
175 unit_normal_vector()
176 {
177 return {{Tensor<1, 3>{{-1, 0, 0}},
178 Tensor<1, 3>{{1, 0, 0}},
179 Tensor<1, 3>{{0, -1, 0}},
180 Tensor<1, 3>{{0, 1, 0}},
181 Tensor<1, 3>{{0, 0, -1}},
182 Tensor<1, 3>{{0, 0, 1}}}};
183 }
184
185 static constexpr ::ndarray<Tensor<1, 3>, 6, 2>
186 unit_tangential_vectors()
187 {
188 return {{{{Tensor<1, 3>{{0, -1, 0}}, Tensor<1, 3>{{0, 0, 1}}}},
189 {{Tensor<1, 3>{{0, 1, 0}}, Tensor<1, 3>{{0, 0, 1}}}},
190 {{Tensor<1, 3>{{0, 0, -1}}, Tensor<1, 3>{{1, 0, 0}}}},
191 {{Tensor<1, 3>{{0, 0, 1}}, Tensor<1, 3>{{1, 0, 0}}}},
192 {{Tensor<1, 3>{{-1, 0, 0}}, Tensor<1, 3>{{0, 1, 0}}}},
193 {{Tensor<1, 3>{{1, 0, 0}}, Tensor<1, 3>{{0, 1, 0}}}}}};
194 }
195
196 static constexpr std::array<unsigned int, 6>
197 opposite_face()
198 {
199 return {{1, 0, 3, 2, 5, 4}};
200 }
201
202 static constexpr std::array<unsigned int, 8>
203 dx_to_deal()
204 {
205 return {{0, 4, 2, 6, 1, 5, 3, 7}};
206 }
207
208 static constexpr ::ndarray<unsigned int, 8, 3>
209 vertex_to_face()
210 {
211 return {{{{0, 2, 4}},
212 {{1, 2, 4}},
213 {{0, 3, 4}},
214 {{1, 3, 4}},
215 {{0, 2, 5}},
216 {{1, 2, 5}},
217 {{0, 3, 5}},
218 {{1, 3, 5}}}};
219 }
220 };
221
222 template <>
223 struct Initializers<4>
224 {
225 static constexpr std::array<unsigned int, 16>
226 ucd_to_deal()
227 {
244 }
245
246 static constexpr std::array<unsigned int, 8>
247 unit_normal_direction()
248 {
249 return {{0, 0, 1, 1, 2, 2, 3, 3}};
250 }
251
252 static constexpr std::array<int, 8>
253 unit_normal_orientation()
254 {
255 return {{-1, 1, -1, 1, -1, 1, -1, 1}};
256 }
257
258 static constexpr std::array<Tensor<1, 4>, 8>
259 unit_normal_vector()
260 {
261 return {{Tensor<1, 4>{{-1, 0, 0, 0}},
262 Tensor<1, 4>{{1, 0, 0, 0}},
263 Tensor<1, 4>{{0, -1, 0, 0}},
264 Tensor<1, 4>{{0, 1, 0, 0}},
265 Tensor<1, 4>{{0, 0, -1, 0}},
266 Tensor<1, 4>{{0, 0, 1, 0}},
267 Tensor<1, 4>{{0, 0, 0, -1}},
268 Tensor<1, 4>{{0, 0, 0, 1}}}};
269 }
270
271 static constexpr ::ndarray<Tensor<1, 4>, 8, 3>
272 unit_tangential_vectors()
273 {
274 return {{{{Tensor<1, 4>{{0, -1, 0, 0}},
275 Tensor<1, 4>{{0, 0, 1, 0}},
276 Tensor<1, 4>{{0, 0, 0, 1}}}},
277 {{Tensor<1, 4>{{0, 1, 0, 0}},
278 Tensor<1, 4>{{0, 0, 1, 0}},
279 Tensor<1, 4>{{0, 0, 0, 1}}}},
280 {{Tensor<1, 4>{{0, 0, -1, 0}},
281 Tensor<1, 4>{{0, 0, 0, 1}},
282 Tensor<1, 4>{{1, 0, 0, 0}}}},
283 {{Tensor<1, 4>{{0, 0, 1, 0}},
284 Tensor<1, 4>{{0, 0, 0, 1}},
285 Tensor<1, 4>{{1, 0, 0, 0}}}},
286 {{Tensor<1, 4>{{0, 0, 0, -1}},
287 Tensor<1, 4>{{1, 0, 0, 0}},
288 Tensor<1, 4>{{0, 1, 0, 0}}}},
289 {{Tensor<1, 4>{{0, 0, 0, 1}},
290 Tensor<1, 4>{{1, 0, 0, 0}},
291 Tensor<1, 4>{{0, 1, 0, 0}}}},
292 {{Tensor<1, 4>{{-1, 0, 0, 0}},
293 Tensor<1, 4>{{0, 1, 0, 0}},
294 Tensor<1, 4>{{0, 0, 1, 0}}}},
295 {{Tensor<1, 4>{{1, 0, 0, 0}},
296 Tensor<1, 4>{{0, 1, 0, 0}},
297 Tensor<1, 4>{{0, 0, 1, 0}}}}}};
298 }
299
300 static constexpr std::array<unsigned int, 8>
301 opposite_face()
302 {
303 return {{1, 0, 3, 2, 5, 4, 7, 6}};
304 }
305
306 static constexpr std::array<unsigned int, 16>
307 dx_to_deal()
308 {
325 }
326
327 static constexpr ::ndarray<unsigned int, 16, 4>
328 vertex_to_face()
329 {
394 }
395 };
396 } // namespace GeometryInfoHelper
397} // namespace internal
398#endif // DOXYGEN
399
400
417{
418public:
426 {
434 line = 1,
438 quad = 2,
442 hex = 3
443 };
444
450
456 GeometryPrimitive(const unsigned int object_dimension);
457
462 operator unsigned int() const;
463
464private:
469};
470
471
472
485template <int dim>
487{
522 enum Possibilities : std::uint8_t
523 {
528
541 };
542};
543
544
545
555template <>
557{
592 enum Possibilities : std::uint8_t
593 {
601 cut_x = 1,
606 };
607};
608
609
610
621template <>
623{
658 enum Possibilities : std::uint8_t
659 {
667 cut_x = 1,
671 cut_y = 2,
675 cut_xy = cut_x | cut_y,
676
680 isotropic_refinement = cut_xy
681 };
682};
683
684
685
696template <>
698{
733 enum Possibilities : std::uint8_t
734 {
742 cut_x = 1,
746 cut_y = 2,
750 cut_xy = cut_x | cut_y,
754 cut_z = 4,
758 cut_xz = cut_x | cut_z,
762 cut_yz = cut_y | cut_z,
766 cut_xyz = cut_x | cut_y | cut_z,
767
771 isotropic_refinement = cut_xyz
772 };
773};
774
775
776
787template <int dim>
789{
790public:
794 static constexpr unsigned int n_refinement_cases = (1 << dim);
795
800
806 const typename RefinementPossibilities<dim>::Possibilities refinement_case);
807
813 explicit RefinementCase(const std::uint8_t refinement_case);
814
826 DEAL_II_HOST_DEVICE operator std::uint8_t() const;
827
833 operator|(const RefinementCase &r) const;
834
840 operator&(const RefinementCase &r) const;
841
850 operator~() const;
851
852
858 static RefinementCase
859 cut_axis(const unsigned int i);
860
869 static std::array<RefinementCase<dim>, n_refinement_cases>
871
875 static std::size_t
877
883 template <class Archive>
884 void
885 serialize(Archive &ar, const unsigned int version);
886
892 int,
893 << "The refinement flags given (" << arg1
894 << ") contain set bits that do not "
895 << "make sense for the space dimension of the object to which they are applied.");
896
897private:
902 std::uint8_t value : (dim > 0 ? dim : 1);
903};
904
905
906
907namespace internal
908{
928 template <int dim>
930 {
935 {
940
944 case_isotropic = static_cast<std::uint8_t>(-1)
945 };
946 };
947
948
957 template <>
959 {
977 };
978
979
980
990 template <>
992 {
1010 };
1011
1012
1013
1024 template <>
1026 {
1034 {
1042 case_x = 1,
1046 case_isotropic = case_x
1048 };
1049
1050
1051
1143 template <>
1145 {
1153 {
1155 case_x = 1,
1156 case_x1y = 2,
1157 case_x2y = 3,
1158 case_x1y2y = 4,
1159 case_y = 5,
1160 case_y1x = 6,
1161 case_y2x = 7,
1162 case_y1x2x = 8,
1163 case_xy = 9,
1164
1165 case_isotropic = case_xy
1167 };
1168
1169
1170
1177 template <int dim>
1179 {
1180 public:
1187 subface_possibility);
1188
1200 operator std::uint8_t() const;
1201
1205 static constexpr std::size_t
1207
1213 int,
1214 << "The subface case given (" << arg1 << ") does not make sense "
1215 << "for the space dimension of the object to which they are applied.");
1216
1217 private:
1222 std::uint8_t value : (dim == 3 ? 4 : 1);
1223 };
1224
1225} // namespace internal
1226
1227
1228
1229template <int dim>
1230struct GeometryInfo;
1231
1232
1233
1253template <>
1255{
1263 static constexpr unsigned int max_children_per_cell = 1;
1264
1268 static constexpr unsigned int faces_per_cell = 0;
1269
1286 static std::array<unsigned int, 0>
1288
1296 static constexpr unsigned int max_children_per_face = 0;
1297
1303 static unsigned int
1304 n_children(const RefinementCase<0> &refinement_case);
1305
1309 static constexpr unsigned int vertices_per_cell = 1;
1310
1329 static std::array<unsigned int, vertices_per_cell>
1331
1355 static unsigned int
1356 face_to_cell_vertices(const unsigned int face,
1357 const unsigned int vertex,
1358 const bool face_orientation = true,
1359 const bool face_flip = false,
1360 const bool face_rotation = false);
1361
1376 static unsigned int
1377 face_to_cell_lines(const unsigned int face,
1378 const unsigned int line,
1379 const bool face_orientation = true,
1380 const bool face_flip = false,
1381 const bool face_rotation = false);
1382
1389 static constexpr unsigned int vertices_per_face = 0;
1390
1394 static constexpr unsigned int lines_per_face = 0;
1395
1399 static constexpr unsigned int quads_per_face = 0;
1400
1404 static constexpr unsigned int lines_per_cell = 0;
1405
1409 static constexpr unsigned int quads_per_cell = 0;
1410
1414 static constexpr unsigned int hexes_per_cell = 0;
1415
1433 static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1434
1448 static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1449};
1450
1451
1452
1984template <int dim>
1986{
1994 static constexpr unsigned int max_children_per_cell = 1 << dim;
1995
1999 static constexpr unsigned int faces_per_cell = 2 * dim;
2000
2019
2027 static constexpr unsigned int max_children_per_face =
2029
2033 static constexpr unsigned int vertices_per_cell = 1 << dim;
2034
2052
2056 static constexpr unsigned int vertices_per_face =
2058
2062 static constexpr unsigned int lines_per_face =
2064
2068 static constexpr unsigned int quads_per_face =
2070
2080 static constexpr unsigned int lines_per_cell =
2081 (2 * GeometryInfo<dim - 1>::lines_per_cell +
2083
2091 static constexpr unsigned int quads_per_cell =
2092 (2 * GeometryInfo<dim - 1>::quads_per_cell +
2094
2098 static constexpr unsigned int hexes_per_cell =
2099 (2 * GeometryInfo<dim - 1>::hexes_per_cell +
2101
2119 static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
2120 internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal();
2121
2135 static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
2136 internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal();
2137
2150 internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face();
2151
2156 static unsigned int
2157 n_children(const RefinementCase<dim> &refinement_case);
2158
2163 static unsigned int
2165
2175 static double
2177 const unsigned int subface_no);
2178
2184 static RefinementCase<dim - 1>
2185 face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2186 const unsigned int face_no,
2187 const bool face_orientation = true,
2188 const bool face_flip = false,
2189 const bool face_rotation = false);
2190
2196 static RefinementCase<dim>
2199 const unsigned int face_no,
2200 const bool face_orientation = true,
2201 const bool face_flip = false,
2202 const bool face_rotation = false);
2203
2208 static RefinementCase<1>
2209 line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2210 const unsigned int line_no);
2211
2216 static RefinementCase<dim>
2218
2265 static unsigned int
2267 const unsigned int face,
2268 const unsigned int subface,
2269 const bool face_orientation = true,
2270 const bool face_flip = false,
2271 const bool face_rotation = false,
2274
2288 static unsigned int
2289 line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
2290
2311 static unsigned int
2312 face_to_cell_vertices(const unsigned int face,
2313 const unsigned int vertex,
2314 const bool face_orientation = true,
2315 const bool face_flip = false,
2316 const bool face_rotation = false);
2317
2329 static unsigned int
2330 face_to_cell_lines(const unsigned int face,
2331 const unsigned int line,
2332 const bool face_orientation = true,
2333 const bool face_flip = false,
2334 const bool face_rotation = false);
2335
2345 static unsigned int
2346 standard_to_real_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 real_to_standard_face_vertex(const unsigned int vertex,
2362 const bool face_orientation = true,
2363 const bool face_flip = false,
2364 const bool face_rotation = false);
2365
2375 static unsigned int
2376 standard_to_real_face_line(const unsigned int line,
2377 const bool face_orientation = true,
2378 const bool face_flip = false,
2379 const bool face_rotation = false);
2380
2386 static unsigned int
2387 standard_to_real_line_vertex(const unsigned int vertex,
2388 const bool line_orientation = true);
2389
2397 static std::array<unsigned int, 2>
2399
2407 static std::array<unsigned int, 2>
2409
2417 static std::array<unsigned int, 2>
2419
2429 static unsigned int
2430 real_to_standard_face_line(const unsigned int line,
2431 const bool face_orientation = true,
2432 const bool face_flip = false,
2433 const bool face_rotation = false);
2434
2440 static Point<dim>
2441 unit_cell_vertex(const unsigned int vertex);
2442
2452 static unsigned int
2454
2462 static Point<dim>
2464 const unsigned int child_index,
2465 const RefinementCase<dim> refine_case =
2467
2473 static Point<dim>
2475 const unsigned int child_index,
2476 const RefinementCase<dim> refine_case =
2478
2483 static bool
2485
2497 static bool
2498 is_inside_unit_cell(const Point<dim> &p, const double eps);
2499
2504 template <typename Number = double>
2505 static Point<dim, Number>
2507
2513 static double
2515
2520 static double
2521 d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2522
2527 static Tensor<1, dim>
2528 d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2529
2581 template <int spacedim>
2582 static void
2584#ifndef DEAL_II_CXX14_CONSTEXPR_BUG
2586 Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2587#else
2588 (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2589#endif
2590
2600 static constexpr std::array<unsigned int, faces_per_cell>
2602 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction();
2603
2620 static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2621 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation();
2622
2633 static constexpr std::array<Tensor<1, dim>, faces_per_cell>
2635 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_vector();
2636
2650 static constexpr ndarray<Tensor<1, dim>, faces_per_cell, dim - 1>
2651 unit_tangential_vectors = internal::GeometryInfoHelper::Initializers<
2653
2659 static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2660 internal::GeometryInfoHelper::Initializers<dim>::opposite_face();
2661
2662
2667 double,
2668 << "The coordinates must satisfy 0 <= x_i <= 1, "
2669 << "but here we have x_i=" << arg1);
2670
2675 int,
2676 int,
2677 int,
2678 << "RefinementCase<dim> " << arg1 << ": face " << arg2
2679 << " has no subface " << arg3);
2680};
2681
2682
2683
2684#ifndef DOXYGEN
2685
2686
2687/* -------------- declaration of explicit specializations ------------- */
2688
2689
2690template <>
2693 const unsigned int i);
2694template <>
2697 const unsigned int i);
2698template <>
2701 const unsigned int i);
2702
2703
2704
2705/* -------------- inline functions ------------- */
2706
2707
2708inline GeometryPrimitive::GeometryPrimitive(const Object object)
2709 : object(object)
2710{}
2711
2712
2713
2714inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2715 : object(static_cast<Object>(object_dimension))
2716{}
2717
2718
2719inline GeometryPrimitive::operator unsigned int() const
2720{
2721 return static_cast<unsigned int>(object);
2722}
2723
2724
2725
2726namespace internal
2727{
2728 template <int dim>
2730 const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2731 : value(subface_possibility)
2732 {}
2733
2734
2735 template <int dim>
2736 inline SubfaceCase<dim>::operator std::uint8_t() const
2737 {
2738 return value;
2739 }
2740
2741
2742} // namespace internal
2743
2744
2745template <int dim>
2747RefinementCase<dim>::cut_axis(const unsigned int)
2748{
2750 return static_cast<std::uint8_t>(-1);
2751}
2752
2753
2754template <>
2755inline RefinementCase<1>
2756RefinementCase<1>::cut_axis(const unsigned int i)
2757{
2758 AssertIndexRange(i, 1);
2759
2761 return options[i];
2762}
2763
2764
2765
2766template <>
2767inline RefinementCase<2>
2768RefinementCase<2>::cut_axis(const unsigned int i)
2769{
2770 AssertIndexRange(i, 2);
2771
2774 return options[i];
2775}
2776
2777
2778
2779template <>
2780inline RefinementCase<3>
2781RefinementCase<3>::cut_axis(const unsigned int i)
2782{
2783 AssertIndexRange(i, 3);
2784
2788 return options[i];
2789}
2790
2791
2792
2793template <>
2794inline std::array<RefinementCase<1>, 2>
2796{
2799}
2800
2801
2802
2803template <>
2804inline std::array<RefinementCase<2>, 4>
2806{
2811}
2812
2813
2814
2815template <>
2816inline std::array<RefinementCase<3>, 8>
2818{
2827}
2828
2829
2830
2831template <int dim>
2833 : value(RefinementPossibilities<dim>::no_refinement)
2834{}
2835
2836
2837
2838template <int dim>
2840 const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2841 : value(refinement_case)
2842{
2843 // check that only those bits of
2844 // the given argument are set that
2845 // make sense for a given space
2846 // dimension
2847 Assert((refinement_case &
2849 refinement_case,
2850 ExcInvalidRefinementCase(refinement_case));
2851}
2852
2853
2854
2855template <int dim>
2856inline RefinementCase<dim>::RefinementCase(const std::uint8_t refinement_case)
2857 : value(refinement_case)
2858{
2859 // check that only those bits of
2860 // the given argument are set that
2861 // make sense for a given space
2862 // dimension
2863 Assert((refinement_case &
2865 refinement_case,
2866 ExcInvalidRefinementCase(refinement_case));
2867}
2868
2869
2870
2871template <int dim>
2872inline DEAL_II_HOST_DEVICE RefinementCase<dim>::operator std::uint8_t() const
2873{
2874 return value;
2875}
2876
2877
2878
2879template <int dim>
2882{
2883 return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2884}
2885
2886
2887
2888template <int dim>
2891{
2892 return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2893}
2894
2895
2896
2897template <int dim>
2900{
2901 return RefinementCase<dim>(static_cast<std::uint8_t>(
2903}
2904
2905
2906
2907template <int dim>
2908inline std::size_t
2910{
2911 return sizeof(RefinementCase<dim>);
2912}
2913
2914
2915
2916template <int dim>
2917template <class Archive>
2918inline void
2919RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2920{
2921 // serialization can't deal with bitfields, so copy from/to a full sized
2922 // std::uint8_t
2923 std::uint8_t uchar_value = value;
2924 ar &uchar_value;
2925 value = uchar_value;
2926}
2927
2928
2929
2930template <>
2931inline Point<1>
2932GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2933{
2934 AssertIndexRange(vertex, vertices_per_cell);
2935
2936 return Point<1>(static_cast<double>(vertex));
2937}
2938
2939
2940
2941template <>
2942inline Point<2>
2943GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2944{
2945 AssertIndexRange(vertex, vertices_per_cell);
2946
2947 return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2948}
2949
2950
2951
2952template <>
2953inline Point<3>
2954GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2955{
2956 AssertIndexRange(vertex, vertices_per_cell);
2957
2958 return {static_cast<double>(vertex % 2),
2959 static_cast<double>(vertex / 2 % 2),
2960 static_cast<double>(vertex / 4)};
2961}
2962
2963
2964
2965inline std::array<unsigned int, 0>
2967{
2968 return {{}};
2969}
2970
2971
2972
2973inline std::array<unsigned int, 1>
2975{
2976 return {{0}};
2977}
2978
2979
2980
2981template <int dim>
2984{
2985 return {0U, faces_per_cell};
2986}
2987
2988
2989
2990template <int dim>
2993{
2994 return {0U, vertices_per_cell};
2995}
2996
2997
2998
2999template <int dim>
3000inline Point<dim>
3001GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
3002{
3004
3005 return {};
3006}
3007
3008
3009
3010template <>
3011inline unsigned int
3013{
3014 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
3015
3016 return (p[0] <= 0.5 ? 0 : 1);
3017}
3018
3019
3020
3021template <>
3022inline unsigned int
3024{
3025 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
3026 Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
3027
3028 return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
3029}
3030
3031
3032
3033template <>
3034inline unsigned int
3036{
3037 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
3038 Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
3039 Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
3040
3041 return (p[0] <= 0.5 ?
3042 (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
3043 (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
3044}
3045
3046
3047template <int dim>
3048inline unsigned int
3050{
3052
3053 return 0;
3054}
3055
3056
3057
3058template <>
3059inline Point<1>
3061 const unsigned int child_index,
3062 const RefinementCase<1> refine_case)
3063
3064{
3065 AssertIndexRange(child_index, 2);
3067 (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3068
3069 return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
3070}
3071
3072
3073
3074template <>
3075inline Point<2>
3077 const unsigned int child_index,
3078 const RefinementCase<2> refine_case)
3079
3080{
3081 AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3082
3083 Point<2> point = p;
3084 switch (refine_case)
3085 {
3087 point[0] *= 2.0;
3088 if (child_index == 1)
3089 point[0] -= 1.0;
3090 break;
3092 point[1] *= 2.0;
3093 if (child_index == 1)
3094 point[1] -= 1.0;
3095 break;
3097 point *= 2.0;
3098 point -= unit_cell_vertex(child_index);
3099 break;
3100 default:
3102 }
3103
3104 return point;
3105}
3106
3107
3108
3109template <>
3110inline Point<3>
3112 const unsigned int child_index,
3113 const RefinementCase<3> refine_case)
3114
3115{
3116 AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3117
3118 Point<3> point = p;
3119 // there might be a cleverer way to do
3120 // this, but since this function is called
3121 // in very few places for initialization
3122 // purposes only, I don't care at the
3123 // moment
3124 switch (refine_case)
3125 {
3127 point[0] *= 2.0;
3128 if (child_index == 1)
3129 point[0] -= 1.0;
3130 break;
3132 point[1] *= 2.0;
3133 if (child_index == 1)
3134 point[1] -= 1.0;
3135 break;
3137 point[2] *= 2.0;
3138 if (child_index == 1)
3139 point[2] -= 1.0;
3140 break;
3142 point[0] *= 2.0;
3143 point[1] *= 2.0;
3144 if (child_index % 2 == 1)
3145 point[0] -= 1.0;
3146 if (child_index / 2 == 1)
3147 point[1] -= 1.0;
3148 break;
3150 // careful, this is slightly
3151 // different from xy and yz due to
3152 // different internal numbering of
3153 // children!
3154 point[0] *= 2.0;
3155 point[2] *= 2.0;
3156 if (child_index / 2 == 1)
3157 point[0] -= 1.0;
3158 if (child_index % 2 == 1)
3159 point[2] -= 1.0;
3160 break;
3162 point[1] *= 2.0;
3163 point[2] *= 2.0;
3164 if (child_index % 2 == 1)
3165 point[1] -= 1.0;
3166 if (child_index / 2 == 1)
3167 point[2] -= 1.0;
3168 break;
3170 point *= 2.0;
3171 point -= unit_cell_vertex(child_index);
3172 break;
3173 default:
3175 }
3176
3177 return point;
3178}
3179
3180
3181
3182template <int dim>
3183inline Point<dim>
3185 const Point<dim> & /*p*/,
3186 const unsigned int /*child_index*/,
3187 const RefinementCase<dim> /*refine_case*/)
3188
3189{
3191 return {};
3192}
3193
3194
3195
3196template <>
3197inline Point<1>
3199 const unsigned int child_index,
3200 const RefinementCase<1> refine_case)
3201
3202{
3203 AssertIndexRange(child_index, 2);
3205 (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3206
3207 return (p + unit_cell_vertex(child_index)) * 0.5;
3208}
3209
3210
3211
3212template <>
3213inline Point<3>
3215 const unsigned int child_index,
3216 const RefinementCase<3> refine_case)
3217
3218{
3219 AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3220
3221 Point<3> point = p;
3222 // there might be a cleverer way to do
3223 // this, but since this function is called
3224 // in very few places for initialization
3225 // purposes only, I don't care at the
3226 // moment
3227 switch (refine_case)
3228 {
3230 if (child_index == 1)
3231 point[0] += 1.0;
3232 point[0] *= 0.5;
3233 break;
3235 if (child_index == 1)
3236 point[1] += 1.0;
3237 point[1] *= 0.5;
3238 break;
3240 if (child_index == 1)
3241 point[2] += 1.0;
3242 point[2] *= 0.5;
3243 break;
3245 if (child_index % 2 == 1)
3246 point[0] += 1.0;
3247 if (child_index / 2 == 1)
3248 point[1] += 1.0;
3249 point[0] *= 0.5;
3250 point[1] *= 0.5;
3251 break;
3253 // careful, this is slightly
3254 // different from xy and yz due to
3255 // different internal numbering of
3256 // children!
3257 if (child_index / 2 == 1)
3258 point[0] += 1.0;
3259 if (child_index % 2 == 1)
3260 point[2] += 1.0;
3261 point[0] *= 0.5;
3262 point[2] *= 0.5;
3263 break;
3265 if (child_index % 2 == 1)
3266 point[1] += 1.0;
3267 if (child_index / 2 == 1)
3268 point[2] += 1.0;
3269 point[1] *= 0.5;
3270 point[2] *= 0.5;
3271 break;
3273 point += unit_cell_vertex(child_index);
3274 point *= 0.5;
3275 break;
3276 default:
3278 }
3279
3280 return point;
3281}
3282
3283
3284
3285template <>
3286inline Point<2>
3288 const unsigned int child_index,
3289 const RefinementCase<2> refine_case)
3290{
3291 AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3292
3293 Point<2> point = p;
3294 switch (refine_case)
3295 {
3297 if (child_index == 1)
3298 point[0] += 1.0;
3299 point[0] *= 0.5;
3300 break;
3302 if (child_index == 1)
3303 point[1] += 1.0;
3304 point[1] *= 0.5;
3305 break;
3307 point += unit_cell_vertex(child_index);
3308 point *= 0.5;
3309 break;
3310 default:
3312 }
3313
3314 return point;
3315}
3316
3317
3318
3319template <int dim>
3320inline Point<dim>
3322 const Point<dim> & /*p*/,
3323 const unsigned int /*child_index*/,
3324 const RefinementCase<dim> /*refine_case*/)
3325{
3327 return {};
3328}
3329
3330
3331
3332template <int dim>
3333inline bool
3335{
3337 return false;
3338}
3339
3340template <>
3341inline bool
3343{
3344 return (p[0] >= 0.) && (p[0] <= 1.);
3345}
3346
3347
3348
3349template <>
3350inline bool
3352{
3353 return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
3354}
3355
3356
3357
3358template <>
3359inline bool
3361{
3362 return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
3363 (p[2] >= 0.) && (p[2] <= 1.);
3364}
3365
3366
3367
3368template <int dim>
3369inline bool
3371{
3373 return false;
3374}
3375
3376template <>
3377inline bool
3378GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
3379{
3380 return (p[0] >= -eps) && (p[0] <= 1. + eps);
3381}
3382
3383
3384
3385template <>
3386inline bool
3387GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
3388{
3389 const double l = -eps, u = 1 + eps;
3390 return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
3391}
3392
3393
3394
3395template <>
3396inline bool
3397GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
3398{
3399 const double l = -eps, u = 1.0 + eps;
3400 return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
3401 (p[2] >= l) && (p[2] <= u);
3402}
3403
3404
3405
3406template <>
3407inline unsigned int
3408GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
3409 const unsigned int vertex)
3410{
3411 (void)line;
3412 AssertIndexRange(line, lines_per_cell);
3413 AssertIndexRange(vertex, 2);
3414
3415 return vertex;
3416}
3417
3418
3419template <>
3420inline unsigned int
3421GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
3422 const unsigned int vertex)
3423{
3424 constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
3425 return cell_vertices[line][vertex];
3426}
3427
3428
3429
3430template <>
3431inline unsigned int
3432GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
3433 const unsigned int vertex)
3434{
3435 AssertIndexRange(line, lines_per_cell);
3436 AssertIndexRange(vertex, 2);
3437
3438 constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3439 {1, 3},
3440 {0, 1},
3441 {2, 3},
3442 {4, 6}, // top face
3443 {5, 7},
3444 {4, 5},
3445 {6, 7},
3446 {0,
3447 4}, // connects of bottom
3448 {1, 5}, // top face
3449 {2, 6},
3450 {3, 7}};
3451 return vertices[line][vertex];
3452}
3453
3454
3455template <>
3456inline unsigned int
3457GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3458{
3461}
3462
3463template <>
3464inline unsigned int
3465GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3466 const bool face_orientation,
3467 const bool face_flip,
3468 const bool face_rotation)
3469{
3471
3472 // set up a table to make sure that
3473 // we handle non-standard faces correctly
3474 //
3475 // so set up a table that for each vertex (of
3476 // a quad in standard position) describes
3477 // which vertex to take
3478 //
3479 // first index: four vertices 0...3
3480 //
3481 // second index: face_orientation; 0:
3482 // opposite normal, 1: standard
3483 //
3484 // third index: face_flip; 0: standard, 1:
3485 // face rotated by 180 degrees
3486 //
3487 // forth index: face_rotation: 0: standard,
3488 // 1: face rotated by 90 degrees
3489
3490 constexpr unsigned int vertex_translation[4][2][2][2] = {
3491 {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3492 // face_rotation=false and true
3493 {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3494 // face_rotation=false and true
3495 {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3496 // face_rotation=false and true
3497 {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3498 // face_rotation=false and true
3499
3500 {{{2, 3}, // vertex 1 ...
3501 {1, 0}},
3502 {{1, 0}, {2, 3}}},
3503
3504 {{{1, 0}, // vertex 2 ...
3505 {2, 3}},
3506 {{2, 3}, {1, 0}}},
3507
3508 {{{3, 1}, // vertex 3 ...
3509 {0, 2}},
3510 {{3, 1}, {0, 2}}}};
3511
3512 return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3513}
3514
3515
3516
3517template <int dim>
3518inline unsigned int
3519GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3520 const bool,
3521 const bool,
3522 const bool)
3523{
3524 Assert(dim > 1, ExcImpossibleInDim(dim));
3526 return vertex;
3527}
3528
3529template <int dim>
3530inline unsigned int
3532{
3533 constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3534 0, 2, 2, 4, 2, 4, 4, 8};
3535
3536 return n_children[ref_case];
3537}
3538
3539
3540
3541template <int dim>
3542inline unsigned int
3544{
3546 return 0;
3547}
3548
3549template <>
3550inline unsigned int
3552{
3553 Assert(false, ExcImpossibleInDim(1));
3554 return 0;
3555}
3556
3557template <>
3558inline unsigned int
3560{
3561 return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3562}
3563
3564
3565
3566template <>
3567inline unsigned int
3569{
3570 const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3571 0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3572 return nsubs[subface_case];
3573}
3574
3575
3576
3577template <int dim>
3578inline double
3580 const unsigned int)
3581{
3583 return 0.;
3584}
3585
3586template <>
3587inline double
3589 const unsigned int)
3590{
3591 return 1;
3592}
3593
3594
3595template <>
3596inline double
3598 const unsigned int)
3599{
3600 double ratio = 1;
3601 switch (subface_case)
3602 {
3604 // Here, an
3605 // Assert(false,ExcInternalError())
3606 // would be the right
3607 // choice, but
3608 // unfortunately the
3609 // current function is
3610 // also called for faces
3611 // without children (see
3612 // tests/fe/mapping.cc).
3613 // Assert(false, ExcMessage("Face has no subfaces."));
3614 // Furthermore, assign
3615 // following value as
3616 // otherwise the
3617 // bits/volume_x tests
3618 // break
3620 break;
3622 ratio = 0.5;
3623 break;
3624 default:
3625 // there should be no
3626 // cases left
3628 break;
3629 }
3630
3631 return ratio;
3632}
3633
3634
3635template <>
3636inline double
3638 const unsigned int subface_no)
3639{
3640 double ratio = 1;
3641 switch (subface_case)
3642 {
3644 // Here, an
3645 // Assert(false,ExcInternalError())
3646 // would be the right
3647 // choice, but
3648 // unfortunately the
3649 // current function is
3650 // also called for faces
3651 // without children (see
3652 // tests/bits/mesh_3d_16.cc). Add
3653 // following switch to
3654 // avoid diffs in
3655 // tests/bits/mesh_3d_16
3657 break;
3660 ratio = 0.5;
3661 break;
3665 ratio = 0.25;
3666 break;
3669 if (subface_no < 2)
3670 ratio = 0.25;
3671 else
3672 ratio = 0.5;
3673 break;
3676 if (subface_no == 0)
3677 ratio = 0.5;
3678 else
3679 ratio = 0.25;
3680 break;
3681 default:
3682 // there should be no
3683 // cases left
3685 break;
3686 }
3687
3688 return ratio;
3689}
3690
3691
3692
3693template <int dim>
3695 const RefinementCase<dim> &,
3696 const unsigned int,
3697 const bool,
3698 const bool,
3699 const bool)
3700{
3702 return RefinementCase<dim - 1>::no_refinement;
3703}
3704
3705template <>
3707 const RefinementCase<1> &,
3708 const unsigned int,
3709 const bool,
3710 const bool,
3711 const bool)
3712{
3713 Assert(false, ExcImpossibleInDim(1));
3714
3716}
3717
3718
3719template <>
3720inline RefinementCase<1>
3722 const RefinementCase<2> &cell_refinement_case,
3723 const unsigned int face_no,
3724 const bool,
3725 const bool,
3726 const bool)
3727{
3728 const unsigned int dim = 2;
3729 AssertIndexRange(cell_refinement_case,
3732
3733 // simple special case
3734 if (cell_refinement_case == RefinementCase<dim>::cut_xy)
3736
3737 const RefinementCase<dim - 1>
3740 {RefinementCase<dim - 1>::no_refinement, // no_refinement
3741 RefinementCase<dim - 1>::no_refinement},
3742
3743 {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3744
3745 {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3746
3747 {RefinementCase<dim - 1>::cut_x, // cut_xy
3748 RefinementCase<dim - 1>::cut_x}};
3749
3750 return ref_cases[cell_refinement_case][face_no / 2];
3751}
3752
3753
3754template <>
3755inline RefinementCase<2>
3757 const RefinementCase<3> &cell_refinement_case,
3758 const unsigned int face_no,
3759 const bool face_orientation,
3760 const bool /*face_flip*/,
3761 const bool face_rotation)
3762{
3763 const unsigned int dim = 3;
3764 AssertIndexRange(cell_refinement_case,
3767
3768 // simple special case
3769 if (cell_refinement_case == RefinementCase<dim>::cut_xyz)
3770 return RefinementCase<dim - 1>::cut_xy;
3771
3772 const RefinementCase<dim - 1>
3775 {RefinementCase<dim - 1>::no_refinement, // no_refinement
3776 RefinementCase<dim - 1>::no_refinement,
3777 RefinementCase<dim - 1>::no_refinement},
3778
3779 {RefinementCase<dim - 1>::no_refinement, // cut_x
3780 RefinementCase<dim - 1>::cut_y,
3781 RefinementCase<dim - 1>::cut_x},
3782
3783 {RefinementCase<dim - 1>::cut_x, // cut_y
3784 RefinementCase<dim - 1>::no_refinement,
3785 RefinementCase<dim - 1>::cut_y},
3786
3787 {RefinementCase<dim - 1>::cut_x, // cut_xy
3788 RefinementCase<dim - 1>::cut_y,
3789 RefinementCase<dim - 1>::cut_xy},
3790
3791 {RefinementCase<dim - 1>::cut_y, // cut_z
3792 RefinementCase<dim - 1>::cut_x,
3793 RefinementCase<dim - 1>::no_refinement},
3794
3795 {RefinementCase<dim - 1>::cut_y, // cut_xz
3796 RefinementCase<dim - 1>::cut_xy,
3797 RefinementCase<dim - 1>::cut_x},
3798
3799 {RefinementCase<dim - 1>::cut_xy, // cut_yz
3800 RefinementCase<dim - 1>::cut_x,
3801 RefinementCase<dim - 1>::cut_y},
3802
3803 {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3804 RefinementCase<dim - 1>::cut_xy,
3805 RefinementCase<dim - 1>::cut_xy},
3806 };
3807
3808 const RefinementCase<dim - 1> ref_case =
3809 ref_cases[cell_refinement_case][face_no / 2];
3810
3811 const RefinementCase<dim - 1> flip[4] = {
3812 RefinementCase<dim - 1>::no_refinement,
3813 RefinementCase<dim - 1>::cut_y,
3814 RefinementCase<dim - 1>::cut_x,
3815 RefinementCase<dim - 1>::cut_xy};
3816
3817 // correct the ref_case for face_orientation
3818 // and face_rotation. for face_orientation,
3819 // 'true' is the default value whereas for
3820 // face_rotation, 'false' is standard. If
3821 // <tt>face_rotation==face_orientation</tt>,
3822 // then one of them is non-standard and we
3823 // have to swap cut_x and cut_y, otherwise no
3824 // change is necessary. face_flip has no
3825 // influence. however, in order to keep the
3826 // interface consistent with other functions,
3827 // we still include it as an argument to this
3828 // function
3829 return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3830}
3831
3832
3833
3834template <int dim>
3835inline RefinementCase<1>
3837 const unsigned int)
3838{
3841}
3842
3843template <>
3844inline RefinementCase<1>
3846 const RefinementCase<1> &cell_refinement_case,
3847 const unsigned int line_no)
3848{
3849 (void)line_no;
3850 const unsigned int dim = 1;
3851 (void)dim;
3852 AssertIndexRange(cell_refinement_case,
3855
3856 return cell_refinement_case;
3857}
3858
3859
3860template <>
3861inline RefinementCase<1>
3863 const RefinementCase<2> &cell_refinement_case,
3864 const unsigned int line_no)
3865{
3866 // Assertions are in face_refinement_case()
3867 return face_refinement_case(cell_refinement_case, line_no);
3868}
3869
3870
3871template <>
3872inline RefinementCase<1>
3874 const RefinementCase<3> &cell_refinement_case,
3875 const unsigned int line_no)
3876{
3877 const unsigned int dim = 3;
3878 AssertIndexRange(cell_refinement_case,
3881
3882 // simple special case
3883 if (cell_refinement_case == RefinementCase<dim>::cut_xyz)
3885
3886 // array indicating, which simple refine
3887 // case cuts a line in direction x, y or
3888 // z. For example, cut_y and everything
3889 // containing cut_y (cut_xy, cut_yz,
3890 // cut_xyz) cuts lines, which are in y
3891 // direction.
3892 const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3895
3896 // order the direction of lines
3897 // 0->x, 1->y, 2->z
3898 const unsigned int direction[lines_per_cell] = {
3899 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3900
3901 return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3903 RefinementCase<1>::no_refinement);
3904}
3905
3906
3907
3908template <int dim>
3912 const unsigned int,
3913 const bool,
3914 const bool,
3915 const bool)
3916{
3918
3920}
3921
3922template <>
3923inline RefinementCase<1>
3925 const RefinementCase<0> &,
3926 const unsigned int,
3927 const bool,
3928 const bool,
3929 const bool)
3930{
3931 const unsigned int dim = 1;
3932 Assert(false, ExcImpossibleInDim(dim));
3933
3935}
3936
3937
3938template <>
3939inline RefinementCase<2>
3941 const RefinementCase<1> &face_refinement_case,
3942 const unsigned int face_no,
3943 const bool,
3944 const bool,
3945 const bool)
3946{
3947 const unsigned int dim = 2;
3948 AssertIndexRange(face_refinement_case,
3951
3952 if (face_refinement_case == RefinementCase<dim>::cut_x)
3953 return (face_no / 2) != 0u ? RefinementCase<dim>::cut_x :
3955 else
3957}
3958
3959
3960template <>
3961inline RefinementCase<3>
3963 const RefinementCase<2> &face_refinement_case,
3964 const unsigned int face_no,
3965 const bool face_orientation,
3966 const bool /*face_flip*/,
3967 const bool face_rotation)
3968{
3969 const unsigned int dim = 3;
3970 AssertIndexRange(face_refinement_case,
3973
3978
3979 // correct the face_refinement_case for
3980 // face_orientation and face_rotation. for
3981 // face_orientation, 'true' is the default
3982 // value whereas for face_rotation, 'false'
3983 // is standard. If
3984 // <tt>face_rotation==face_orientation</tt>,
3985 // then one of them is non-standard and we
3986 // have to swap cut_x and cut_y, otherwise no
3987 // change is necessary. face_flip has no
3988 // influence. however, in order to keep the
3989 // interface consistent with other functions,
3990 // we still include it as an argument to this
3991 // function
3992 const RefinementCase<dim - 1> std_face_ref =
3993 (face_orientation == face_rotation) ? flip[face_refinement_case] :
3994 face_refinement_case;
3995
3996 const RefinementCase<dim> face_to_cell[3][4] = {
3997 {RefinementCase<dim>::no_refinement, // faces 0 and 1
3998 RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
4001
4002 {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
4003 // "exchanged on faces 2 and 3")
4007
4008 {RefinementCase<dim>::no_refinement, // faces 4 and 5
4012
4013 return face_to_cell[face_no / 2][std_face_ref];
4014}
4015
4016
4017
4018template <int dim>
4021 const unsigned int)
4022{
4024
4026}
4027
4028template <>
4029inline RefinementCase<1>
4031 const unsigned int line_no)
4032{
4033 (void)line_no;
4034 AssertIndexRange(line_no, 1);
4035
4037}
4038
4039
4040template <>
4041inline RefinementCase<2>
4043 const unsigned int line_no)
4044{
4045 const unsigned int dim = 2;
4046 (void)dim;
4048
4049 return (line_no / 2) != 0u ? RefinementCase<2>::cut_x :
4051}
4052
4053
4054template <>
4055inline RefinementCase<3>
4057 const unsigned int line_no)
4058{
4059 const unsigned int dim = 3;
4061
4062 const RefinementCase<dim> ref_cases[6] = {
4063 RefinementCase<dim>::cut_y, // lines 0 and 1
4064 RefinementCase<dim>::cut_x, // lines 2 and 3
4065 RefinementCase<dim>::cut_y, // lines 4 and 5
4066 RefinementCase<dim>::cut_x, // lines 6 and 7
4067 RefinementCase<dim>::cut_z, // lines 8 and 9
4068 RefinementCase<dim>::cut_z}; // lines 10 and 11
4069
4070 return ref_cases[line_no / 2];
4071}
4072
4073
4074
4075template <>
4076inline unsigned int
4077GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
4078 const bool face_orientation,
4079 const bool face_flip,
4080 const bool face_rotation)
4081{
4083
4084 // set up a table to make sure that
4085 // we handle non-standard faces correctly
4086 //
4087 // so set up a table that for each vertex (of
4088 // a quad in standard position) describes
4089 // which vertex to take
4090 //
4091 // first index: four vertices 0...3
4092 //
4093 // second index: face_orientation; 0:
4094 // opposite normal, 1: standard
4095 //
4096 // third index: face_flip; 0: standard, 1:
4097 // face rotated by 180 degrees
4098 //
4099 // forth index: face_rotation: 0: standard,
4100 // 1: face rotated by 90 degrees
4101
4102 const unsigned int vertex_translation[4][2][2][2] = {
4103 {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
4104 // face_rotation=false and true
4105 {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
4106 // face_rotation=false and true
4107 {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
4108 // face_rotation=false and true
4109 {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
4110 // face_rotation=false and true
4111
4112 {{{2, 3}, // vertex 1 ...
4113 {1, 0}},
4114 {{1, 3}, {2, 0}}},
4115
4116 {{{1, 0}, // vertex 2 ...
4117 {2, 3}},
4118 {{2, 0}, {1, 3}}},
4119
4120 {{{3, 1}, // vertex 3 ...
4121 {0, 2}},
4122 {{3, 2}, {0, 1}}}};
4123
4124 return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
4125}
4126
4127
4128
4129template <int dim>
4130inline unsigned int
4131GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
4132 const bool,
4133 const bool,
4134 const bool)
4135{
4136 Assert(dim > 1, ExcImpossibleInDim(dim));
4138 return vertex;
4139}
4140
4141
4142
4143template <>
4144inline unsigned int
4145GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
4146 const bool face_orientation,
4147 const bool face_flip,
4148 const bool face_rotation)
4149{
4151
4152
4153 // make sure we handle
4154 // non-standard faces correctly
4155 //
4156 // so set up a table that for each line (of a
4157 // quad) describes which line to take
4158 //
4159 // first index: four lines 0...3
4160 //
4161 // second index: face_orientation; 0:
4162 // opposite normal, 1: standard
4163 //
4164 // third index: face_flip; 0: standard, 1:
4165 // face rotated by 180 degrees
4166 //
4167 // forth index: face_rotation: 0: standard,
4168 // 1: face rotated by 90 degrees
4169
4170 const unsigned int line_translation[4][2][2][2] = {
4171 {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4172 // face_rotation=false and true
4173 {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4174 // face_rotation=false and true
4175 {{0, 3}, // line 0, face_orientation=true, face_flip=false,
4176 // face_rotation=false and true
4177 {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
4178 // face_rotation=false and true
4179
4180 {{{3, 1}, // line 1 ...
4181 {2, 0}},
4182 {{1, 2}, {0, 3}}},
4183
4184 {{{0, 3}, // line 2 ...
4185 {1, 2}},
4186 {{2, 0}, {3, 1}}},
4187
4188 {{{1, 2}, // line 3 ...
4189 {0, 3}},
4190 {{3, 1}, {2, 0}}}};
4191
4192 return line_translation[line][face_orientation][face_flip][face_rotation];
4193}
4194
4195
4196
4197template <int dim>
4198inline unsigned int
4199GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
4200 const bool,
4201 const bool,
4202 const bool)
4203{
4205 return line;
4206}
4207
4208
4209
4210template <>
4211inline unsigned int
4212GeometryInfo<2>::standard_to_real_line_vertex(const unsigned int vertex,
4213 const bool line_orientation)
4214{
4215 return line_orientation ? vertex : (1 - vertex);
4216}
4217
4218
4219
4220template <int dim>
4221inline unsigned int
4222GeometryInfo<dim>::standard_to_real_line_vertex(const unsigned int vertex,
4223 const bool)
4224{
4226 return vertex;
4227}
4228
4229
4230
4231template <>
4232inline std::array<unsigned int, 2>
4234 const unsigned int vertex)
4235{
4236 return {{vertex % 2, vertex / 2}};
4237}
4238
4239
4240
4241template <int dim>
4242inline std::array<unsigned int, 2>
4244 const unsigned int vertex)
4245{
4247 (void)vertex;
4248 return {{0, 0}};
4249}
4250
4251
4252
4253template <>
4254inline std::array<unsigned int, 2>
4256{
4257 // set up a table that for each
4258 // line describes a) from which
4259 // quad to take it, b) which line
4260 // therein it is if the face is
4261 // oriented correctly
4262 static const unsigned int lookup_table[GeometryInfo<3>::lines_per_cell][2] = {
4263 {4, 0}, // take first four lines from bottom face
4264 {4, 1},
4265 {4, 2},
4266 {4, 3},
4267
4268 {5, 0}, // second four lines from top face
4269 {5, 1},
4270 {5, 2},
4271 {5, 3},
4272
4273 {0, 0}, // the rest randomly
4274 {1, 0},
4275 {0, 1},
4276 {1, 1}};
4277
4278 return {{lookup_table[i][0], lookup_table[i][1]}};
4279}
4280
4281
4282
4283template <int dim>
4284inline std::array<unsigned int, 2>
4286{
4288 (void)line;
4289 return {{0, 0}};
4290}
4291
4292
4293
4294template <>
4295inline std::array<unsigned int, 2>
4297 const unsigned int vertex)
4298{
4299 // get the corner indices by asking either the bottom or the top face for its
4300 // vertices. handle non-standard faces by calling the vertex reordering
4301 // function from GeometryInfo
4302
4303 // bottom face (4) for first four vertices, top face (5) for the rest
4304 return {{4 + vertex / 4, vertex % 4}};
4305}
4306
4307
4308
4309template <int dim>
4310inline std::array<unsigned int, 2>
4312 const unsigned int vertex)
4313{
4315 (void)vertex;
4316 return {{0, 0}};
4317}
4318
4319
4320
4321template <>
4322inline unsigned int
4323GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
4324 const bool face_orientation,
4325 const bool face_flip,
4326 const bool face_rotation)
4327{
4329
4330
4331 // make sure we handle
4332 // non-standard faces correctly
4333 //
4334 // so set up a table that for each line (of a
4335 // quad) describes which line to take
4336 //
4337 // first index: four lines 0...3
4338 //
4339 // second index: face_orientation; 0:
4340 // opposite normal, 1: standard
4341 //
4342 // third index: face_flip; 0: standard, 1:
4343 // face rotated by 180 degrees
4344 //
4345 // forth index: face_rotation: 0: standard,
4346 // 1: face rotated by 90 degrees
4347
4348 const unsigned int line_translation[4][2][2][2] = {
4349 {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4350 // face_rotation=false and true
4351 {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4352 // face_rotation=false and true
4353 {{0, 2}, // line 0, face_orientation=true, face_flip=false,
4354 // face_rotation=false and true
4355 {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
4356 // face_rotation=false and true
4357
4358 {{{3, 1}, // line 1 ...
4359 {2, 0}},
4360 {{1, 3}, {0, 2}}},
4361
4362 {{{0, 3}, // line 2 ...
4363 {1, 2}},
4364 {{2, 1}, {3, 0}}},
4365
4366 {{{1, 2}, // line 3 ...
4367 {0, 3}},
4368 {{3, 0}, {2, 1}}}};
4369
4370 return line_translation[line][face_orientation][face_flip][face_rotation];
4371}
4372
4373
4374
4375template <int dim>
4376inline unsigned int
4377GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
4378 const bool,
4379 const bool,
4380 const bool)
4381{
4383 return line;
4384}
4385
4386
4387
4388template <>
4389inline unsigned int
4391 const unsigned int face,
4392 const unsigned int subface,
4393 const bool,
4394 const bool,
4395 const bool,
4396 const RefinementCase<0> &)
4397{
4398 (void)subface;
4399 AssertIndexRange(face, faces_per_cell);
4400 AssertIndexRange(subface, max_children_per_face);
4401
4402 return face;
4403}
4404
4405
4406
4407template <>
4408inline unsigned int
4410 const unsigned int face,
4411 const unsigned int subface,
4412 const bool face_orientation,
4413 const bool /*face_flip*/,
4414 const bool /*face_rotation*/,
4415 const RefinementCase<1> &)
4416{
4417 AssertIndexRange(face, faces_per_cell);
4418 AssertIndexRange(subface, max_children_per_face);
4419
4420 // always return the child adjacent to the specified
4421 // subface. if the face of a cell is not refined, don't
4422 // throw an assertion but deliver the child adjacent to
4423 // the face nevertheless, i.e. deliver the child of
4424 // this cell adjacent to the subface of a possibly
4425 // refined neighbor. this simplifies setting neighbor
4426 // information in execute_refinement.
4427 constexpr unsigned int
4428 subcells[/* possible face orientation */ 2]
4429 [/* number of different ways to refine a cell */ 4]
4430 [/* faces_per_cell */ 4][/* max_children_per_face */ 2] = {
4431 {
4432 // Flipped orientation (face_orientation = false)
4433 {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
4434 {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
4435 {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
4436 },
4437 {
4438 // Normal orientation (face_orientation = true)
4439 {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
4440 {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
4441 {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
4442 }};
4443
4444 return subcells[face_orientation][ref_case - 1][face][subface];
4445}
4446
4447
4448
4449template <>
4450inline unsigned int
4452 const unsigned int face,
4453 const unsigned int subface,
4454 const bool face_orientation,
4455 const bool face_flip,
4456 const bool face_rotation,
4457 const RefinementCase<2> &face_ref_case)
4458{
4459 const unsigned int dim = 3;
4460
4462 ExcMessage("Cell has no children."));
4463 AssertIndexRange(face, faces_per_cell);
4464 if (!(subface == 0 &&
4466 {
4467 AssertIndexRange(subface,
4468 GeometryInfo<dim - 1>::n_children(face_ref_case));
4469 }
4470
4471 // invalid number used for invalid cases,
4472 // e.g. when the children are more refined at
4473 // a given face than the face itself
4474 const unsigned int e = numbers::invalid_unsigned_int;
4475
4476 // the whole process of finding a child cell
4477 // at a given subface considering the
4478 // possibly anisotropic refinement cases of
4479 // the cell and the face as well as
4480 // orientation, flip and rotation of the face
4481 // is quite complicated. thus, we break it
4482 // down into several steps.
4483
4484 // first step: convert the given face refine
4485 // case to a face refine case concerning the
4486 // face in standard orientation (, flip and
4487 // rotation). This only affects cut_x and
4488 // cut_y
4489 const RefinementCase<dim - 1> flip[4] = {
4490 RefinementCase<dim - 1>::no_refinement,
4491 RefinementCase<dim - 1>::cut_y,
4492 RefinementCase<dim - 1>::cut_x,
4493 RefinementCase<dim - 1>::cut_xy};
4494 // for face_orientation, 'true' is the
4495 // default value whereas for face_rotation,
4496 // 'false' is standard. If
4497 // <tt>face_rotation==face_orientation</tt>,
4498 // then one of them is non-standard and we
4499 // have to swap cut_x and cut_y, otherwise no
4500 // change is necessary.
4501 const RefinementCase<dim - 1> std_face_ref =
4502 (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4503
4504 // second step: convert the given subface
4505 // index to the one for a standard face
4506 // respecting face_orientation, face_flip and
4507 // face_rotation
4508
4509 // first index: face_ref_case
4510 // second index: face_orientation
4511 // third index: face_flip
4512 // forth index: face_rotation
4513 // fifth index: subface index
4514 const unsigned int subface_exchange[4][2][2][2][4] = {
4515 // no_refinement (subface 0 stays 0,
4516 // all others are invalid)
4517 {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4518 {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4519 // cut_x (here, if the face is only
4520 // rotated OR only falsely oriented,
4521 // then subface 0 of the non-standard
4522 // face does NOT correspond to one of
4523 // the subfaces of a standard
4524 // face. Thus we indicate the subface
4525 // which is located at the lower left
4526 // corner (the origin of the face's
4527 // local coordinate system) with
4528 // '0'. The rest of this issue is
4529 // taken care of using the above
4530 // conversion to a 'standard face
4531 // refine case')
4532 {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4533 {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4534 // cut_y (the same applies as for
4535 // cut_x)
4536 {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4537 {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4538 // cut_xyz: this information is
4539 // identical to the information
4540 // returned by
4541 // GeometryInfo<3>::real_to_standard_face_vertex()
4542 {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4543 // face_rotation=false, subfaces 0,1,2,3
4544 {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4545 // face_rotation=true, subfaces 0,1,2,3
4546 {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4547 // face_rotation=false, subfaces 0,1,2,3
4548 {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4549 // face_rotation=true, subfaces 0,1,2,3
4550 {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4551 // face_rotation=false, subfaces 0,1,2,3
4552 {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4553 // face_rotation=true, subfaces 0,1,2,3
4554 {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4555 // face_rotation=false, subfaces 0,1,2,3
4556 {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4557 // face_rotation=true, subfaces 0,1,2,3
4558
4559 const unsigned int std_subface =
4560 subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4561 [subface];
4562 Assert(std_subface != e, ExcInternalError());
4563
4564 // third step: these are the children, which
4565 // can be found at the given subfaces of an
4566 // isotropically refined (standard) face
4567 //
4568 // first index: (refinement_case-1)
4569 // second index: face_index
4570 // third index: subface_index (isotropic refinement)
4571 const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4572 [max_children_per_face] = {
4573 // cut_x
4574 {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4575 {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4576 {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4577 {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4578 {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4579 {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4580 // cut_y
4581 {{0, 1, 0, 1},
4582 {0, 1, 0, 1},
4583 {0, 0, 0, 0},
4584 {1, 1, 1, 1},
4585 {0, 0, 1, 1},
4586 {0, 0, 1, 1}},
4587 // cut_xy
4588 {{0, 2, 0, 2},
4589 {1, 3, 1, 3},
4590 {0, 0, 1, 1},
4591 {2, 2, 3, 3},
4592 {0, 1, 2, 3},
4593 {0, 1, 2, 3}},
4594 // cut_z
4595 {{0, 0, 1, 1},
4596 {0, 0, 1, 1},
4597 {0, 1, 0, 1},
4598 {0, 1, 0, 1},
4599 {0, 0, 0, 0},
4600 {1, 1, 1, 1}},
4601 // cut_xz
4602 {{0, 0, 1, 1},
4603 {2, 2, 3, 3},
4604 {0, 1, 2, 3},
4605 {0, 1, 2, 3},
4606 {0, 2, 0, 2},
4607 {1, 3, 1, 3}},
4608 // cut_yz
4609 {{0, 1, 2, 3},
4610 {0, 1, 2, 3},
4611 {0, 2, 0, 2},
4612 {1, 3, 1, 3},
4613 {0, 0, 1, 1},
4614 {2, 2, 3, 3}},
4615 // cut_xyz
4616 {{0, 2, 4, 6},
4617 {1, 3, 5, 7},
4618 {0, 4, 1, 5},
4619 {2, 6, 3, 7},
4620 {0, 1, 2, 3},
4621 {4, 5, 6, 7}}};
4622
4623 // forth step: check, whether the given face
4624 // refine case is valid for the given cell
4625 // refine case. this is the case, if the
4626 // given face refine case is at least as
4627 // refined as the face is for the given cell
4628 // refine case
4629
4630 // note, that we are considering standard
4631 // face refinement cases here and thus must
4632 // not pass the given orientation, flip and
4633 // rotation flags
4634 if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4635 face_refinement_case(ref_case, face))
4636 {
4637 // all is fine. for anisotropic face
4638 // refine cases, select one of the
4639 // isotropic subfaces which neighbors the
4640 // same child
4641
4642 // first index: (standard) face refine case
4643 // second index: subface index
4644 const unsigned int equivalent_iso_subface[4][4] = {
4645 {0, e, e, e}, // no_refinement
4646 {0, 3, e, e}, // cut_x
4647 {0, 3, e, e}, // cut_y
4648 {0, 1, 2, 3}}; // cut_xy
4649
4650 const unsigned int equ_std_subface =
4651 equivalent_iso_subface[std_face_ref][std_subface];
4652 Assert(equ_std_subface != e, ExcInternalError());
4653
4654 return iso_children[ref_case - 1][face][equ_std_subface];
4655 }
4656 else
4657 {
4658 // the face_ref_case was too coarse,
4659 // throw an error
4660 Assert(false,
4661 ExcMessage("The face RefineCase is too coarse "
4662 "for the given cell RefineCase."));
4663 }
4664 // we only get here in case of an error
4665 return e;
4666}
4667
4668
4669
4670template <>
4671inline unsigned int
4673 const unsigned int,
4674 const unsigned int,
4675 const bool,
4676 const bool,
4677 const bool,
4678 const RefinementCase<3> &)
4679{
4682}
4683
4684
4685
4686template <>
4687inline unsigned int
4688GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4689 const unsigned int line,
4690 const bool,
4691 const bool,
4692 const bool)
4693{
4694 (void)line;
4695 AssertIndexRange(face, faces_per_cell);
4696 AssertIndexRange(line, lines_per_face);
4697
4698 // The face is a line itself.
4699 return face;
4700}
4701
4702
4703
4704template <>
4705inline unsigned int
4706GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4707 const unsigned int line,
4708 const bool face_orientation,
4709 const bool face_flip,
4710 const bool face_rotation)
4711{
4712 AssertIndexRange(face, faces_per_cell);
4713 AssertIndexRange(line, lines_per_face);
4714
4715 const unsigned lines[faces_per_cell][lines_per_face] = {
4716 {8, 10, 0, 4}, // left face
4717 {9, 11, 1, 5}, // right face
4718 {2, 6, 8, 9}, // front face
4719 {3, 7, 10, 11}, // back face
4720 {0, 1, 2, 3}, // bottom face
4721 {4, 5, 6, 7}}; // top face
4722 return lines[face][real_to_standard_face_line(
4723 line, face_orientation, face_flip, face_rotation)];
4724}
4725
4726
4727
4728inline unsigned int
4729GeometryInfo<0>::face_to_cell_lines(const unsigned int,
4730 const unsigned int,
4731 const bool,
4732 const bool,
4733 const bool)
4734{
4737}
4738
4739
4740
4741template <int dim>
4742inline unsigned int
4743GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4744 const unsigned int,
4745 const bool,
4746 const bool,
4747 const bool)
4748{
4751}
4752
4753
4754
4755template <int dim>
4756inline unsigned int
4757GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4758 const unsigned int vertex,
4759 const bool face_orientation,
4760 const bool face_flip,
4761 const bool face_rotation)
4762{
4763 return child_cell_on_face(RefinementCase<dim>::isotropic_refinement,
4764 face,
4765 vertex,
4766 face_orientation,
4767 face_flip,
4768 face_rotation);
4769}
4770
4771
4772
4773inline unsigned int
4775 const unsigned int,
4776 const bool,
4777 const bool,
4778 const bool)
4779{
4782}
4783
4784
4785
4786template <int dim>
4787template <typename Number>
4788inline Point<dim, Number>
4790{
4792 for (unsigned int i = 0; i < dim; ++i)
4793 p[i] = std::min(std::max(q[i], Number(0.)), Number(1.));
4794
4795 return p;
4796}
4797
4798
4799
4800template <int dim>
4801inline double
4803{
4804 double result = 0.0;
4805
4806 for (unsigned int i = 0; i < dim; ++i)
4807 {
4808 result = std::max(result, -p[i]);
4809 result = std::max(result, p[i] - 1.);
4810 }
4811
4812 return result;
4813}
4814
4815
4816
4817template <int dim>
4818inline double
4820 const unsigned int i)
4821{
4823
4824 switch (dim)
4825 {
4826 case 1:
4827 {
4828 const double x = xi[0];
4829 switch (i)
4830 {
4831 case 0:
4832 return 1 - x;
4833 case 1:
4834 return x;
4835 }
4836 break;
4837 }
4838
4839 case 2:
4840 {
4841 const double x = xi[0];
4842 const double y = xi[1];
4843 switch (i)
4844 {
4845 case 0:
4846 return (1 - x) * (1 - y);
4847 case 1:
4848 return x * (1 - y);
4849 case 2:
4850 return (1 - x) * y;
4851 case 3:
4852 return x * y;
4853 }
4854 break;
4855 }
4856
4857 case 3:
4858 {
4859 const double x = xi[0];
4860 const double y = xi[1];
4861 const double z = xi[2];
4862 switch (i)
4863 {
4864 case 0:
4865 return (1 - x) * (1 - y) * (1 - z);
4866 case 1:
4867 return x * (1 - y) * (1 - z);
4868 case 2:
4869 return (1 - x) * y * (1 - z);
4870 case 3:
4871 return x * y * (1 - z);
4872 case 4:
4873 return (1 - x) * (1 - y) * z;
4874 case 5:
4875 return x * (1 - y) * z;
4876 case 6:
4877 return (1 - x) * y * z;
4878 case 7:
4879 return x * y * z;
4880 }
4881 break;
4882 }
4883
4884 default:
4886 }
4887 return -1e9;
4888}
4889
4890
4891
4892template <>
4894 const Point<1> &,
4895 const unsigned int i)
4896{
4898
4899 switch (i)
4900 {
4901 case 0:
4902 return Point<1>(-1.);
4903 case 1:
4904 return Point<1>(1.);
4905 }
4906
4907 return Point<1>(-1e9);
4908}
4909
4910
4911
4912template <>
4914 const Point<2> &xi,
4915 const unsigned int i)
4916{
4918
4919 const double x = xi[0];
4920 const double y = xi[1];
4921 switch (i)
4922 {
4923 case 0:
4924 return Point<2>(-(1 - y), -(1 - x));
4925 case 1:
4926 return Point<2>(1 - y, -x);
4927 case 2:
4928 return Point<2>(-y, 1 - x);
4929 case 3:
4930 return Point<2>(y, x);
4931 }
4932 return Point<2>(-1e9, -1e9);
4933}
4934
4935
4936
4937template <>
4939 const Point<3> &xi,
4940 const unsigned int i)
4941{
4943
4944 const double x = xi[0];
4945 const double y = xi[1];
4946 const double z = xi[2];
4947 switch (i)
4948 {
4949 case 0:
4950 return Point<3>(-(1 - y) * (1 - z),
4951 -(1 - x) * (1 - z),
4952 -(1 - x) * (1 - y));
4953 case 1:
4954 return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4955 case 2:
4956 return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4957 case 3:
4958 return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4959 case 4:
4960 return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4961 case 5:
4962 return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4963 case 6:
4964 return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4965 case 7:
4966 return Point<3>(y * z, x * z, x * y);
4967 }
4968
4969 return Point<3>(-1e9, -1e9, -1e9);
4970}
4971
4972
4973
4974template <int dim>
4975inline Tensor<1, dim>
4977 const unsigned int)
4978{
4980 return Tensor<1, dim>();
4981}
4982
4983
4984unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4985{
4986 return 0;
4987}
4988
4989
4990namespace internal
4991{
4992 namespace GeometryInfoHelper
4993 {
4994 // wedge product of a single
4995 // vector in 2d: we just have to
4996 // rotate it by 90 degrees to the
4997 // right
4998 inline Tensor<1, 2>
4999 wedge_product(const Tensor<1, 2> (&derivative)[1])
5000 {
5001 Tensor<1, 2> result;
5002 result[0] = derivative[0][1];
5003 result[1] = -derivative[0][0];
5004
5005 return result;
5006 }
5007
5008
5009 // wedge product of 2 vectors in
5010 // 3d is the cross product
5011 inline Tensor<1, 3>
5012 wedge_product(const Tensor<1, 3> (&derivative)[2])
5013 {
5014 return cross_product_3d(derivative[0], derivative[1]);
5015 }
5016
5017
5018 // wedge product of dim vectors
5019 // in dim-d: that's the
5020 // determinant of the matrix
5021 template <int dim>
5022 inline Tensor<0, dim>
5023 wedge_product(const Tensor<1, dim> (&derivative)[dim])
5024 {
5025 Tensor<2, dim> jacobian;
5026 for (unsigned int i = 0; i < dim; ++i)
5027 jacobian[i] = derivative[i];
5028
5029 return determinant(jacobian);
5030 }
5031 } // namespace GeometryInfoHelper
5032} // namespace internal
5033
5034
5035template <int dim>
5036template <int spacedim>
5037inline void
5039# ifndef DEAL_II_CXX14_CONSTEXPR_BUG
5040 (const Point<spacedim> (&vertices)[vertices_per_cell],
5041 Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell])
5042# else
5043 (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms)
5044# endif
5045{
5046 // for each of the vertices,
5047 // compute the alternating form
5048 // of the mapped unit
5049 // vectors. consider for
5050 // example the case of a quad
5051 // in spacedim==3: to do so, we
5052 // need to see how the
5053 // infinitesimal vectors
5054 // (d\xi_1,0) and (0,d\xi_2)
5055 // are transformed into
5056 // spacedim-dimensional space
5057 // and then form their cross
5058 // product (i.e. the wedge product
5059 // of two vectors). to this end, note
5060 // that
5061 // \vec x = sum_i \vec v_i phi_i(\vec xi)
5062 // so the transformed vectors are
5063 // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
5064 // and
5065 // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
5066 // which boils down to the columns
5067 // of the 3x2 matrix \grad_\xi x(\xi)
5068 //
5069 // a similar reasoning would
5070 // hold for all dim,spacedim
5071 // pairs -- we only have to
5072 // compute the wedge product of
5073 // the columns of the
5074 // derivatives
5075 for (unsigned int i = 0; i < vertices_per_cell; ++i)
5076 {
5077 Tensor<1, spacedim> derivatives[dim];
5078
5079 for (unsigned int j = 0; j < vertices_per_cell; ++j)
5080 {
5081 const Tensor<1, dim> grad_phi_j =
5082 d_linear_shape_function_gradient(unit_cell_vertex(i), j);
5083 for (unsigned int l = 0; l < dim; ++l)
5084 derivatives[l] += vertices[j] * grad_phi_j[l];
5085 }
5086
5087 forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
5088 }
5089}
5090
5091#endif // DOXYGEN
5093
5094#endif
GeometryPrimitive(const unsigned int object_dimension)
GeometryPrimitive(const Object object)
Definition point.h:111
RefinementCase operator|(const RefinementCase &r) const
static std::size_t memory_consumption()
RefinementCase operator~() const
RefinementCase operator&(const RefinementCase &r) const
RefinementCase(const typename RefinementPossibilities< dim >::Possibilities refinement_case)
static constexpr unsigned int n_refinement_cases
static RefinementCase cut_axis(const unsigned int i)
std::uint8_t value
static std::array< RefinementCase< dim >, n_refinement_cases > all_refinement_cases()
void serialize(Archive &ar, const unsigned int version)
RefinementCase(const std::uint8_t refinement_case)
static constexpr std::size_t memory_consumption()
SubfaceCase(const typename SubfacePossibilities< dim >::Possibilities subface_possibility)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 3 > vertices[4]
#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)
Definition exceptions.h:562
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcInvalidSubfaceCase(int arg1)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
spacedim const Point< spacedim > & p
Definition grid_tools.h:990
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
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)
static const unsigned int invalid_unsigned_int
Definition types.h:220
boost::integer_range< IncrementableType > iota_view
Definition iota_view.h:45
::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:107
#define DEAL_II_HOST_DEVICE
Definition numbers.h:34
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
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)