deal.II version GIT relicensing-1927-g3de9220933 2024-10-03 08:40:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
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
483template <int dim>
485{
520 enum Possibilities : std::uint8_t
521 {
526
539 };
540};
541
542
543
551template <>
553{
588 enum Possibilities : std::uint8_t
589 {
597 cut_x = 1,
602 };
603};
604
605
606
615template <>
617{
652 enum Possibilities : std::uint8_t
653 {
661 cut_x = 1,
665 cut_y = 2,
669 cut_xy = cut_x | cut_y,
670
674 isotropic_refinement = cut_xy
675 };
676};
677
678
679
688template <>
690{
725 enum Possibilities : std::uint8_t
726 {
734 cut_x = 1,
738 cut_y = 2,
742 cut_xy = cut_x | cut_y,
746 cut_z = 4,
750 cut_xz = cut_x | cut_z,
754 cut_yz = cut_y | cut_z,
758 cut_xyz = cut_x | cut_y | cut_z,
759
763 isotropic_refinement = cut_xyz
764 };
765};
766
767
768
777template <int dim>
779{
780public:
784 static constexpr unsigned int n_refinement_cases = (1 << dim);
785
790
796 const typename RefinementPossibilities<dim>::Possibilities refinement_case);
797
803 explicit RefinementCase(const std::uint8_t refinement_case);
804
816 DEAL_II_HOST_DEVICE operator std::uint8_t() const;
817
823 operator|(const RefinementCase &r) const;
824
830 operator&(const RefinementCase &r) const;
831
840 operator~() const;
841
842
848 static RefinementCase
849 cut_axis(const unsigned int i);
850
859 static std::array<RefinementCase<dim>, n_refinement_cases>
861
865 static std::size_t
867
873 template <class Archive>
874 void
875 serialize(Archive &ar, const unsigned int version);
876
882 int,
883 << "The refinement flags given (" << arg1
884 << ") contain set bits that do not "
885 << "make sense for the space dimension of the object to which they are applied.");
886
887private:
892 std::uint8_t value : (dim > 0 ? dim : 1);
893};
894
895
896
897namespace internal
898{
916 template <int dim>
918 {
923 {
928
932 case_isotropic = static_cast<std::uint8_t>(-1)
933 };
934 };
935
936
943 template <>
945 {
963 };
964
965
966
974 template <>
976 {
994 };
995
996
997
1006 template <>
1008 {
1016 {
1024 case_x = 1,
1028 case_isotropic = case_x
1030 };
1031
1032
1033
1123 template <>
1125 {
1133 {
1135 case_x = 1,
1136 case_x1y = 2,
1137 case_x2y = 3,
1138 case_x1y2y = 4,
1139 case_y = 5,
1140 case_y1x = 6,
1141 case_y2x = 7,
1142 case_y1x2x = 8,
1143 case_xy = 9,
1144
1145 case_isotropic = case_xy
1147 };
1148
1149
1150
1155 template <int dim>
1157 {
1158 public:
1165 subface_possibility);
1166
1178 operator std::uint8_t() const;
1179
1183 static constexpr std::size_t
1185
1191 int,
1192 << "The subface case given (" << arg1 << ") does not make sense "
1193 << "for the space dimension of the object to which they are applied.");
1194
1195 private:
1200 std::uint8_t value : (dim == 3 ? 4 : 1);
1201 };
1202
1203} // namespace internal
1204
1205
1206
1207template <int dim>
1208struct GeometryInfo;
1209
1210
1211
1231template <>
1233{
1241 static constexpr unsigned int max_children_per_cell = 1;
1242
1246 static constexpr unsigned int faces_per_cell = 0;
1247
1264 static std::array<unsigned int, 0>
1266
1274 static constexpr unsigned int max_children_per_face = 0;
1275
1281 static unsigned int
1282 n_children(const RefinementCase<0> &refinement_case);
1283
1287 static constexpr unsigned int vertices_per_cell = 1;
1288
1307 static std::array<unsigned int, vertices_per_cell>
1309
1333 static unsigned int
1334 face_to_cell_vertices(const unsigned int face,
1335 const unsigned int vertex,
1336 const bool face_orientation = true,
1337 const bool face_flip = false,
1338 const bool face_rotation = false);
1339
1354 static unsigned int
1355 face_to_cell_lines(const unsigned int face,
1356 const unsigned int line,
1357 const bool face_orientation = true,
1358 const bool face_flip = false,
1359 const bool face_rotation = false);
1360
1367 static constexpr unsigned int vertices_per_face = 0;
1368
1372 static constexpr unsigned int lines_per_face = 0;
1373
1377 static constexpr unsigned int quads_per_face = 0;
1378
1382 static constexpr unsigned int lines_per_cell = 0;
1383
1387 static constexpr unsigned int quads_per_cell = 0;
1388
1392 static constexpr unsigned int hexes_per_cell = 0;
1393
1411 static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1412
1426 static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1427};
1428
1429
1430
1962template <int dim>
1964{
1972 static constexpr unsigned int max_children_per_cell = 1 << dim;
1973
1977 static constexpr unsigned int faces_per_cell = 2 * dim;
1978
1997
2005 static constexpr unsigned int max_children_per_face =
2007
2011 static constexpr unsigned int vertices_per_cell = 1 << dim;
2012
2030
2034 static constexpr unsigned int vertices_per_face =
2036
2040 static constexpr unsigned int lines_per_face =
2042
2046 static constexpr unsigned int quads_per_face =
2048
2058 static constexpr unsigned int lines_per_cell =
2059 (2 * GeometryInfo<dim - 1>::lines_per_cell +
2061
2069 static constexpr unsigned int quads_per_cell =
2070 (2 * GeometryInfo<dim - 1>::quads_per_cell +
2072
2076 static constexpr unsigned int hexes_per_cell =
2077 (2 * GeometryInfo<dim - 1>::hexes_per_cell +
2079
2097 static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
2098 internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal();
2099
2113 static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
2114 internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal();
2115
2128 internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face();
2129
2134 static unsigned int
2135 n_children(const RefinementCase<dim> &refinement_case);
2136
2141 static unsigned int
2143
2153 static double
2155 const unsigned int subface_no);
2156
2162 static RefinementCase<dim - 1>
2163 face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2164 const unsigned int face_no,
2165 const bool face_orientation = true,
2166 const bool face_flip = false,
2167 const bool face_rotation = false);
2168
2174 static RefinementCase<dim>
2177 const unsigned int face_no,
2178 const bool face_orientation = true,
2179 const bool face_flip = false,
2180 const bool face_rotation = false);
2181
2186 static RefinementCase<1>
2187 line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2188 const unsigned int line_no);
2189
2194 static RefinementCase<dim>
2196
2243 static unsigned int
2245 const unsigned int face,
2246 const unsigned int subface,
2247 const bool face_orientation = true,
2248 const bool face_flip = false,
2249 const bool face_rotation = false,
2252
2266 static unsigned int
2267 line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
2268
2289 static unsigned int
2290 face_to_cell_vertices(const unsigned int face,
2291 const unsigned int vertex,
2292 const bool face_orientation = true,
2293 const bool face_flip = false,
2294 const bool face_rotation = false);
2295
2307 static unsigned int
2308 face_to_cell_lines(const unsigned int face,
2309 const unsigned int line,
2310 const bool face_orientation = true,
2311 const bool face_flip = false,
2312 const bool face_rotation = false);
2313
2323 static unsigned int
2324 standard_to_real_face_vertex(const unsigned int vertex,
2325 const bool face_orientation = true,
2326 const bool face_flip = false,
2327 const bool face_rotation = false);
2328
2338 static unsigned int
2339 real_to_standard_face_vertex(const unsigned int vertex,
2340 const bool face_orientation = true,
2341 const bool face_flip = false,
2342 const bool face_rotation = false);
2343
2353 static unsigned int
2354 standard_to_real_face_line(const unsigned int line,
2355 const bool face_orientation = true,
2356 const bool face_flip = false,
2357 const bool face_rotation = false);
2358
2364 static unsigned int
2365 standard_to_real_line_vertex(const unsigned int vertex,
2366 const bool line_orientation = true);
2367
2375 static std::array<unsigned int, 2>
2377
2385 static std::array<unsigned int, 2>
2387
2395 static std::array<unsigned int, 2>
2397
2407 static unsigned int
2408 real_to_standard_face_line(const unsigned int line,
2409 const bool face_orientation = true,
2410 const bool face_flip = false,
2411 const bool face_rotation = false);
2412
2418 static Point<dim>
2419 unit_cell_vertex(const unsigned int vertex);
2420
2430 static unsigned int
2432
2440 static Point<dim>
2442 const unsigned int child_index,
2443 const RefinementCase<dim> refine_case =
2445
2451 static Point<dim>
2453 const unsigned int child_index,
2454 const RefinementCase<dim> refine_case =
2456
2461 static bool
2463
2475 static bool
2476 is_inside_unit_cell(const Point<dim> &p, const double eps);
2477
2482 template <typename Number = double>
2483 static Point<dim, Number>
2485
2491 static double
2493
2498 static double
2499 d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2500
2505 static Tensor<1, dim>
2506 d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2507
2559 template <int spacedim>
2560 static void
2562#ifndef DEAL_II_CXX14_CONSTEXPR_BUG
2564 Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2565#else
2566 (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2567#endif
2568
2578 static constexpr std::array<unsigned int, faces_per_cell>
2580 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction();
2581
2598 static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2599 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation();
2600
2611 static constexpr std::array<Tensor<1, dim>, faces_per_cell>
2613 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_vector();
2614
2628 static constexpr ndarray<Tensor<1, dim>, faces_per_cell, dim - 1>
2629 unit_tangential_vectors = internal::GeometryInfoHelper::Initializers<
2631
2637 static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2638 internal::GeometryInfoHelper::Initializers<dim>::opposite_face();
2639
2640
2645 double,
2646 << "The coordinates must satisfy 0 <= x_i <= 1, "
2647 << "but here we have x_i=" << arg1);
2648
2653 int,
2654 int,
2655 int,
2656 << "RefinementCase<dim> " << arg1 << ": face " << arg2
2657 << " has no subface " << arg3);
2658};
2659
2660
2661
2662#ifndef DOXYGEN
2663
2664
2665/* -------------- declaration of explicit specializations ------------- */
2666
2667
2668template <>
2671 const unsigned int i);
2672template <>
2675 const unsigned int i);
2676template <>
2679 const unsigned int i);
2680
2681
2682
2683/* -------------- inline functions ------------- */
2684
2685
2686inline GeometryPrimitive::GeometryPrimitive(const Object object)
2687 : object(object)
2688{}
2689
2690
2691
2692inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2693 : object(static_cast<Object>(object_dimension))
2694{}
2695
2696
2697inline GeometryPrimitive::operator unsigned int() const
2698{
2699 return static_cast<unsigned int>(object);
2700}
2701
2702
2703
2704namespace internal
2705{
2706 template <int dim>
2708 const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2709 : value(subface_possibility)
2710 {}
2711
2712
2713 template <int dim>
2714 inline SubfaceCase<dim>::operator std::uint8_t() const
2715 {
2716 return value;
2717 }
2718
2719
2720} // namespace internal
2721
2722
2723template <int dim>
2725RefinementCase<dim>::cut_axis(const unsigned int)
2726{
2728 return static_cast<std::uint8_t>(-1);
2729}
2730
2731
2732template <>
2733inline RefinementCase<1>
2734RefinementCase<1>::cut_axis(const unsigned int i)
2735{
2736 AssertIndexRange(i, 1);
2737
2739 return options[i];
2740}
2741
2742
2743
2744template <>
2745inline RefinementCase<2>
2746RefinementCase<2>::cut_axis(const unsigned int i)
2747{
2748 AssertIndexRange(i, 2);
2749
2752 return options[i];
2753}
2754
2755
2756
2757template <>
2758inline RefinementCase<3>
2759RefinementCase<3>::cut_axis(const unsigned int i)
2760{
2761 AssertIndexRange(i, 3);
2762
2766 return options[i];
2767}
2768
2769
2770
2771template <>
2772inline std::array<RefinementCase<1>, 2>
2774{
2777}
2778
2779
2780
2781template <>
2782inline std::array<RefinementCase<2>, 4>
2784{
2789}
2790
2791
2792
2793template <>
2794inline std::array<RefinementCase<3>, 8>
2796{
2805}
2806
2807
2808
2809template <int dim>
2811 : value(RefinementPossibilities<dim>::no_refinement)
2812{}
2813
2814
2815
2816template <int dim>
2818 const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2819 : value(refinement_case)
2820{
2821 // check that only those bits of
2822 // the given argument are set that
2823 // make sense for a given space
2824 // dimension
2825 Assert((refinement_case &
2827 refinement_case,
2828 ExcInvalidRefinementCase(refinement_case));
2829}
2830
2831
2832
2833template <int dim>
2834inline RefinementCase<dim>::RefinementCase(const std::uint8_t refinement_case)
2835 : value(refinement_case)
2836{
2837 // check that only those bits of
2838 // the given argument are set that
2839 // make sense for a given space
2840 // dimension
2841 Assert((refinement_case &
2843 refinement_case,
2844 ExcInvalidRefinementCase(refinement_case));
2845}
2846
2847
2848
2849template <int dim>
2850inline DEAL_II_HOST_DEVICE RefinementCase<dim>::operator std::uint8_t() const
2851{
2852 return value;
2853}
2854
2855
2856
2857template <int dim>
2860{
2861 return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2862}
2863
2864
2865
2866template <int dim>
2869{
2870 return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2871}
2872
2873
2874
2875template <int dim>
2878{
2879 return RefinementCase<dim>(static_cast<std::uint8_t>(
2881}
2882
2883
2884
2885template <int dim>
2886inline std::size_t
2888{
2889 return sizeof(RefinementCase<dim>);
2890}
2891
2892
2893
2894template <int dim>
2895template <class Archive>
2896inline void
2897RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2898{
2899 // serialization can't deal with bitfields, so copy from/to a full sized
2900 // std::uint8_t
2901 std::uint8_t uchar_value = value;
2902 ar &uchar_value;
2903 value = uchar_value;
2904}
2905
2906
2907
2908template <>
2909inline Point<1>
2910GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2911{
2912 AssertIndexRange(vertex, vertices_per_cell);
2913
2914 return Point<1>(static_cast<double>(vertex));
2915}
2916
2917
2918
2919template <>
2920inline Point<2>
2921GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2922{
2923 AssertIndexRange(vertex, vertices_per_cell);
2924
2925 return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2926}
2927
2928
2929
2930template <>
2931inline Point<3>
2932GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2933{
2934 AssertIndexRange(vertex, vertices_per_cell);
2935
2936 return {static_cast<double>(vertex % 2),
2937 static_cast<double>(vertex / 2 % 2),
2938 static_cast<double>(vertex / 4)};
2939}
2940
2941
2942
2943inline std::array<unsigned int, 0>
2945{
2946 return {{}};
2947}
2948
2949
2950
2951inline std::array<unsigned int, 1>
2953{
2954 return {{0}};
2955}
2956
2957
2958
2959template <int dim>
2962{
2963 return {0U, faces_per_cell};
2964}
2965
2966
2967
2968template <int dim>
2971{
2972 return {0U, vertices_per_cell};
2973}
2974
2975
2976
2977template <int dim>
2978inline Point<dim>
2979GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2980{
2982
2983 return {};
2984}
2985
2986
2987
2988template <>
2989inline unsigned int
2991{
2992 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2993
2994 return (p[0] <= 0.5 ? 0 : 1);
2995}
2996
2997
2998
2999template <>
3000inline unsigned int
3002{
3003 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
3004 Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
3005
3006 return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
3007}
3008
3009
3010
3011template <>
3012inline unsigned int
3014{
3015 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
3016 Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
3017 Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
3018
3019 return (p[0] <= 0.5 ?
3020 (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
3021 (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
3022}
3023
3024
3025template <int dim>
3026inline unsigned int
3028{
3030
3031 return 0;
3032}
3033
3034
3035
3036template <>
3037inline Point<1>
3039 const unsigned int child_index,
3040 const RefinementCase<1> refine_case)
3041
3042{
3043 AssertIndexRange(child_index, 2);
3045 (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3046
3047 return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
3048}
3049
3050
3051
3052template <>
3053inline Point<2>
3055 const unsigned int child_index,
3056 const RefinementCase<2> refine_case)
3057
3058{
3059 AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3060
3061 Point<2> point = p;
3062 switch (refine_case)
3063 {
3065 point[0] *= 2.0;
3066 if (child_index == 1)
3067 point[0] -= 1.0;
3068 break;
3070 point[1] *= 2.0;
3071 if (child_index == 1)
3072 point[1] -= 1.0;
3073 break;
3075 point *= 2.0;
3076 point -= unit_cell_vertex(child_index);
3077 break;
3078 default:
3080 }
3081
3082 return point;
3083}
3084
3085
3086
3087template <>
3088inline Point<3>
3090 const unsigned int child_index,
3091 const RefinementCase<3> refine_case)
3092
3093{
3094 AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3095
3096 Point<3> point = p;
3097 // there might be a cleverer way to do
3098 // this, but since this function is called
3099 // in very few places for initialization
3100 // purposes only, I don't care at the
3101 // moment
3102 switch (refine_case)
3103 {
3105 point[0] *= 2.0;
3106 if (child_index == 1)
3107 point[0] -= 1.0;
3108 break;
3110 point[1] *= 2.0;
3111 if (child_index == 1)
3112 point[1] -= 1.0;
3113 break;
3115 point[2] *= 2.0;
3116 if (child_index == 1)
3117 point[2] -= 1.0;
3118 break;
3120 point[0] *= 2.0;
3121 point[1] *= 2.0;
3122 if (child_index % 2 == 1)
3123 point[0] -= 1.0;
3124 if (child_index / 2 == 1)
3125 point[1] -= 1.0;
3126 break;
3128 // careful, this is slightly
3129 // different from xy and yz due to
3130 // different internal numbering of
3131 // children!
3132 point[0] *= 2.0;
3133 point[2] *= 2.0;
3134 if (child_index / 2 == 1)
3135 point[0] -= 1.0;
3136 if (child_index % 2 == 1)
3137 point[2] -= 1.0;
3138 break;
3140 point[1] *= 2.0;
3141 point[2] *= 2.0;
3142 if (child_index % 2 == 1)
3143 point[1] -= 1.0;
3144 if (child_index / 2 == 1)
3145 point[2] -= 1.0;
3146 break;
3148 point *= 2.0;
3149 point -= unit_cell_vertex(child_index);
3150 break;
3151 default:
3153 }
3154
3155 return point;
3156}
3157
3158
3159
3160template <int dim>
3161inline Point<dim>
3163 const Point<dim> & /*p*/,
3164 const unsigned int /*child_index*/,
3165 const RefinementCase<dim> /*refine_case*/)
3166
3167{
3169 return {};
3170}
3171
3172
3173
3174template <>
3175inline Point<1>
3177 const unsigned int child_index,
3178 const RefinementCase<1> refine_case)
3179
3180{
3181 AssertIndexRange(child_index, 2);
3183 (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3184
3185 return (p + unit_cell_vertex(child_index)) * 0.5;
3186}
3187
3188
3189
3190template <>
3191inline Point<3>
3193 const unsigned int child_index,
3194 const RefinementCase<3> refine_case)
3195
3196{
3197 AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3198
3199 Point<3> point = p;
3200 // there might be a cleverer way to do
3201 // this, but since this function is called
3202 // in very few places for initialization
3203 // purposes only, I don't care at the
3204 // moment
3205 switch (refine_case)
3206 {
3208 if (child_index == 1)
3209 point[0] += 1.0;
3210 point[0] *= 0.5;
3211 break;
3213 if (child_index == 1)
3214 point[1] += 1.0;
3215 point[1] *= 0.5;
3216 break;
3218 if (child_index == 1)
3219 point[2] += 1.0;
3220 point[2] *= 0.5;
3221 break;
3223 if (child_index % 2 == 1)
3224 point[0] += 1.0;
3225 if (child_index / 2 == 1)
3226 point[1] += 1.0;
3227 point[0] *= 0.5;
3228 point[1] *= 0.5;
3229 break;
3231 // careful, this is slightly
3232 // different from xy and yz due to
3233 // different internal numbering of
3234 // children!
3235 if (child_index / 2 == 1)
3236 point[0] += 1.0;
3237 if (child_index % 2 == 1)
3238 point[2] += 1.0;
3239 point[0] *= 0.5;
3240 point[2] *= 0.5;
3241 break;
3243 if (child_index % 2 == 1)
3244 point[1] += 1.0;
3245 if (child_index / 2 == 1)
3246 point[2] += 1.0;
3247 point[1] *= 0.5;
3248 point[2] *= 0.5;
3249 break;
3251 point += unit_cell_vertex(child_index);
3252 point *= 0.5;
3253 break;
3254 default:
3256 }
3257
3258 return point;
3259}
3260
3261
3262
3263template <>
3264inline Point<2>
3266 const unsigned int child_index,
3267 const RefinementCase<2> refine_case)
3268{
3269 AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3270
3271 Point<2> point = p;
3272 switch (refine_case)
3273 {
3275 if (child_index == 1)
3276 point[0] += 1.0;
3277 point[0] *= 0.5;
3278 break;
3280 if (child_index == 1)
3281 point[1] += 1.0;
3282 point[1] *= 0.5;
3283 break;
3285 point += unit_cell_vertex(child_index);
3286 point *= 0.5;
3287 break;
3288 default:
3290 }
3291
3292 return point;
3293}
3294
3295
3296
3297template <int dim>
3298inline Point<dim>
3300 const Point<dim> & /*p*/,
3301 const unsigned int /*child_index*/,
3302 const RefinementCase<dim> /*refine_case*/)
3303{
3305 return {};
3306}
3307
3308
3309
3310template <int dim>
3311inline bool
3313{
3315 return false;
3316}
3317
3318template <>
3319inline bool
3321{
3322 return (p[0] >= 0.) && (p[0] <= 1.);
3323}
3324
3325
3326
3327template <>
3328inline bool
3330{
3331 return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
3332}
3333
3334
3335
3336template <>
3337inline bool
3339{
3340 return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
3341 (p[2] >= 0.) && (p[2] <= 1.);
3342}
3343
3344
3345
3346template <int dim>
3347inline bool
3349{
3351 return false;
3352}
3353
3354template <>
3355inline bool
3356GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
3357{
3358 return (p[0] >= -eps) && (p[0] <= 1. + eps);
3359}
3360
3361
3362
3363template <>
3364inline bool
3365GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
3366{
3367 const double l = -eps, u = 1 + eps;
3368 return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
3369}
3370
3371
3372
3373template <>
3374inline bool
3375GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
3376{
3377 const double l = -eps, u = 1.0 + eps;
3378 return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
3379 (p[2] >= l) && (p[2] <= u);
3380}
3381
3382
3383
3384template <>
3385inline unsigned int
3386GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
3387 const unsigned int vertex)
3388{
3389 (void)line;
3390 AssertIndexRange(line, lines_per_cell);
3391 AssertIndexRange(vertex, 2);
3392
3393 return vertex;
3394}
3395
3396
3397template <>
3398inline unsigned int
3399GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
3400 const unsigned int vertex)
3401{
3402 constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
3403 return cell_vertices[line][vertex];
3404}
3405
3406
3407
3408template <>
3409inline unsigned int
3410GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
3411 const unsigned int vertex)
3412{
3413 AssertIndexRange(line, lines_per_cell);
3414 AssertIndexRange(vertex, 2);
3415
3416 constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3417 {1, 3},
3418 {0, 1},
3419 {2, 3},
3420 {4, 6}, // top face
3421 {5, 7},
3422 {4, 5},
3423 {6, 7},
3424 {0,
3425 4}, // connects of bottom
3426 {1, 5}, // top face
3427 {2, 6},
3428 {3, 7}};
3429 return vertices[line][vertex];
3430}
3431
3432
3433template <>
3434inline unsigned int
3435GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3436{
3439}
3440
3441template <>
3442inline unsigned int
3443GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3444 const bool face_orientation,
3445 const bool face_flip,
3446 const bool face_rotation)
3447{
3449
3450 // set up a table to make sure that
3451 // we handle non-standard faces correctly
3452 //
3453 // so set up a table that for each vertex (of
3454 // a quad in standard position) describes
3455 // which vertex to take
3456 //
3457 // first index: four vertices 0...3
3458 //
3459 // second index: face_orientation; 0:
3460 // opposite normal, 1: standard
3461 //
3462 // third index: face_flip; 0: standard, 1:
3463 // face rotated by 180 degrees
3464 //
3465 // forth index: face_rotation: 0: standard,
3466 // 1: face rotated by 90 degrees
3467
3468 constexpr unsigned int vertex_translation[4][2][2][2] = {
3469 {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3470 // face_rotation=false and true
3471 {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3472 // face_rotation=false and true
3473 {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3474 // face_rotation=false and true
3475 {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3476 // face_rotation=false and true
3477
3478 {{{2, 3}, // vertex 1 ...
3479 {1, 0}},
3480 {{1, 0}, {2, 3}}},
3481
3482 {{{1, 0}, // vertex 2 ...
3483 {2, 3}},
3484 {{2, 3}, {1, 0}}},
3485
3486 {{{3, 1}, // vertex 3 ...
3487 {0, 2}},
3488 {{3, 1}, {0, 2}}}};
3489
3490 return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3491}
3492
3493
3494
3495template <int dim>
3496inline unsigned int
3497GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3498 const bool,
3499 const bool,
3500 const bool)
3501{
3502 Assert(dim > 1, ExcImpossibleInDim(dim));
3504 return vertex;
3505}
3506
3507template <int dim>
3508inline unsigned int
3510{
3511 constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3512 0, 2, 2, 4, 2, 4, 4, 8};
3513
3514 return n_children[ref_case];
3515}
3516
3517
3518
3519template <int dim>
3520inline unsigned int
3522{
3524 return 0;
3525}
3526
3527template <>
3528inline unsigned int
3530{
3531 Assert(false, ExcImpossibleInDim(1));
3532 return 0;
3533}
3534
3535template <>
3536inline unsigned int
3538{
3539 return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3540}
3541
3542
3543
3544template <>
3545inline unsigned int
3547{
3548 const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3549 0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3550 return nsubs[subface_case];
3551}
3552
3553
3554
3555template <int dim>
3556inline double
3558 const unsigned int)
3559{
3561 return 0.;
3562}
3563
3564template <>
3565inline double
3567 const unsigned int)
3568{
3569 return 1;
3570}
3571
3572
3573template <>
3574inline double
3576 const unsigned int)
3577{
3578 double ratio = 1;
3579 switch (subface_case)
3580 {
3582 // Here, an
3583 // Assert(false,ExcInternalError())
3584 // would be the right
3585 // choice, but
3586 // unfortunately the
3587 // current function is
3588 // also called for faces
3589 // without children (see
3590 // tests/fe/mapping.cc).
3591 // Assert(false, ExcMessage("Face has no subfaces."));
3592 // Furthermore, assign
3593 // following value as
3594 // otherwise the
3595 // bits/volume_x tests
3596 // break
3598 break;
3600 ratio = 0.5;
3601 break;
3602 default:
3603 // there should be no
3604 // cases left
3606 break;
3607 }
3608
3609 return ratio;
3610}
3611
3612
3613template <>
3614inline double
3616 const unsigned int subface_no)
3617{
3618 double ratio = 1;
3619 switch (subface_case)
3620 {
3622 // Here, an
3623 // Assert(false,ExcInternalError())
3624 // would be the right
3625 // choice, but
3626 // unfortunately the
3627 // current function is
3628 // also called for faces
3629 // without children (see
3630 // tests/bits/mesh_3d_16.cc). Add
3631 // following switch to
3632 // avoid diffs in
3633 // tests/bits/mesh_3d_16
3635 break;
3638 ratio = 0.5;
3639 break;
3643 ratio = 0.25;
3644 break;
3647 if (subface_no < 2)
3648 ratio = 0.25;
3649 else
3650 ratio = 0.5;
3651 break;
3654 if (subface_no == 0)
3655 ratio = 0.5;
3656 else
3657 ratio = 0.25;
3658 break;
3659 default:
3660 // there should be no
3661 // cases left
3663 break;
3664 }
3665
3666 return ratio;
3667}
3668
3669
3670
3671template <int dim>
3673 const RefinementCase<dim> &,
3674 const unsigned int,
3675 const bool,
3676 const bool,
3677 const bool)
3678{
3680 return RefinementCase<dim - 1>::no_refinement;
3681}
3682
3683template <>
3685 const RefinementCase<1> &,
3686 const unsigned int,
3687 const bool,
3688 const bool,
3689 const bool)
3690{
3691 Assert(false, ExcImpossibleInDim(1));
3692
3694}
3695
3696
3697template <>
3698inline RefinementCase<1>
3700 const RefinementCase<2> &cell_refinement_case,
3701 const unsigned int face_no,
3702 const bool,
3703 const bool,
3704 const bool)
3705{
3706 const unsigned int dim = 2;
3707 AssertIndexRange(cell_refinement_case,
3710
3711 // simple special case
3712 if (cell_refinement_case == RefinementCase<dim>::cut_xy)
3714
3715 const RefinementCase<dim - 1>
3718 {RefinementCase<dim - 1>::no_refinement, // no_refinement
3719 RefinementCase<dim - 1>::no_refinement},
3720
3721 {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3722
3723 {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3724
3725 {RefinementCase<dim - 1>::cut_x, // cut_xy
3726 RefinementCase<dim - 1>::cut_x}};
3727
3728 return ref_cases[cell_refinement_case][face_no / 2];
3729}
3730
3731
3732template <>
3733inline RefinementCase<2>
3735 const RefinementCase<3> &cell_refinement_case,
3736 const unsigned int face_no,
3737 const bool face_orientation,
3738 const bool /*face_flip*/,
3739 const bool face_rotation)
3740{
3741 const unsigned int dim = 3;
3742 AssertIndexRange(cell_refinement_case,
3745
3746 // simple special case
3747 if (cell_refinement_case == RefinementCase<dim>::cut_xyz)
3748 return RefinementCase<dim - 1>::cut_xy;
3749
3750 const RefinementCase<dim - 1>
3753 {RefinementCase<dim - 1>::no_refinement, // no_refinement
3754 RefinementCase<dim - 1>::no_refinement,
3755 RefinementCase<dim - 1>::no_refinement},
3756
3757 {RefinementCase<dim - 1>::no_refinement, // cut_x
3758 RefinementCase<dim - 1>::cut_y,
3759 RefinementCase<dim - 1>::cut_x},
3760
3761 {RefinementCase<dim - 1>::cut_x, // cut_y
3762 RefinementCase<dim - 1>::no_refinement,
3763 RefinementCase<dim - 1>::cut_y},
3764
3765 {RefinementCase<dim - 1>::cut_x, // cut_xy
3766 RefinementCase<dim - 1>::cut_y,
3767 RefinementCase<dim - 1>::cut_xy},
3768
3769 {RefinementCase<dim - 1>::cut_y, // cut_z
3770 RefinementCase<dim - 1>::cut_x,
3771 RefinementCase<dim - 1>::no_refinement},
3772
3773 {RefinementCase<dim - 1>::cut_y, // cut_xz
3774 RefinementCase<dim - 1>::cut_xy,
3775 RefinementCase<dim - 1>::cut_x},
3776
3777 {RefinementCase<dim - 1>::cut_xy, // cut_yz
3778 RefinementCase<dim - 1>::cut_x,
3779 RefinementCase<dim - 1>::cut_y},
3780
3781 {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3782 RefinementCase<dim - 1>::cut_xy,
3783 RefinementCase<dim - 1>::cut_xy},
3784 };
3785
3786 const RefinementCase<dim - 1> ref_case =
3787 ref_cases[cell_refinement_case][face_no / 2];
3788
3789 const RefinementCase<dim - 1> flip[4] = {
3790 RefinementCase<dim - 1>::no_refinement,
3791 RefinementCase<dim - 1>::cut_y,
3792 RefinementCase<dim - 1>::cut_x,
3793 RefinementCase<dim - 1>::cut_xy};
3794
3795 // correct the ref_case for face_orientation
3796 // and face_rotation. for face_orientation,
3797 // 'true' is the default value whereas for
3798 // face_rotation, 'false' is standard. If
3799 // <tt>face_rotation==face_orientation</tt>,
3800 // then one of them is non-standard and we
3801 // have to swap cut_x and cut_y, otherwise no
3802 // change is necessary. face_flip has no
3803 // influence. however, in order to keep the
3804 // interface consistent with other functions,
3805 // we still include it as an argument to this
3806 // function
3807 return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3808}
3809
3810
3811
3812template <int dim>
3813inline RefinementCase<1>
3815 const unsigned int)
3816{
3819}
3820
3821template <>
3822inline RefinementCase<1>
3824 const RefinementCase<1> &cell_refinement_case,
3825 const unsigned int line_no)
3826{
3827 (void)line_no;
3828 const unsigned int dim = 1;
3829 (void)dim;
3830 AssertIndexRange(cell_refinement_case,
3833
3834 return cell_refinement_case;
3835}
3836
3837
3838template <>
3839inline RefinementCase<1>
3841 const RefinementCase<2> &cell_refinement_case,
3842 const unsigned int line_no)
3843{
3844 // Assertions are in face_refinement_case()
3845 return face_refinement_case(cell_refinement_case, line_no);
3846}
3847
3848
3849template <>
3850inline RefinementCase<1>
3852 const RefinementCase<3> &cell_refinement_case,
3853 const unsigned int line_no)
3854{
3855 const unsigned int dim = 3;
3856 AssertIndexRange(cell_refinement_case,
3859
3860 // simple special case
3861 if (cell_refinement_case == RefinementCase<dim>::cut_xyz)
3863
3864 // array indicating, which simple refine
3865 // case cuts a line in direction x, y or
3866 // z. For example, cut_y and everything
3867 // containing cut_y (cut_xy, cut_yz,
3868 // cut_xyz) cuts lines, which are in y
3869 // direction.
3870 const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3873
3874 // order the direction of lines
3875 // 0->x, 1->y, 2->z
3876 const unsigned int direction[lines_per_cell] = {
3877 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3878
3879 return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3881 RefinementCase<1>::no_refinement);
3882}
3883
3884
3885
3886template <int dim>
3890 const unsigned int,
3891 const bool,
3892 const bool,
3893 const bool)
3894{
3896
3898}
3899
3900template <>
3901inline RefinementCase<1>
3903 const RefinementCase<0> &,
3904 const unsigned int,
3905 const bool,
3906 const bool,
3907 const bool)
3908{
3909 const unsigned int dim = 1;
3910 Assert(false, ExcImpossibleInDim(dim));
3911
3913}
3914
3915
3916template <>
3917inline RefinementCase<2>
3919 const RefinementCase<1> &face_refinement_case,
3920 const unsigned int face_no,
3921 const bool,
3922 const bool,
3923 const bool)
3924{
3925 const unsigned int dim = 2;
3926 AssertIndexRange(face_refinement_case,
3929
3930 if (face_refinement_case == RefinementCase<dim>::cut_x)
3931 return (face_no / 2) != 0u ? RefinementCase<dim>::cut_x :
3933 else
3935}
3936
3937
3938template <>
3939inline RefinementCase<3>
3941 const RefinementCase<2> &face_refinement_case,
3942 const unsigned int face_no,
3943 const bool face_orientation,
3944 const bool /*face_flip*/,
3945 const bool face_rotation)
3946{
3947 const unsigned int dim = 3;
3948 AssertIndexRange(face_refinement_case,
3951
3956
3957 // correct the face_refinement_case for
3958 // face_orientation and face_rotation. for
3959 // face_orientation, 'true' is the default
3960 // value whereas for face_rotation, 'false'
3961 // is standard. If
3962 // <tt>face_rotation==face_orientation</tt>,
3963 // then one of them is non-standard and we
3964 // have to swap cut_x and cut_y, otherwise no
3965 // change is necessary. face_flip has no
3966 // influence. however, in order to keep the
3967 // interface consistent with other functions,
3968 // we still include it as an argument to this
3969 // function
3970 const RefinementCase<dim - 1> std_face_ref =
3971 (face_orientation == face_rotation) ? flip[face_refinement_case] :
3972 face_refinement_case;
3973
3974 const RefinementCase<dim> face_to_cell[3][4] = {
3975 {RefinementCase<dim>::no_refinement, // faces 0 and 1
3976 RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3979
3980 {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3981 // "exchanged on faces 2 and 3")
3985
3986 {RefinementCase<dim>::no_refinement, // faces 4 and 5
3990
3991 return face_to_cell[face_no / 2][std_face_ref];
3992}
3993
3994
3995
3996template <int dim>
3999 const unsigned int)
4000{
4002
4004}
4005
4006template <>
4007inline RefinementCase<1>
4009 const unsigned int line_no)
4010{
4011 (void)line_no;
4012 AssertIndexRange(line_no, 1);
4013
4015}
4016
4017
4018template <>
4019inline RefinementCase<2>
4021 const unsigned int line_no)
4022{
4023 const unsigned int dim = 2;
4024 (void)dim;
4026
4027 return (line_no / 2) != 0u ? RefinementCase<2>::cut_x :
4029}
4030
4031
4032template <>
4033inline RefinementCase<3>
4035 const unsigned int line_no)
4036{
4037 const unsigned int dim = 3;
4039
4040 const RefinementCase<dim> ref_cases[6] = {
4041 RefinementCase<dim>::cut_y, // lines 0 and 1
4042 RefinementCase<dim>::cut_x, // lines 2 and 3
4043 RefinementCase<dim>::cut_y, // lines 4 and 5
4044 RefinementCase<dim>::cut_x, // lines 6 and 7
4045 RefinementCase<dim>::cut_z, // lines 8 and 9
4046 RefinementCase<dim>::cut_z}; // lines 10 and 11
4047
4048 return ref_cases[line_no / 2];
4049}
4050
4051
4052
4053template <>
4054inline unsigned int
4055GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
4056 const bool face_orientation,
4057 const bool face_flip,
4058 const bool face_rotation)
4059{
4061
4062 // set up a table to make sure that
4063 // we handle non-standard faces correctly
4064 //
4065 // so set up a table that for each vertex (of
4066 // a quad in standard position) describes
4067 // which vertex to take
4068 //
4069 // first index: four vertices 0...3
4070 //
4071 // second index: face_orientation; 0:
4072 // opposite normal, 1: standard
4073 //
4074 // third index: face_flip; 0: standard, 1:
4075 // face rotated by 180 degrees
4076 //
4077 // forth index: face_rotation: 0: standard,
4078 // 1: face rotated by 90 degrees
4079
4080 const unsigned int vertex_translation[4][2][2][2] = {
4081 {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
4082 // face_rotation=false and true
4083 {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
4084 // face_rotation=false and true
4085 {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
4086 // face_rotation=false and true
4087 {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
4088 // face_rotation=false and true
4089
4090 {{{2, 3}, // vertex 1 ...
4091 {1, 0}},
4092 {{1, 3}, {2, 0}}},
4093
4094 {{{1, 0}, // vertex 2 ...
4095 {2, 3}},
4096 {{2, 0}, {1, 3}}},
4097
4098 {{{3, 1}, // vertex 3 ...
4099 {0, 2}},
4100 {{3, 2}, {0, 1}}}};
4101
4102 return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
4103}
4104
4105
4106
4107template <int dim>
4108inline unsigned int
4109GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
4110 const bool,
4111 const bool,
4112 const bool)
4113{
4114 Assert(dim > 1, ExcImpossibleInDim(dim));
4116 return vertex;
4117}
4118
4119
4120
4121template <>
4122inline unsigned int
4123GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
4124 const bool face_orientation,
4125 const bool face_flip,
4126 const bool face_rotation)
4127{
4129
4130
4131 // make sure we handle
4132 // non-standard faces correctly
4133 //
4134 // so set up a table that for each line (of a
4135 // quad) describes which line to take
4136 //
4137 // first index: four lines 0...3
4138 //
4139 // second index: face_orientation; 0:
4140 // opposite normal, 1: standard
4141 //
4142 // third index: face_flip; 0: standard, 1:
4143 // face rotated by 180 degrees
4144 //
4145 // forth index: face_rotation: 0: standard,
4146 // 1: face rotated by 90 degrees
4147
4148 const unsigned int line_translation[4][2][2][2] = {
4149 {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4150 // face_rotation=false and true
4151 {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4152 // face_rotation=false and true
4153 {{0, 3}, // line 0, face_orientation=true, face_flip=false,
4154 // face_rotation=false and true
4155 {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
4156 // face_rotation=false and true
4157
4158 {{{3, 1}, // line 1 ...
4159 {2, 0}},
4160 {{1, 2}, {0, 3}}},
4161
4162 {{{0, 3}, // line 2 ...
4163 {1, 2}},
4164 {{2, 0}, {3, 1}}},
4165
4166 {{{1, 2}, // line 3 ...
4167 {0, 3}},
4168 {{3, 1}, {2, 0}}}};
4169
4170 return line_translation[line][face_orientation][face_flip][face_rotation];
4171}
4172
4173
4174
4175template <int dim>
4176inline unsigned int
4177GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
4178 const bool,
4179 const bool,
4180 const bool)
4181{
4183 return line;
4184}
4185
4186
4187
4188template <>
4189inline unsigned int
4190GeometryInfo<2>::standard_to_real_line_vertex(const unsigned int vertex,
4191 const bool line_orientation)
4192{
4193 return line_orientation ? vertex : (1 - vertex);
4194}
4195
4196
4197
4198template <int dim>
4199inline unsigned int
4200GeometryInfo<dim>::standard_to_real_line_vertex(const unsigned int vertex,
4201 const bool)
4202{
4204 return vertex;
4205}
4206
4207
4208
4209template <>
4210inline std::array<unsigned int, 2>
4212 const unsigned int vertex)
4213{
4214 return {{vertex % 2, vertex / 2}};
4215}
4216
4217
4218
4219template <int dim>
4220inline std::array<unsigned int, 2>
4222 const unsigned int vertex)
4223{
4225 (void)vertex;
4226 return {{0, 0}};
4227}
4228
4229
4230
4231template <>
4232inline std::array<unsigned int, 2>
4234{
4235 // set up a table that for each
4236 // line describes a) from which
4237 // quad to take it, b) which line
4238 // therein it is if the face is
4239 // oriented correctly
4240 static const unsigned int lookup_table[GeometryInfo<3>::lines_per_cell][2] = {
4241 {4, 0}, // take first four lines from bottom face
4242 {4, 1},
4243 {4, 2},
4244 {4, 3},
4245
4246 {5, 0}, // second four lines from top face
4247 {5, 1},
4248 {5, 2},
4249 {5, 3},
4250
4251 {0, 0}, // the rest randomly
4252 {1, 0},
4253 {0, 1},
4254 {1, 1}};
4255
4256 return {{lookup_table[i][0], lookup_table[i][1]}};
4257}
4258
4259
4260
4261template <int dim>
4262inline std::array<unsigned int, 2>
4264{
4266 (void)line;
4267 return {{0, 0}};
4268}
4269
4270
4271
4272template <>
4273inline std::array<unsigned int, 2>
4275 const unsigned int vertex)
4276{
4277 // get the corner indices by asking either the bottom or the top face for its
4278 // vertices. handle non-standard faces by calling the vertex reordering
4279 // function from GeometryInfo
4280
4281 // bottom face (4) for first four vertices, top face (5) for the rest
4282 return {{4 + vertex / 4, vertex % 4}};
4283}
4284
4285
4286
4287template <int dim>
4288inline std::array<unsigned int, 2>
4290 const unsigned int vertex)
4291{
4293 (void)vertex;
4294 return {{0, 0}};
4295}
4296
4297
4298
4299template <>
4300inline unsigned int
4301GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
4302 const bool face_orientation,
4303 const bool face_flip,
4304 const bool face_rotation)
4305{
4307
4308
4309 // make sure we handle
4310 // non-standard faces correctly
4311 //
4312 // so set up a table that for each line (of a
4313 // quad) describes which line to take
4314 //
4315 // first index: four lines 0...3
4316 //
4317 // second index: face_orientation; 0:
4318 // opposite normal, 1: standard
4319 //
4320 // third index: face_flip; 0: standard, 1:
4321 // face rotated by 180 degrees
4322 //
4323 // forth index: face_rotation: 0: standard,
4324 // 1: face rotated by 90 degrees
4325
4326 const unsigned int line_translation[4][2][2][2] = {
4327 {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4328 // face_rotation=false and true
4329 {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4330 // face_rotation=false and true
4331 {{0, 2}, // line 0, face_orientation=true, face_flip=false,
4332 // face_rotation=false and true
4333 {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
4334 // face_rotation=false and true
4335
4336 {{{3, 1}, // line 1 ...
4337 {2, 0}},
4338 {{1, 3}, {0, 2}}},
4339
4340 {{{0, 3}, // line 2 ...
4341 {1, 2}},
4342 {{2, 1}, {3, 0}}},
4343
4344 {{{1, 2}, // line 3 ...
4345 {0, 3}},
4346 {{3, 0}, {2, 1}}}};
4347
4348 return line_translation[line][face_orientation][face_flip][face_rotation];
4349}
4350
4351
4352
4353template <int dim>
4354inline unsigned int
4355GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
4356 const bool,
4357 const bool,
4358 const bool)
4359{
4361 return line;
4362}
4363
4364
4365
4366template <>
4367inline unsigned int
4369 const unsigned int face,
4370 const unsigned int subface,
4371 const bool,
4372 const bool,
4373 const bool,
4374 const RefinementCase<0> &)
4375{
4376 (void)subface;
4377 AssertIndexRange(face, faces_per_cell);
4378 AssertIndexRange(subface, max_children_per_face);
4379
4380 return face;
4381}
4382
4383
4384
4385template <>
4386inline unsigned int
4388 const unsigned int face,
4389 const unsigned int subface,
4390 const bool face_orientation,
4391 const bool /*face_flip*/,
4392 const bool /*face_rotation*/,
4393 const RefinementCase<1> &)
4394{
4395 AssertIndexRange(face, faces_per_cell);
4396 AssertIndexRange(subface, max_children_per_face);
4397
4398 // always return the child adjacent to the specified
4399 // subface. if the face of a cell is not refined, don't
4400 // throw an assertion but deliver the child adjacent to
4401 // the face nevertheless, i.e. deliver the child of
4402 // this cell adjacent to the subface of a possibly
4403 // refined neighbor. this simplifies setting neighbor
4404 // information in execute_refinement.
4405 constexpr unsigned int
4406 subcells[/* possible face orientation */ 2]
4407 [/* number of different ways to refine a cell */ 4]
4408 [/* faces_per_cell */ 4][/* max_children_per_face */ 2] = {
4409 {
4410 // Flipped orientation (face_orientation = false)
4411 {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
4412 {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
4413 {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
4414 },
4415 {
4416 // Normal orientation (face_orientation = true)
4417 {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
4418 {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
4419 {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
4420 }};
4421
4422 return subcells[face_orientation][ref_case - 1][face][subface];
4423}
4424
4425
4426
4427template <>
4428inline unsigned int
4430 const unsigned int face,
4431 const unsigned int subface,
4432 const bool face_orientation,
4433 const bool face_flip,
4434 const bool face_rotation,
4435 const RefinementCase<2> &face_ref_case)
4436{
4437 const unsigned int dim = 3;
4438
4440 ExcMessage("Cell has no children."));
4441 AssertIndexRange(face, faces_per_cell);
4442 if (!(subface == 0 &&
4444 {
4445 AssertIndexRange(subface,
4446 GeometryInfo<dim - 1>::n_children(face_ref_case));
4447 }
4448
4449 // invalid number used for invalid cases,
4450 // e.g. when the children are more refined at
4451 // a given face than the face itself
4452 const unsigned int e = numbers::invalid_unsigned_int;
4453
4454 // the whole process of finding a child cell
4455 // at a given subface considering the
4456 // possibly anisotropic refinement cases of
4457 // the cell and the face as well as
4458 // orientation, flip and rotation of the face
4459 // is quite complicated. thus, we break it
4460 // down into several steps.
4461
4462 // first step: convert the given face refine
4463 // case to a face refine case concerning the
4464 // face in standard orientation (, flip and
4465 // rotation). This only affects cut_x and
4466 // cut_y
4467 const RefinementCase<dim - 1> flip[4] = {
4468 RefinementCase<dim - 1>::no_refinement,
4469 RefinementCase<dim - 1>::cut_y,
4470 RefinementCase<dim - 1>::cut_x,
4471 RefinementCase<dim - 1>::cut_xy};
4472 // for face_orientation, 'true' is the
4473 // default value whereas for face_rotation,
4474 // 'false' is standard. If
4475 // <tt>face_rotation==face_orientation</tt>,
4476 // then one of them is non-standard and we
4477 // have to swap cut_x and cut_y, otherwise no
4478 // change is necessary.
4479 const RefinementCase<dim - 1> std_face_ref =
4480 (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4481
4482 // second step: convert the given subface
4483 // index to the one for a standard face
4484 // respecting face_orientation, face_flip and
4485 // face_rotation
4486
4487 // first index: face_ref_case
4488 // second index: face_orientation
4489 // third index: face_flip
4490 // forth index: face_rotation
4491 // fifth index: subface index
4492 const unsigned int subface_exchange[4][2][2][2][4] = {
4493 // no_refinement (subface 0 stays 0,
4494 // all others are invalid)
4495 {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4496 {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4497 // cut_x (here, if the face is only
4498 // rotated OR only falsely oriented,
4499 // then subface 0 of the non-standard
4500 // face does NOT correspond to one of
4501 // the subfaces of a standard
4502 // face. Thus we indicate the subface
4503 // which is located at the lower left
4504 // corner (the origin of the face's
4505 // local coordinate system) with
4506 // '0'. The rest of this issue is
4507 // taken care of using the above
4508 // conversion to a 'standard face
4509 // refine case')
4510 {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4511 {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4512 // cut_y (the same applies as for
4513 // cut_x)
4514 {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4515 {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4516 // cut_xyz: this information is
4517 // identical to the information
4518 // returned by
4519 // GeometryInfo<3>::real_to_standard_face_vertex()
4520 {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4521 // face_rotation=false, subfaces 0,1,2,3
4522 {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4523 // face_rotation=true, subfaces 0,1,2,3
4524 {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4525 // face_rotation=false, subfaces 0,1,2,3
4526 {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4527 // face_rotation=true, subfaces 0,1,2,3
4528 {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4529 // face_rotation=false, subfaces 0,1,2,3
4530 {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4531 // face_rotation=true, subfaces 0,1,2,3
4532 {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4533 // face_rotation=false, subfaces 0,1,2,3
4534 {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4535 // face_rotation=true, subfaces 0,1,2,3
4536
4537 const unsigned int std_subface =
4538 subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4539 [subface];
4540 Assert(std_subface != e, ExcInternalError());
4541
4542 // third step: these are the children, which
4543 // can be found at the given subfaces of an
4544 // isotropically refined (standard) face
4545 //
4546 // first index: (refinement_case-1)
4547 // second index: face_index
4548 // third index: subface_index (isotropic refinement)
4549 const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4550 [max_children_per_face] = {
4551 // cut_x
4552 {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4553 {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4554 {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4555 {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4556 {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4557 {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4558 // cut_y
4559 {{0, 1, 0, 1},
4560 {0, 1, 0, 1},
4561 {0, 0, 0, 0},
4562 {1, 1, 1, 1},
4563 {0, 0, 1, 1},
4564 {0, 0, 1, 1}},
4565 // cut_xy
4566 {{0, 2, 0, 2},
4567 {1, 3, 1, 3},
4568 {0, 0, 1, 1},
4569 {2, 2, 3, 3},
4570 {0, 1, 2, 3},
4571 {0, 1, 2, 3}},
4572 // cut_z
4573 {{0, 0, 1, 1},
4574 {0, 0, 1, 1},
4575 {0, 1, 0, 1},
4576 {0, 1, 0, 1},
4577 {0, 0, 0, 0},
4578 {1, 1, 1, 1}},
4579 // cut_xz
4580 {{0, 0, 1, 1},
4581 {2, 2, 3, 3},
4582 {0, 1, 2, 3},
4583 {0, 1, 2, 3},
4584 {0, 2, 0, 2},
4585 {1, 3, 1, 3}},
4586 // cut_yz
4587 {{0, 1, 2, 3},
4588 {0, 1, 2, 3},
4589 {0, 2, 0, 2},
4590 {1, 3, 1, 3},
4591 {0, 0, 1, 1},
4592 {2, 2, 3, 3}},
4593 // cut_xyz
4594 {{0, 2, 4, 6},
4595 {1, 3, 5, 7},
4596 {0, 4, 1, 5},
4597 {2, 6, 3, 7},
4598 {0, 1, 2, 3},
4599 {4, 5, 6, 7}}};
4600
4601 // forth step: check, whether the given face
4602 // refine case is valid for the given cell
4603 // refine case. this is the case, if the
4604 // given face refine case is at least as
4605 // refined as the face is for the given cell
4606 // refine case
4607
4608 // note, that we are considering standard
4609 // face refinement cases here and thus must
4610 // not pass the given orientation, flip and
4611 // rotation flags
4612 if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4613 face_refinement_case(ref_case, face))
4614 {
4615 // all is fine. for anisotropic face
4616 // refine cases, select one of the
4617 // isotropic subfaces which neighbors the
4618 // same child
4619
4620 // first index: (standard) face refine case
4621 // second index: subface index
4622 const unsigned int equivalent_iso_subface[4][4] = {
4623 {0, e, e, e}, // no_refinement
4624 {0, 3, e, e}, // cut_x
4625 {0, 3, e, e}, // cut_y
4626 {0, 1, 2, 3}}; // cut_xy
4627
4628 const unsigned int equ_std_subface =
4629 equivalent_iso_subface[std_face_ref][std_subface];
4630 Assert(equ_std_subface != e, ExcInternalError());
4631
4632 return iso_children[ref_case - 1][face][equ_std_subface];
4633 }
4634 else
4635 {
4636 // the face_ref_case was too coarse,
4637 // throw an error
4638 Assert(false,
4639 ExcMessage("The face RefineCase is too coarse "
4640 "for the given cell RefineCase."));
4641 }
4642 // we only get here in case of an error
4643 return e;
4644}
4645
4646
4647
4648template <>
4649inline unsigned int
4651 const unsigned int,
4652 const unsigned int,
4653 const bool,
4654 const bool,
4655 const bool,
4656 const RefinementCase<3> &)
4657{
4660}
4661
4662
4663
4664template <>
4665inline unsigned int
4666GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4667 const unsigned int line,
4668 const bool,
4669 const bool,
4670 const bool)
4671{
4672 (void)line;
4673 AssertIndexRange(face, faces_per_cell);
4674 AssertIndexRange(line, lines_per_face);
4675
4676 // The face is a line itself.
4677 return face;
4678}
4679
4680
4681
4682template <>
4683inline unsigned int
4684GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4685 const unsigned int line,
4686 const bool face_orientation,
4687 const bool face_flip,
4688 const bool face_rotation)
4689{
4690 AssertIndexRange(face, faces_per_cell);
4691 AssertIndexRange(line, lines_per_face);
4692
4693 const unsigned lines[faces_per_cell][lines_per_face] = {
4694 {8, 10, 0, 4}, // left face
4695 {9, 11, 1, 5}, // right face
4696 {2, 6, 8, 9}, // front face
4697 {3, 7, 10, 11}, // back face
4698 {0, 1, 2, 3}, // bottom face
4699 {4, 5, 6, 7}}; // top face
4700 return lines[face][real_to_standard_face_line(
4701 line, face_orientation, face_flip, face_rotation)];
4702}
4703
4704
4705
4706inline unsigned int
4707GeometryInfo<0>::face_to_cell_lines(const unsigned int,
4708 const unsigned int,
4709 const bool,
4710 const bool,
4711 const bool)
4712{
4715}
4716
4717
4718
4719template <int dim>
4720inline unsigned int
4721GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4722 const unsigned int,
4723 const bool,
4724 const bool,
4725 const bool)
4726{
4729}
4730
4731
4732
4733template <int dim>
4734inline unsigned int
4735GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4736 const unsigned int vertex,
4737 const bool face_orientation,
4738 const bool face_flip,
4739 const bool face_rotation)
4740{
4741 return child_cell_on_face(RefinementCase<dim>::isotropic_refinement,
4742 face,
4743 vertex,
4744 face_orientation,
4745 face_flip,
4746 face_rotation);
4747}
4748
4749
4750
4751inline unsigned int
4753 const unsigned int,
4754 const bool,
4755 const bool,
4756 const bool)
4757{
4760}
4761
4762
4763
4764template <int dim>
4765template <typename Number>
4766inline Point<dim, Number>
4768{
4770 for (unsigned int i = 0; i < dim; ++i)
4771 p[i] = std::min(std::max(q[i], Number(0.)), Number(1.));
4772
4773 return p;
4774}
4775
4776
4777
4778template <int dim>
4779inline double
4781{
4782 double result = 0.0;
4783
4784 for (unsigned int i = 0; i < dim; ++i)
4785 {
4786 result = std::max(result, -p[i]);
4787 result = std::max(result, p[i] - 1.);
4788 }
4789
4790 return result;
4791}
4792
4793
4794
4795template <int dim>
4796inline double
4798 const unsigned int i)
4799{
4801
4802 switch (dim)
4803 {
4804 case 1:
4805 {
4806 const double x = xi[0];
4807 switch (i)
4808 {
4809 case 0:
4810 return 1 - x;
4811 case 1:
4812 return x;
4813 }
4814 break;
4815 }
4816
4817 case 2:
4818 {
4819 const double x = xi[0];
4820 const double y = xi[1];
4821 switch (i)
4822 {
4823 case 0:
4824 return (1 - x) * (1 - y);
4825 case 1:
4826 return x * (1 - y);
4827 case 2:
4828 return (1 - x) * y;
4829 case 3:
4830 return x * y;
4831 }
4832 break;
4833 }
4834
4835 case 3:
4836 {
4837 const double x = xi[0];
4838 const double y = xi[1];
4839 const double z = xi[2];
4840 switch (i)
4841 {
4842 case 0:
4843 return (1 - x) * (1 - y) * (1 - z);
4844 case 1:
4845 return x * (1 - y) * (1 - z);
4846 case 2:
4847 return (1 - x) * y * (1 - z);
4848 case 3:
4849 return x * y * (1 - z);
4850 case 4:
4851 return (1 - x) * (1 - y) * z;
4852 case 5:
4853 return x * (1 - y) * z;
4854 case 6:
4855 return (1 - x) * y * z;
4856 case 7:
4857 return x * y * z;
4858 }
4859 break;
4860 }
4861
4862 default:
4864 }
4865 return -1e9;
4866}
4867
4868
4869
4870template <>
4872 const Point<1> &,
4873 const unsigned int i)
4874{
4876
4877 switch (i)
4878 {
4879 case 0:
4880 return Point<1>(-1.);
4881 case 1:
4882 return Point<1>(1.);
4883 }
4884
4885 return Point<1>(-1e9);
4886}
4887
4888
4889
4890template <>
4892 const Point<2> &xi,
4893 const unsigned int i)
4894{
4896
4897 const double x = xi[0];
4898 const double y = xi[1];
4899 switch (i)
4900 {
4901 case 0:
4902 return Point<2>(-(1 - y), -(1 - x));
4903 case 1:
4904 return Point<2>(1 - y, -x);
4905 case 2:
4906 return Point<2>(-y, 1 - x);
4907 case 3:
4908 return Point<2>(y, x);
4909 }
4910 return Point<2>(-1e9, -1e9);
4911}
4912
4913
4914
4915template <>
4917 const Point<3> &xi,
4918 const unsigned int i)
4919{
4921
4922 const double x = xi[0];
4923 const double y = xi[1];
4924 const double z = xi[2];
4925 switch (i)
4926 {
4927 case 0:
4928 return Point<3>(-(1 - y) * (1 - z),
4929 -(1 - x) * (1 - z),
4930 -(1 - x) * (1 - y));
4931 case 1:
4932 return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4933 case 2:
4934 return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4935 case 3:
4936 return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4937 case 4:
4938 return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4939 case 5:
4940 return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4941 case 6:
4942 return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4943 case 7:
4944 return Point<3>(y * z, x * z, x * y);
4945 }
4946
4947 return Point<3>(-1e9, -1e9, -1e9);
4948}
4949
4950
4951
4952template <int dim>
4953inline Tensor<1, dim>
4955 const unsigned int)
4956{
4958 return Tensor<1, dim>();
4959}
4960
4961
4962unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4963{
4964 return 0;
4965}
4966
4967
4968namespace internal
4969{
4970 namespace GeometryInfoHelper
4971 {
4972 // wedge product of a single
4973 // vector in 2d: we just have to
4974 // rotate it by 90 degrees to the
4975 // right
4976 inline Tensor<1, 2>
4977 wedge_product(const Tensor<1, 2> (&derivative)[1])
4978 {
4979 Tensor<1, 2> result;
4980 result[0] = derivative[0][1];
4981 result[1] = -derivative[0][0];
4982
4983 return result;
4984 }
4985
4986
4987 // wedge product of 2 vectors in
4988 // 3d is the cross product
4989 inline Tensor<1, 3>
4990 wedge_product(const Tensor<1, 3> (&derivative)[2])
4991 {
4992 return cross_product_3d(derivative[0], derivative[1]);
4993 }
4994
4995
4996 // wedge product of dim vectors
4997 // in dim-d: that's the
4998 // determinant of the matrix
4999 template <int dim>
5000 inline Tensor<0, dim>
5001 wedge_product(const Tensor<1, dim> (&derivative)[dim])
5002 {
5003 Tensor<2, dim> jacobian;
5004 for (unsigned int i = 0; i < dim; ++i)
5005 jacobian[i] = derivative[i];
5006
5007 return determinant(jacobian);
5008 }
5009 } // namespace GeometryInfoHelper
5010} // namespace internal
5011
5012
5013template <int dim>
5014template <int spacedim>
5015inline void
5017# ifndef DEAL_II_CXX14_CONSTEXPR_BUG
5018 (const Point<spacedim> (&vertices)[vertices_per_cell],
5019 Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell])
5020# else
5021 (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms)
5022# endif
5023{
5024 // for each of the vertices,
5025 // compute the alternating form
5026 // of the mapped unit
5027 // vectors. consider for
5028 // example the case of a quad
5029 // in spacedim==3: to do so, we
5030 // need to see how the
5031 // infinitesimal vectors
5032 // (d\xi_1,0) and (0,d\xi_2)
5033 // are transformed into
5034 // spacedim-dimensional space
5035 // and then form their cross
5036 // product (i.e. the wedge product
5037 // of two vectors). to this end, note
5038 // that
5039 // \vec x = sum_i \vec v_i phi_i(\vec xi)
5040 // so the transformed vectors are
5041 // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
5042 // and
5043 // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
5044 // which boils down to the columns
5045 // of the 3x2 matrix \grad_\xi x(\xi)
5046 //
5047 // a similar reasoning would
5048 // hold for all dim,spacedim
5049 // pairs -- we only have to
5050 // compute the wedge product of
5051 // the columns of the
5052 // derivatives
5053 for (unsigned int i = 0; i < vertices_per_cell; ++i)
5054 {
5055 Tensor<1, spacedim> derivatives[dim];
5056
5057 for (unsigned int j = 0; j < vertices_per_cell; ++j)
5058 {
5059 const Tensor<1, dim> grad_phi_j =
5060 d_linear_shape_function_gradient(unit_cell_vertex(i), j);
5061 for (unsigned int l = 0; l < dim; ++l)
5062 derivatives[l] += vertices[j] * grad_phi_j[l];
5063 }
5064
5065 forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
5066 }
5067}
5068
5069#endif // DOXYGEN
5071
5072#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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
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:557
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:511
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcInvalidSubfaceCase(int arg1)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
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:30
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 > &)