deal.II version GIT relicensing-2888-ge8e3b84039 2025-03-21 18:30: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#include <deal.II/base/tensor.h>
26
27#include <array>
28#include <cstdint>
29
30
31
33
34#ifndef DOXYGEN
35namespace internal
36{
37 namespace GeometryInfoHelper
38 {
39 // A struct that holds the values for all the arrays we want to initialize
40 // in GeometryInfo
41 template <int dim>
42 struct Initializers;
43
44 template <>
45 struct Initializers<1>
46 {
47 static constexpr std::array<unsigned int, 2>
48 ucd_to_deal()
49 {
50 return {{0, 1}};
51 }
52
53 static constexpr std::array<unsigned int, 2>
54 unit_normal_direction()
55 {
56 return {{0, 0}};
57 }
58
59 static constexpr std::array<int, 2>
60 unit_normal_orientation()
61 {
62 return {{-1, 1}};
63 }
64
65 static constexpr std::array<Tensor<1, 1>, 2>
66 unit_normal_vector()
67 {
68 return {{Tensor<1, 1>{{-1}}, Tensor<1, 1>{{1}}}};
69 }
70
71 static constexpr ::ndarray<Tensor<1, 1>, 2, 0>
72 unit_tangential_vectors()
73 {
74 return {{{{}}, {{}}}};
75 }
77 static constexpr std::array<unsigned int, 2>
78 opposite_face()
79 {
80 return {{1, 0}};
81 }
82
83 static constexpr std::array<unsigned int, 2>
84 dx_to_deal()
85 {
86 return {{0, 1}};
87 }
88
89 static constexpr ::ndarray<unsigned int, 2, 1>
90 vertex_to_face()
91 {
92 return {{{{0}}, {{1}}}};
93 }
94 };
95
96 template <>
97 struct Initializers<2>
98 {
99 static constexpr std::array<unsigned int, 4>
100 ucd_to_deal()
101 {
102 return {{0, 1, 3, 2}};
103 }
104
105 static constexpr std::array<unsigned int, 4>
106 unit_normal_direction()
107 {
108 return {{0, 0, 1, 1}};
109 }
110
111 static constexpr std::array<int, 4>
112 unit_normal_orientation()
113 {
114 return {{-1, 1, -1, 1}};
115 }
116
117 static constexpr std::array<Tensor<1, 2>, 4>
118 unit_normal_vector()
119 {
120 return {{Tensor<1, 2>{{-1., 0.}},
121 Tensor<1, 2>{{1., 0.}},
122 Tensor<1, 2>{{0., -1.}},
123 Tensor<1, 2>{{0., 1.}}}};
124 }
125
126 static constexpr ::ndarray<Tensor<1, 2>, 4, 1>
127 unit_tangential_vectors()
128 {
129 return {{{{Tensor<1, 2>{{0, -1}}}},
130 {{Tensor<1, 2>{{0, 1}}}},
131 {{Tensor<1, 2>{{1, 0}}}},
132 {{Tensor<1, 2>{{-1, 0}}}}}};
133 }
134
135 static constexpr std::array<unsigned int, 4>
136 opposite_face()
137 {
138 return {{1, 0, 3, 2}};
139 }
140
141 static constexpr std::array<unsigned int, 4>
142 dx_to_deal()
143 {
144 return {{0, 2, 1, 3}};
145 }
146
147 static constexpr ::ndarray<unsigned int, 4, 2>
148 vertex_to_face()
149 {
150 return {{{{0, 2}}, {{1, 2}}, {{0, 3}}, {{1, 3}}}};
151 }
152 };
153
154 template <>
155 struct Initializers<3>
156 {
157 static constexpr std::array<unsigned int, 8>
158 ucd_to_deal()
159 {
160 return {{0, 4, 5, 1, 2, 6, 7, 3}};
161 }
162
163 static constexpr std::array<unsigned int, 6>
164 unit_normal_direction()
165 {
166 return {{0, 0, 1, 1, 2, 2}};
167 }
168
169 static constexpr std::array<int, 6>
170 unit_normal_orientation()
171 {
172 return {{-1, 1, -1, 1, -1, 1}};
173 }
174
175 static constexpr std::array<Tensor<1, 3>, 6>
176 unit_normal_vector()
177 {
178 return {{Tensor<1, 3>{{-1, 0, 0}},
179 Tensor<1, 3>{{1, 0, 0}},
180 Tensor<1, 3>{{0, -1, 0}},
181 Tensor<1, 3>{{0, 1, 0}},
182 Tensor<1, 3>{{0, 0, -1}},
183 Tensor<1, 3>{{0, 0, 1}}}};
184 }
185
186 static constexpr ::ndarray<Tensor<1, 3>, 6, 2>
187 unit_tangential_vectors()
188 {
189 return {{{{Tensor<1, 3>{{0, -1, 0}}, Tensor<1, 3>{{0, 0, 1}}}},
190 {{Tensor<1, 3>{{0, 1, 0}}, Tensor<1, 3>{{0, 0, 1}}}},
191 {{Tensor<1, 3>{{0, 0, -1}}, Tensor<1, 3>{{1, 0, 0}}}},
192 {{Tensor<1, 3>{{0, 0, 1}}, Tensor<1, 3>{{1, 0, 0}}}},
193 {{Tensor<1, 3>{{-1, 0, 0}}, Tensor<1, 3>{{0, 1, 0}}}},
194 {{Tensor<1, 3>{{1, 0, 0}}, Tensor<1, 3>{{0, 1, 0}}}}}};
195 }
196
197 static constexpr std::array<unsigned int, 6>
198 opposite_face()
199 {
200 return {{1, 0, 3, 2, 5, 4}};
201 }
202
203 static constexpr std::array<unsigned int, 8>
204 dx_to_deal()
205 {
206 return {{0, 4, 2, 6, 1, 5, 3, 7}};
207 }
208
209 static constexpr ::ndarray<unsigned int, 8, 3>
210 vertex_to_face()
211 {
212 return {{{{0, 2, 4}},
213 {{1, 2, 4}},
214 {{0, 3, 4}},
215 {{1, 3, 4}},
216 {{0, 2, 5}},
217 {{1, 2, 5}},
218 {{0, 3, 5}},
219 {{1, 3, 5}}}};
220 }
221 };
222
223 template <>
224 struct Initializers<4>
225 {
226 static constexpr std::array<unsigned int, 16>
227 ucd_to_deal()
228 {
245 }
246
247 static constexpr std::array<unsigned int, 8>
248 unit_normal_direction()
249 {
250 return {{0, 0, 1, 1, 2, 2, 3, 3}};
251 }
252
253 static constexpr std::array<int, 8>
254 unit_normal_orientation()
255 {
256 return {{-1, 1, -1, 1, -1, 1, -1, 1}};
257 }
258
259 static constexpr std::array<Tensor<1, 4>, 8>
260 unit_normal_vector()
261 {
262 return {{Tensor<1, 4>{{-1, 0, 0, 0}},
263 Tensor<1, 4>{{1, 0, 0, 0}},
264 Tensor<1, 4>{{0, -1, 0, 0}},
265 Tensor<1, 4>{{0, 1, 0, 0}},
266 Tensor<1, 4>{{0, 0, -1, 0}},
267 Tensor<1, 4>{{0, 0, 1, 0}},
268 Tensor<1, 4>{{0, 0, 0, -1}},
269 Tensor<1, 4>{{0, 0, 0, 1}}}};
270 }
271
272 static constexpr ::ndarray<Tensor<1, 4>, 8, 3>
273 unit_tangential_vectors()
274 {
275 return {{{{Tensor<1, 4>{{0, -1, 0, 0}},
276 Tensor<1, 4>{{0, 0, 1, 0}},
277 Tensor<1, 4>{{0, 0, 0, 1}}}},
278 {{Tensor<1, 4>{{0, 1, 0, 0}},
279 Tensor<1, 4>{{0, 0, 1, 0}},
280 Tensor<1, 4>{{0, 0, 0, 1}}}},
281 {{Tensor<1, 4>{{0, 0, -1, 0}},
282 Tensor<1, 4>{{0, 0, 0, 1}},
283 Tensor<1, 4>{{1, 0, 0, 0}}}},
284 {{Tensor<1, 4>{{0, 0, 1, 0}},
285 Tensor<1, 4>{{0, 0, 0, 1}},
286 Tensor<1, 4>{{1, 0, 0, 0}}}},
287 {{Tensor<1, 4>{{0, 0, 0, -1}},
288 Tensor<1, 4>{{1, 0, 0, 0}},
289 Tensor<1, 4>{{0, 1, 0, 0}}}},
290 {{Tensor<1, 4>{{0, 0, 0, 1}},
291 Tensor<1, 4>{{1, 0, 0, 0}},
292 Tensor<1, 4>{{0, 1, 0, 0}}}},
293 {{Tensor<1, 4>{{-1, 0, 0, 0}},
294 Tensor<1, 4>{{0, 1, 0, 0}},
295 Tensor<1, 4>{{0, 0, 1, 0}}}},
296 {{Tensor<1, 4>{{1, 0, 0, 0}},
297 Tensor<1, 4>{{0, 1, 0, 0}},
298 Tensor<1, 4>{{0, 0, 1, 0}}}}}};
299 }
300
301 static constexpr std::array<unsigned int, 8>
302 opposite_face()
303 {
304 return {{1, 0, 3, 2, 5, 4, 7, 6}};
305 }
306
307 static constexpr std::array<unsigned int, 16>
308 dx_to_deal()
309 {
326 }
327
328 static constexpr ::ndarray<unsigned int, 16, 4>
329 vertex_to_face()
330 {
395 }
396 };
397 } // namespace GeometryInfoHelper
398} // namespace internal
399#endif // DOXYGEN
400
401
418{
419public:
427 {
435 line = 1,
439 quad = 2,
443 hex = 3
444 };
445
451
457 GeometryPrimitive(const unsigned int object_dimension);
458
463 operator unsigned int() const;
464
465private:
470};
471
472
473
484template <int dim>
486{
521 enum Possibilities : std::uint8_t
522 {
527
540 };
541};
542
543
544
552template <>
554{
589 enum Possibilities : std::uint8_t
590 {
598 cut_x = 1,
603 };
604};
605
606
607
616template <>
618{
653 enum Possibilities : std::uint8_t
654 {
662 cut_x = 1,
666 cut_y = 2,
670 cut_xy = cut_x | cut_y,
671
675 isotropic_refinement = cut_xy
676 };
677};
678
679
680
689template <>
691{
726 enum Possibilities : std::uint8_t
727 {
735 cut_x = 1,
739 cut_y = 2,
743 cut_xy = cut_x | cut_y,
747 cut_z = 4,
751 cut_xz = cut_x | cut_z,
755 cut_yz = cut_y | cut_z,
759 cut_xyz = cut_x | cut_y | cut_z,
760
764 isotropic_refinement = cut_xyz
765 };
766};
767
768
769
778template <int dim>
780{
781public:
785 static constexpr unsigned int n_refinement_cases = (1 << dim);
786
791
797 const typename RefinementPossibilities<dim>::Possibilities refinement_case);
798
804 explicit RefinementCase(const std::uint8_t refinement_case);
805
818 operator std::uint8_t() const;
819
825 operator|(const RefinementCase &r) const;
826
832 operator&(const RefinementCase &r) const;
833
842 operator~() const;
843
844
850 static RefinementCase
851 cut_axis(const unsigned int i);
852
861 static std::array<RefinementCase<dim>, n_refinement_cases>
863
867 static std::size_t
869
875 template <class Archive>
876 void
877 serialize(Archive &ar, const unsigned int version);
878
884 int,
885 << "The refinement flags given (" << arg1
886 << ") contain set bits that do not "
887 << "make sense for the space dimension of the object to which they are applied.");
888
889private:
894 std::uint8_t value : (dim > 0 ? dim : 1);
895};
896
897
898
899namespace internal
900{
918 template <int dim>
920 {
925 {
930
934 case_isotropic = static_cast<std::uint8_t>(-1)
935 };
936 };
937
938
945 template <>
947 {
965 };
966
967
968
976 template <>
978 {
996 };
997
998
999
1008 template <>
1010 {
1018 {
1026 case_x = 1,
1030 case_isotropic = case_x
1032 };
1033
1034
1035
1125 template <>
1127 {
1135 {
1137 case_x = 1,
1138 case_x1y = 2,
1139 case_x2y = 3,
1140 case_x1y2y = 4,
1141 case_y = 5,
1142 case_y1x = 6,
1143 case_y2x = 7,
1144 case_y1x2x = 8,
1145 case_xy = 9,
1146
1147 case_isotropic = case_xy
1149 };
1150
1151
1152
1157 template <int dim>
1159 {
1160 public:
1167 subface_possibility);
1168
1180 operator std::uint8_t() const;
1181
1185 static constexpr std::size_t
1187
1193 int,
1194 << "The subface case given (" << arg1 << ") does not make sense "
1195 << "for the space dimension of the object to which they are applied.");
1196
1197 private:
1202 std::uint8_t value : (dim == 3 ? 4 : 1);
1203 };
1204
1205} // namespace internal
1206
1207
1208
1209template <int dim>
1210struct GeometryInfo;
1211
1212
1213
1233template <>
1235{
1243 static constexpr unsigned int max_children_per_cell = 1;
1244
1248 static constexpr unsigned int faces_per_cell = 0;
1249
1266 static std::array<unsigned int, 0>
1268
1276 static constexpr unsigned int max_children_per_face = 0;
1277
1283 static unsigned int
1284 n_children(const RefinementCase<0> &refinement_case);
1285
1289 static constexpr unsigned int vertices_per_cell = 1;
1290
1309 static std::array<unsigned int, vertices_per_cell>
1311
1335 static unsigned int
1336 face_to_cell_vertices(const unsigned int face,
1337 const unsigned int vertex,
1338 const bool face_orientation = true,
1339 const bool face_flip = false,
1340 const bool face_rotation = false);
1341
1356 static unsigned int
1357 face_to_cell_lines(const unsigned int face,
1358 const unsigned int line,
1359 const bool face_orientation = true,
1360 const bool face_flip = false,
1361 const bool face_rotation = false);
1362
1369 static constexpr unsigned int vertices_per_face = 0;
1370
1374 static constexpr unsigned int lines_per_face = 0;
1375
1379 static constexpr unsigned int quads_per_face = 0;
1380
1384 static constexpr unsigned int lines_per_cell = 0;
1385
1389 static constexpr unsigned int quads_per_cell = 0;
1390
1394 static constexpr unsigned int hexes_per_cell = 0;
1395
1413 static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1414
1428 static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1429};
1430
1431
1432
1964template <int dim>
1966{
1974 static constexpr unsigned int max_children_per_cell = 1 << dim;
1975
1979 static constexpr unsigned int faces_per_cell = 2 * dim;
1980
1999
2007 static constexpr unsigned int max_children_per_face =
2009
2013 static constexpr unsigned int vertices_per_cell = 1 << dim;
2014
2032
2036 static constexpr unsigned int vertices_per_face =
2038
2042 static constexpr unsigned int lines_per_face =
2044
2048 static constexpr unsigned int quads_per_face =
2050
2060 static constexpr unsigned int lines_per_cell =
2061 (2 * GeometryInfo<dim - 1>::lines_per_cell +
2063
2071 static constexpr unsigned int quads_per_cell =
2072 (2 * GeometryInfo<dim - 1>::quads_per_cell +
2074
2078 static constexpr unsigned int hexes_per_cell =
2079 (2 * GeometryInfo<dim - 1>::hexes_per_cell +
2081
2099 static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
2100 internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal();
2101
2115 static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
2116 internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal();
2117
2130 internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face();
2131
2136 static unsigned int
2137 n_children(const RefinementCase<dim> &refinement_case);
2138
2143 static unsigned int
2145
2155 static double
2157 const unsigned int subface_no);
2158
2164 static RefinementCase<dim - 1>
2165 face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2166 const unsigned int face_no,
2167 const bool face_orientation = true,
2168 const bool face_flip = false,
2169 const bool face_rotation = false);
2170
2176 static RefinementCase<dim>
2179 const unsigned int face_no,
2180 const bool face_orientation = true,
2181 const bool face_flip = false,
2182 const bool face_rotation = false);
2183
2188 static RefinementCase<1>
2189 line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2190 const unsigned int line_no);
2191
2196 static RefinementCase<dim>
2198
2245 static unsigned int
2247 const unsigned int face,
2248 const unsigned int subface,
2249 const bool face_orientation = true,
2250 const bool face_flip = false,
2251 const bool face_rotation = false,
2254
2268 static unsigned int
2269 line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
2270
2291 static unsigned int
2292 face_to_cell_vertices(const unsigned int face,
2293 const unsigned int vertex,
2294 const bool face_orientation = true,
2295 const bool face_flip = false,
2296 const bool face_rotation = false);
2297
2309 static unsigned int
2310 face_to_cell_lines(const unsigned int face,
2311 const unsigned int line,
2312 const bool face_orientation = true,
2313 const bool face_flip = false,
2314 const bool face_rotation = false);
2315
2325 static unsigned int
2326 standard_to_real_face_vertex(const unsigned int vertex,
2327 const bool face_orientation = true,
2328 const bool face_flip = false,
2329 const bool face_rotation = false);
2330
2340 static unsigned int
2341 real_to_standard_face_vertex(const unsigned int vertex,
2342 const bool face_orientation = true,
2343 const bool face_flip = false,
2344 const bool face_rotation = false);
2345
2355 static unsigned int
2356 standard_to_real_face_line(const unsigned int line,
2357 const bool face_orientation = true,
2358 const bool face_flip = false,
2359 const bool face_rotation = false);
2360
2366 static unsigned int
2367 standard_to_real_line_vertex(const unsigned int vertex,
2368 const bool line_orientation = true);
2369
2377 static std::array<unsigned int, 2>
2379
2387 static std::array<unsigned int, 2>
2389
2397 static std::array<unsigned int, 2>
2399
2409 static unsigned int
2410 real_to_standard_face_line(const unsigned int line,
2411 const bool face_orientation = true,
2412 const bool face_flip = false,
2413 const bool face_rotation = false);
2414
2420 static Point<dim>
2421 unit_cell_vertex(const unsigned int vertex);
2422
2432 static unsigned int
2434
2442 static Point<dim>
2444 const unsigned int child_index,
2445 const RefinementCase<dim> refine_case =
2447
2453 static Point<dim>
2455 const unsigned int child_index,
2456 const RefinementCase<dim> refine_case =
2458
2463 static bool
2465
2477 static bool
2478 is_inside_unit_cell(const Point<dim> &p, const double eps);
2479
2484 template <typename Number = double>
2485 static Point<dim, Number>
2487
2493 static double
2495
2500 static double
2501 d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2502
2507 static Tensor<1, dim>
2508 d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2509
2561 template <int spacedim>
2562 static void
2564#ifndef DEAL_II_CXX14_CONSTEXPR_BUG
2565 (const Point<spacedim> (&vertices)[vertices_per_cell],
2566 Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2567#else
2568 (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2569#endif
2570
2580 static constexpr std::array<unsigned int, faces_per_cell>
2582 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction();
2583
2600 static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2601 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation();
2602
2613 static constexpr std::array<Tensor<1, dim>, faces_per_cell>
2615 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_vector();
2616
2630 static constexpr ndarray<Tensor<1, dim>, faces_per_cell, dim - 1>
2631 unit_tangential_vectors = internal::GeometryInfoHelper::Initializers<
2633
2639 static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2640 internal::GeometryInfoHelper::Initializers<dim>::opposite_face();
2641
2642
2647 double,
2648 << "The coordinates must satisfy 0 <= x_i <= 1, "
2649 << "but here we have x_i=" << arg1);
2650
2655 int,
2656 int,
2657 int,
2658 << "RefinementCase<dim> " << arg1 << ": face " << arg2
2659 << " has no subface " << arg3);
2660};
2661
2662
2663
2664#ifndef DOXYGEN
2665
2666
2667/* -------------- declaration of explicit specializations ------------- */
2668
2669
2670template <>
2673 const unsigned int i);
2674template <>
2677 const unsigned int i);
2678template <>
2681 const unsigned int i);
2682
2683
2684
2685/* -------------- inline functions ------------- */
2686
2687
2688inline GeometryPrimitive::GeometryPrimitive(const Object object)
2689 : object(object)
2690{}
2691
2692
2693
2694inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2695 : object(static_cast<Object>(object_dimension))
2696{}
2697
2698
2699inline GeometryPrimitive::operator unsigned int() const
2700{
2701 return static_cast<unsigned int>(object);
2702}
2703
2704
2705
2706namespace internal
2707{
2708 template <int dim>
2710 const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2711 : value(subface_possibility)
2712 {}
2713
2714
2715 template <int dim>
2716 inline SubfaceCase<dim>::operator std::uint8_t() const
2717 {
2718 return value;
2719 }
2720
2721
2722} // namespace internal
2723
2724
2725template <int dim>
2727RefinementCase<dim>::cut_axis(const unsigned int)
2728{
2730 return static_cast<std::uint8_t>(-1);
2731}
2732
2733
2734template <>
2735inline RefinementCase<1>
2736RefinementCase<1>::cut_axis(const unsigned int i)
2737{
2738 AssertIndexRange(i, 1);
2739
2741 return options[i];
2742}
2743
2744
2745
2746template <>
2747inline RefinementCase<2>
2748RefinementCase<2>::cut_axis(const unsigned int i)
2749{
2750 AssertIndexRange(i, 2);
2751
2754 return options[i];
2755}
2756
2757
2758
2759template <>
2760inline RefinementCase<3>
2761RefinementCase<3>::cut_axis(const unsigned int i)
2762{
2763 AssertIndexRange(i, 3);
2764
2768 return options[i];
2769}
2770
2771
2772
2773template <>
2774inline std::array<RefinementCase<1>, 2>
2776{
2779}
2780
2781
2782
2783template <>
2784inline std::array<RefinementCase<2>, 4>
2786{
2791}
2792
2793
2794
2795template <>
2796inline std::array<RefinementCase<3>, 8>
2798{
2807}
2808
2809
2810
2811template <int dim>
2813 : value(RefinementPossibilities<dim>::no_refinement)
2814{}
2815
2816
2817
2818template <int dim>
2820 const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2821 : value(refinement_case)
2822{
2823 // check that only those bits of
2824 // the given argument are set that
2825 // make sense for a given space
2826 // dimension
2827 Assert((refinement_case &
2829 refinement_case,
2830 ExcInvalidRefinementCase(refinement_case));
2831}
2832
2833
2834
2835template <int dim>
2836inline RefinementCase<dim>::RefinementCase(const std::uint8_t refinement_case)
2837 : value(refinement_case)
2838{
2839 // check that only those bits of
2840 // the given argument are set that
2841 // make sense for a given space
2842 // dimension
2843 Assert((refinement_case &
2845 refinement_case,
2846 ExcInvalidRefinementCase(refinement_case));
2847}
2848
2849
2850
2851template <int dim>
2852inline DEAL_II_HOST_DEVICE RefinementCase<dim>::operator std::uint8_t() const
2853{
2854 return value;
2855}
2856
2857
2858
2859template <int dim>
2862{
2863 return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2864}
2865
2866
2867
2868template <int dim>
2871{
2872 return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2873}
2874
2875
2876
2877template <int dim>
2880{
2881 return RefinementCase<dim>(static_cast<std::uint8_t>(
2883}
2884
2885
2886
2887template <int dim>
2888inline std::size_t
2890{
2891 return sizeof(RefinementCase<dim>);
2892}
2893
2894
2895
2896template <int dim>
2897template <class Archive>
2898inline void
2899RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2900{
2901 // serialization can't deal with bitfields, so copy from/to a full sized
2902 // std::uint8_t
2903 std::uint8_t uchar_value = value;
2904 ar &uchar_value;
2905 value = uchar_value;
2906}
2907
2908
2909
2910template <>
2911inline Point<1>
2912GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2913{
2914 AssertIndexRange(vertex, vertices_per_cell);
2915
2916 return Point<1>(static_cast<double>(vertex));
2917}
2918
2919
2920
2921template <>
2922inline Point<2>
2923GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2924{
2925 AssertIndexRange(vertex, vertices_per_cell);
2926
2927 return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2928}
2929
2930
2931
2932template <>
2933inline Point<3>
2934GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2935{
2936 AssertIndexRange(vertex, vertices_per_cell);
2937
2938 return {static_cast<double>(vertex % 2),
2939 static_cast<double>(vertex / 2 % 2),
2940 static_cast<double>(vertex / 4)};
2941}
2942
2943
2944
2945inline std::array<unsigned int, 0>
2947{
2948 return {{}};
2949}
2950
2951
2952
2953inline std::array<unsigned int, 1>
2955{
2956 return {{0}};
2957}
2958
2959
2960
2961template <int dim>
2964{
2966 0U, faces_per_cell);
2967}
2968
2969
2970
2971template <int dim>
2974{
2976 0U, vertices_per_cell);
2977}
2978
2979
2980
2981template <int dim>
2982inline Point<dim>
2983GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2984{
2986
2987 return {};
2988}
2989
2990
2991
2992template <>
2993inline unsigned int
2995{
2996 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2997
2998 return (p[0] <= 0.5 ? 0 : 1);
2999}
3000
3001
3002
3003template <>
3004inline unsigned int
3006{
3007 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
3008 Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
3009
3010 return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
3011}
3012
3013
3014
3015template <>
3016inline unsigned int
3018{
3019 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
3020 Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
3021 Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
3022
3023 return (p[0] <= 0.5 ?
3024 (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
3025 (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
3026}
3027
3028
3029template <int dim>
3030inline unsigned int
3032{
3034
3035 return 0;
3036}
3037
3038
3039
3040template <>
3041inline Point<1>
3043 const unsigned int child_index,
3044 const RefinementCase<1> refine_case)
3045
3046{
3047 AssertIndexRange(child_index, 2);
3049 (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3050
3051 return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
3052}
3053
3054
3055
3056template <>
3057inline Point<2>
3059 const unsigned int child_index,
3060 const RefinementCase<2> refine_case)
3061
3062{
3063 AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3064
3065 Point<2> point = p;
3066 switch (refine_case)
3067 {
3069 point[0] *= 2.0;
3070 if (child_index == 1)
3071 point[0] -= 1.0;
3072 break;
3074 point[1] *= 2.0;
3075 if (child_index == 1)
3076 point[1] -= 1.0;
3077 break;
3079 point *= 2.0;
3080 point -= unit_cell_vertex(child_index);
3081 break;
3082 default:
3084 }
3085
3086 return point;
3087}
3088
3089
3090
3091template <>
3092inline Point<3>
3094 const unsigned int child_index,
3095 const RefinementCase<3> refine_case)
3096
3097{
3098 AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3099
3100 Point<3> point = p;
3101 // there might be a cleverer way to do
3102 // this, but since this function is called
3103 // in very few places for initialization
3104 // purposes only, I don't care at the
3105 // moment
3106 switch (refine_case)
3107 {
3109 point[0] *= 2.0;
3110 if (child_index == 1)
3111 point[0] -= 1.0;
3112 break;
3114 point[1] *= 2.0;
3115 if (child_index == 1)
3116 point[1] -= 1.0;
3117 break;
3119 point[2] *= 2.0;
3120 if (child_index == 1)
3121 point[2] -= 1.0;
3122 break;
3124 point[0] *= 2.0;
3125 point[1] *= 2.0;
3126 if (child_index % 2 == 1)
3127 point[0] -= 1.0;
3128 if (child_index / 2 == 1)
3129 point[1] -= 1.0;
3130 break;
3132 // careful, this is slightly
3133 // different from xy and yz due to
3134 // different internal numbering of
3135 // children!
3136 point[0] *= 2.0;
3137 point[2] *= 2.0;
3138 if (child_index / 2 == 1)
3139 point[0] -= 1.0;
3140 if (child_index % 2 == 1)
3141 point[2] -= 1.0;
3142 break;
3144 point[1] *= 2.0;
3145 point[2] *= 2.0;
3146 if (child_index % 2 == 1)
3147 point[1] -= 1.0;
3148 if (child_index / 2 == 1)
3149 point[2] -= 1.0;
3150 break;
3152 point *= 2.0;
3153 point -= unit_cell_vertex(child_index);
3154 break;
3155 default:
3157 }
3158
3159 return point;
3160}
3161
3162
3163
3164template <int dim>
3165inline Point<dim>
3167 const Point<dim> & /*p*/,
3168 const unsigned int /*child_index*/,
3169 const RefinementCase<dim> /*refine_case*/)
3170
3171{
3173 return {};
3174}
3175
3176
3177
3178template <>
3179inline Point<1>
3181 const unsigned int child_index,
3182 const RefinementCase<1> refine_case)
3183
3184{
3185 AssertIndexRange(child_index, 2);
3187 (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3188
3189 return (p + unit_cell_vertex(child_index)) * 0.5;
3190}
3191
3192
3193
3194template <>
3195inline Point<3>
3197 const unsigned int child_index,
3198 const RefinementCase<3> refine_case)
3199
3200{
3201 AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3202
3203 Point<3> point = p;
3204 // there might be a cleverer way to do
3205 // this, but since this function is called
3206 // in very few places for initialization
3207 // purposes only, I don't care at the
3208 // moment
3209 switch (refine_case)
3210 {
3212 if (child_index == 1)
3213 point[0] += 1.0;
3214 point[0] *= 0.5;
3215 break;
3217 if (child_index == 1)
3218 point[1] += 1.0;
3219 point[1] *= 0.5;
3220 break;
3222 if (child_index == 1)
3223 point[2] += 1.0;
3224 point[2] *= 0.5;
3225 break;
3227 if (child_index % 2 == 1)
3228 point[0] += 1.0;
3229 if (child_index / 2 == 1)
3230 point[1] += 1.0;
3231 point[0] *= 0.5;
3232 point[1] *= 0.5;
3233 break;
3235 // careful, this is slightly
3236 // different from xy and yz due to
3237 // different internal numbering of
3238 // children!
3239 if (child_index / 2 == 1)
3240 point[0] += 1.0;
3241 if (child_index % 2 == 1)
3242 point[2] += 1.0;
3243 point[0] *= 0.5;
3244 point[2] *= 0.5;
3245 break;
3247 if (child_index % 2 == 1)
3248 point[1] += 1.0;
3249 if (child_index / 2 == 1)
3250 point[2] += 1.0;
3251 point[1] *= 0.5;
3252 point[2] *= 0.5;
3253 break;
3255 point += unit_cell_vertex(child_index);
3256 point *= 0.5;
3257 break;
3258 default:
3260 }
3261
3262 return point;
3263}
3264
3265
3266
3267template <>
3268inline Point<2>
3270 const unsigned int child_index,
3271 const RefinementCase<2> refine_case)
3272{
3273 AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3274
3275 Point<2> point = p;
3276 switch (refine_case)
3277 {
3279 if (child_index == 1)
3280 point[0] += 1.0;
3281 point[0] *= 0.5;
3282 break;
3284 if (child_index == 1)
3285 point[1] += 1.0;
3286 point[1] *= 0.5;
3287 break;
3289 point += unit_cell_vertex(child_index);
3290 point *= 0.5;
3291 break;
3292 default:
3294 }
3295
3296 return point;
3297}
3298
3299
3300
3301template <int dim>
3302inline Point<dim>
3304 const Point<dim> & /*p*/,
3305 const unsigned int /*child_index*/,
3306 const RefinementCase<dim> /*refine_case*/)
3307{
3309 return {};
3310}
3311
3312
3313
3314template <int dim>
3315inline bool
3317{
3319 return false;
3320}
3321
3322template <>
3323inline bool
3325{
3326 return (p[0] >= 0.) && (p[0] <= 1.);
3327}
3328
3329
3330
3331template <>
3332inline bool
3334{
3335 return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
3336}
3337
3338
3339
3340template <>
3341inline bool
3343{
3344 return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
3345 (p[2] >= 0.) && (p[2] <= 1.);
3346}
3347
3348
3349
3350template <int dim>
3351inline bool
3353{
3355 return false;
3356}
3357
3358template <>
3359inline bool
3360GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
3361{
3362 return (p[0] >= -eps) && (p[0] <= 1. + eps);
3363}
3364
3365
3366
3367template <>
3368inline bool
3369GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
3370{
3371 const double l = -eps, u = 1 + eps;
3372 return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
3373}
3374
3375
3376
3377template <>
3378inline bool
3379GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
3380{
3381 const double l = -eps, u = 1.0 + eps;
3382 return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
3383 (p[2] >= l) && (p[2] <= u);
3384}
3385
3386
3387
3388template <>
3389inline unsigned int
3390GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
3391 const unsigned int vertex)
3392{
3393 (void)line;
3394 AssertIndexRange(line, lines_per_cell);
3395 AssertIndexRange(vertex, 2);
3396
3397 return vertex;
3398}
3399
3400
3401template <>
3402inline unsigned int
3403GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
3404 const unsigned int vertex)
3405{
3406 constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
3407 return cell_vertices[line][vertex];
3408}
3409
3410
3411
3412template <>
3413inline unsigned int
3414GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
3415 const unsigned int vertex)
3416{
3417 AssertIndexRange(line, lines_per_cell);
3418 AssertIndexRange(vertex, 2);
3419
3420 constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3421 {1, 3},
3422 {0, 1},
3423 {2, 3},
3424 {4, 6}, // top face
3425 {5, 7},
3426 {4, 5},
3427 {6, 7},
3428 {0,
3429 4}, // connects of bottom
3430 {1, 5}, // top face
3431 {2, 6},
3432 {3, 7}};
3433 return vertices[line][vertex];
3434}
3435
3436
3437template <>
3438inline unsigned int
3439GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3440{
3443}
3444
3445template <>
3446inline unsigned int
3447GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3448 const bool face_orientation,
3449 const bool face_flip,
3450 const bool face_rotation)
3451{
3453
3454 // set up a table to make sure that
3455 // we handle non-standard faces correctly
3456 //
3457 // so set up a table that for each vertex (of
3458 // a quad in standard position) describes
3459 // which vertex to take
3460 //
3461 // first index: four vertices 0...3
3462 //
3463 // second index: face_orientation; 0:
3464 // opposite normal, 1: standard
3465 //
3466 // third index: face_flip; 0: standard, 1:
3467 // face rotated by 180 degrees
3468 //
3469 // forth index: face_rotation: 0: standard,
3470 // 1: face rotated by 90 degrees
3471
3472 constexpr unsigned int vertex_translation[4][2][2][2] = {
3473 {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3474 // face_rotation=false and true
3475 {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3476 // face_rotation=false and true
3477 {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3478 // face_rotation=false and true
3479 {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3480 // face_rotation=false and true
3481
3482 {{{2, 3}, // vertex 1 ...
3483 {1, 0}},
3484 {{1, 0}, {2, 3}}},
3485
3486 {{{1, 0}, // vertex 2 ...
3487 {2, 3}},
3488 {{2, 3}, {1, 0}}},
3489
3490 {{{3, 1}, // vertex 3 ...
3491 {0, 2}},
3492 {{3, 1}, {0, 2}}}};
3493
3494 return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3495}
3496
3497
3498
3499template <int dim>
3500inline unsigned int
3501GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3502 const bool,
3503 const bool,
3504 const bool)
3505{
3506 Assert(dim > 1, ExcImpossibleInDim(dim));
3508 return vertex;
3509}
3510
3511template <int dim>
3512inline unsigned int
3514{
3515 constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3516 0, 2, 2, 4, 2, 4, 4, 8};
3517
3518 return n_children[ref_case];
3519}
3520
3521
3522
3523template <int dim>
3524inline unsigned int
3526{
3528 return 0;
3529}
3530
3531template <>
3532inline unsigned int
3534{
3535 Assert(false, ExcImpossibleInDim(1));
3536 return 0;
3537}
3538
3539template <>
3540inline unsigned int
3542{
3543 return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3544}
3545
3546
3547
3548template <>
3549inline unsigned int
3551{
3552 const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3553 0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3554 return nsubs[subface_case];
3555}
3556
3557
3558
3559template <int dim>
3560inline double
3562 const unsigned int)
3563{
3565 return 0.;
3566}
3567
3568template <>
3569inline double
3571 const unsigned int)
3572{
3573 return 1;
3574}
3575
3576
3577template <>
3578inline double
3580 const unsigned int)
3581{
3582 double ratio = 1;
3583 switch (subface_case)
3584 {
3586 // Here, an
3587 // Assert(false,ExcInternalError())
3588 // would be the right
3589 // choice, but
3590 // unfortunately the
3591 // current function is
3592 // also called for faces
3593 // without children (see
3594 // tests/fe/mapping.cc).
3595 // Assert(false, ExcMessage("Face has no subfaces."));
3596 // Furthermore, assign
3597 // following value as
3598 // otherwise the
3599 // bits/volume_x tests
3600 // break
3602 break;
3604 ratio = 0.5;
3605 break;
3606 default:
3607 // there should be no
3608 // cases left
3610 break;
3611 }
3612
3613 return ratio;
3614}
3615
3616
3617template <>
3618inline double
3620 const unsigned int subface_no)
3621{
3622 double ratio = 1;
3623 switch (subface_case)
3624 {
3626 // Here, an
3627 // Assert(false,ExcInternalError())
3628 // would be the right
3629 // choice, but
3630 // unfortunately the
3631 // current function is
3632 // also called for faces
3633 // without children (see
3634 // tests/bits/mesh_3d_16.cc). Add
3635 // following switch to
3636 // avoid diffs in
3637 // tests/bits/mesh_3d_16
3639 break;
3642 ratio = 0.5;
3643 break;
3647 ratio = 0.25;
3648 break;
3651 if (subface_no < 2)
3652 ratio = 0.25;
3653 else
3654 ratio = 0.5;
3655 break;
3658 if (subface_no == 0)
3659 ratio = 0.5;
3660 else
3661 ratio = 0.25;
3662 break;
3663 default:
3664 // there should be no
3665 // cases left
3667 break;
3668 }
3669
3670 return ratio;
3671}
3672
3673
3674
3675template <int dim>
3677 const RefinementCase<dim> &,
3678 const unsigned int,
3679 const bool,
3680 const bool,
3681 const bool)
3682{
3684 return RefinementCase<dim - 1>::no_refinement;
3685}
3686
3687template <>
3689 const RefinementCase<1> &,
3690 const unsigned int,
3691 const bool,
3692 const bool,
3693 const bool)
3694{
3695 Assert(false, ExcImpossibleInDim(1));
3696
3698}
3699
3700
3701template <>
3702inline RefinementCase<1>
3704 const RefinementCase<2> &cell_refinement_case,
3705 const unsigned int face_no,
3706 const bool,
3707 const bool,
3708 const bool)
3709{
3710 const unsigned int dim = 2;
3711 AssertIndexRange(cell_refinement_case,
3714
3715 // simple special case
3716 if (cell_refinement_case == RefinementCase<dim>::cut_xy)
3718
3719 const RefinementCase<dim - 1>
3722 {RefinementCase<dim - 1>::no_refinement, // no_refinement
3723 RefinementCase<dim - 1>::no_refinement},
3724
3725 {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3726
3727 {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3728
3729 {RefinementCase<dim - 1>::cut_x, // cut_xy
3730 RefinementCase<dim - 1>::cut_x}};
3731
3732 return ref_cases[cell_refinement_case][face_no / 2];
3733}
3734
3735
3736template <>
3737inline RefinementCase<2>
3739 const RefinementCase<3> &cell_refinement_case,
3740 const unsigned int face_no,
3741 const bool face_orientation,
3742 const bool /*face_flip*/,
3743 const bool face_rotation)
3744{
3745 const unsigned int dim = 3;
3746 AssertIndexRange(cell_refinement_case,
3749
3750 // simple special case
3751 if (cell_refinement_case == RefinementCase<dim>::cut_xyz)
3752 return RefinementCase<dim - 1>::cut_xy;
3753
3754 const RefinementCase<dim - 1>
3757 {RefinementCase<dim - 1>::no_refinement, // no_refinement
3758 RefinementCase<dim - 1>::no_refinement,
3759 RefinementCase<dim - 1>::no_refinement},
3760
3761 {RefinementCase<dim - 1>::no_refinement, // cut_x
3762 RefinementCase<dim - 1>::cut_y,
3763 RefinementCase<dim - 1>::cut_x},
3764
3765 {RefinementCase<dim - 1>::cut_x, // cut_y
3766 RefinementCase<dim - 1>::no_refinement,
3767 RefinementCase<dim - 1>::cut_y},
3768
3769 {RefinementCase<dim - 1>::cut_x, // cut_xy
3770 RefinementCase<dim - 1>::cut_y,
3771 RefinementCase<dim - 1>::cut_xy},
3772
3773 {RefinementCase<dim - 1>::cut_y, // cut_z
3774 RefinementCase<dim - 1>::cut_x,
3775 RefinementCase<dim - 1>::no_refinement},
3776
3777 {RefinementCase<dim - 1>::cut_y, // cut_xz
3778 RefinementCase<dim - 1>::cut_xy,
3779 RefinementCase<dim - 1>::cut_x},
3780
3781 {RefinementCase<dim - 1>::cut_xy, // cut_yz
3782 RefinementCase<dim - 1>::cut_x,
3783 RefinementCase<dim - 1>::cut_y},
3784
3785 {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3786 RefinementCase<dim - 1>::cut_xy,
3787 RefinementCase<dim - 1>::cut_xy},
3788 };
3789
3790 const RefinementCase<dim - 1> ref_case =
3791 ref_cases[cell_refinement_case][face_no / 2];
3792
3793 const RefinementCase<dim - 1> flip[4] = {
3794 RefinementCase<dim - 1>::no_refinement,
3795 RefinementCase<dim - 1>::cut_y,
3796 RefinementCase<dim - 1>::cut_x,
3797 RefinementCase<dim - 1>::cut_xy};
3798
3799 // correct the ref_case for face_orientation
3800 // and face_rotation. for face_orientation,
3801 // 'true' is the default value whereas for
3802 // face_rotation, 'false' is standard. If
3803 // <tt>face_rotation==face_orientation</tt>,
3804 // then one of them is non-standard and we
3805 // have to swap cut_x and cut_y, otherwise no
3806 // change is necessary. face_flip has no
3807 // influence. however, in order to keep the
3808 // interface consistent with other functions,
3809 // we still include it as an argument to this
3810 // function
3811 return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3812}
3813
3814
3815
3816template <int dim>
3817inline RefinementCase<1>
3819 const unsigned int)
3820{
3823}
3824
3825template <>
3826inline RefinementCase<1>
3828 const RefinementCase<1> &cell_refinement_case,
3829 const unsigned int line_no)
3830{
3831 (void)line_no;
3832 const unsigned int dim = 1;
3833 (void)dim;
3834 AssertIndexRange(cell_refinement_case,
3837
3838 return cell_refinement_case;
3839}
3840
3841
3842template <>
3843inline RefinementCase<1>
3845 const RefinementCase<2> &cell_refinement_case,
3846 const unsigned int line_no)
3847{
3848 // Assertions are in face_refinement_case()
3849 return face_refinement_case(cell_refinement_case, line_no);
3850}
3851
3852
3853template <>
3854inline RefinementCase<1>
3856 const RefinementCase<3> &cell_refinement_case,
3857 const unsigned int line_no)
3858{
3859 const unsigned int dim = 3;
3860 AssertIndexRange(cell_refinement_case,
3863
3864 // simple special case
3865 if (cell_refinement_case == RefinementCase<dim>::cut_xyz)
3867
3868 // array indicating, which simple refine
3869 // case cuts a line in direction x, y or
3870 // z. For example, cut_y and everything
3871 // containing cut_y (cut_xy, cut_yz,
3872 // cut_xyz) cuts lines, which are in y
3873 // direction.
3874 const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3877
3878 // order the direction of lines
3879 // 0->x, 1->y, 2->z
3880 const unsigned int direction[lines_per_cell] = {
3881 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3882
3883 return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3885 RefinementCase<1>::no_refinement);
3886}
3887
3888
3889
3890template <int dim>
3894 const unsigned int,
3895 const bool,
3896 const bool,
3897 const bool)
3898{
3900
3902}
3903
3904template <>
3905inline RefinementCase<1>
3907 const RefinementCase<0> &,
3908 const unsigned int,
3909 const bool,
3910 const bool,
3911 const bool)
3912{
3913 const unsigned int dim = 1;
3914 Assert(false, ExcImpossibleInDim(dim));
3915
3917}
3918
3919
3920template <>
3921inline RefinementCase<2>
3923 const RefinementCase<1> &face_refinement_case,
3924 const unsigned int face_no,
3925 const bool,
3926 const bool,
3927 const bool)
3928{
3929 const unsigned int dim = 2;
3930 AssertIndexRange(face_refinement_case,
3933
3934 if (face_refinement_case == RefinementCase<dim>::cut_x)
3935 return (face_no / 2) != 0u ? RefinementCase<dim>::cut_x :
3937 else
3939}
3940
3941
3942template <>
3943inline RefinementCase<3>
3945 const RefinementCase<2> &face_refinement_case,
3946 const unsigned int face_no,
3947 const bool face_orientation,
3948 const bool /*face_flip*/,
3949 const bool face_rotation)
3950{
3951 const unsigned int dim = 3;
3952 AssertIndexRange(face_refinement_case,
3955
3960
3961 // correct the face_refinement_case for
3962 // face_orientation and face_rotation. for
3963 // face_orientation, 'true' is the default
3964 // value whereas for face_rotation, 'false'
3965 // is standard. If
3966 // <tt>face_rotation==face_orientation</tt>,
3967 // then one of them is non-standard and we
3968 // have to swap cut_x and cut_y, otherwise no
3969 // change is necessary. face_flip has no
3970 // influence. however, in order to keep the
3971 // interface consistent with other functions,
3972 // we still include it as an argument to this
3973 // function
3974 const RefinementCase<dim - 1> std_face_ref =
3975 (face_orientation == face_rotation) ? flip[face_refinement_case] :
3976 face_refinement_case;
3977
3978 const RefinementCase<dim> face_to_cell[3][4] = {
3979 {RefinementCase<dim>::no_refinement, // faces 0 and 1
3980 RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3983
3984 {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3985 // "exchanged on faces 2 and 3")
3989
3990 {RefinementCase<dim>::no_refinement, // faces 4 and 5
3994
3995 return face_to_cell[face_no / 2][std_face_ref];
3996}
3997
3998
3999
4000template <int dim>
4003 const unsigned int)
4004{
4006
4008}
4009
4010template <>
4011inline RefinementCase<1>
4013 const unsigned int line_no)
4014{
4015 (void)line_no;
4016 AssertIndexRange(line_no, 1);
4017
4019}
4020
4021
4022template <>
4023inline RefinementCase<2>
4025 const unsigned int line_no)
4026{
4027 const unsigned int dim = 2;
4028 (void)dim;
4030
4031 return (line_no / 2) != 0u ? RefinementCase<2>::cut_x :
4033}
4034
4035
4036template <>
4037inline RefinementCase<3>
4039 const unsigned int line_no)
4040{
4041 const unsigned int dim = 3;
4043
4044 const RefinementCase<dim> ref_cases[6] = {
4045 RefinementCase<dim>::cut_y, // lines 0 and 1
4046 RefinementCase<dim>::cut_x, // lines 2 and 3
4047 RefinementCase<dim>::cut_y, // lines 4 and 5
4048 RefinementCase<dim>::cut_x, // lines 6 and 7
4049 RefinementCase<dim>::cut_z, // lines 8 and 9
4050 RefinementCase<dim>::cut_z}; // lines 10 and 11
4051
4052 return ref_cases[line_no / 2];
4053}
4054
4055
4056
4057template <>
4058inline unsigned int
4059GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
4060 const bool face_orientation,
4061 const bool face_flip,
4062 const bool face_rotation)
4063{
4065
4066 // set up a table to make sure that
4067 // we handle non-standard faces correctly
4068 //
4069 // so set up a table that for each vertex (of
4070 // a quad in standard position) describes
4071 // which vertex to take
4072 //
4073 // first index: four vertices 0...3
4074 //
4075 // second index: face_orientation; 0:
4076 // opposite normal, 1: standard
4077 //
4078 // third index: face_flip; 0: standard, 1:
4079 // face rotated by 180 degrees
4080 //
4081 // forth index: face_rotation: 0: standard,
4082 // 1: face rotated by 90 degrees
4083
4084 const unsigned int vertex_translation[4][2][2][2] = {
4085 {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
4086 // face_rotation=false and true
4087 {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
4088 // face_rotation=false and true
4089 {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
4090 // face_rotation=false and true
4091 {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
4092 // face_rotation=false and true
4093
4094 {{{2, 3}, // vertex 1 ...
4095 {1, 0}},
4096 {{1, 3}, {2, 0}}},
4097
4098 {{{1, 0}, // vertex 2 ...
4099 {2, 3}},
4100 {{2, 0}, {1, 3}}},
4101
4102 {{{3, 1}, // vertex 3 ...
4103 {0, 2}},
4104 {{3, 2}, {0, 1}}}};
4105
4106 return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
4107}
4108
4109
4110
4111template <int dim>
4112inline unsigned int
4113GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
4114 const bool,
4115 const bool,
4116 const bool)
4117{
4118 Assert(dim > 1, ExcImpossibleInDim(dim));
4120 return vertex;
4121}
4122
4123
4124
4125template <>
4126inline unsigned int
4127GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
4128 const bool face_orientation,
4129 const bool face_flip,
4130 const bool face_rotation)
4131{
4133
4134
4135 // make sure we handle
4136 // non-standard faces correctly
4137 //
4138 // so set up a table that for each line (of a
4139 // quad) describes which line to take
4140 //
4141 // first index: four lines 0...3
4142 //
4143 // second index: face_orientation; 0:
4144 // opposite normal, 1: standard
4145 //
4146 // third index: face_flip; 0: standard, 1:
4147 // face rotated by 180 degrees
4148 //
4149 // forth index: face_rotation: 0: standard,
4150 // 1: face rotated by 90 degrees
4151
4152 const unsigned int line_translation[4][2][2][2] = {
4153 {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4154 // face_rotation=false and true
4155 {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4156 // face_rotation=false and true
4157 {{0, 3}, // line 0, face_orientation=true, face_flip=false,
4158 // face_rotation=false and true
4159 {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
4160 // face_rotation=false and true
4161
4162 {{{3, 1}, // line 1 ...
4163 {2, 0}},
4164 {{1, 2}, {0, 3}}},
4165
4166 {{{0, 3}, // line 2 ...
4167 {1, 2}},
4168 {{2, 0}, {3, 1}}},
4169
4170 {{{1, 2}, // line 3 ...
4171 {0, 3}},
4172 {{3, 1}, {2, 0}}}};
4173
4174 return line_translation[line][face_orientation][face_flip][face_rotation];
4175}
4176
4177
4178
4179template <int dim>
4180inline unsigned int
4181GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
4182 const bool,
4183 const bool,
4184 const bool)
4185{
4187 return line;
4188}
4189
4190
4191
4192template <>
4193inline unsigned int
4194GeometryInfo<2>::standard_to_real_line_vertex(const unsigned int vertex,
4195 const bool line_orientation)
4196{
4197 return line_orientation ? vertex : (1 - vertex);
4198}
4199
4200
4201
4202template <int dim>
4203inline unsigned int
4204GeometryInfo<dim>::standard_to_real_line_vertex(const unsigned int vertex,
4205 const bool)
4206{
4208 return vertex;
4209}
4210
4211
4212
4213template <>
4214inline std::array<unsigned int, 2>
4216 const unsigned int vertex)
4217{
4218 return {{vertex % 2, vertex / 2}};
4219}
4220
4221
4222
4223template <int dim>
4224inline std::array<unsigned int, 2>
4226 const unsigned int vertex)
4227{
4229 (void)vertex;
4230 return {{0, 0}};
4231}
4232
4233
4234
4235template <>
4236inline std::array<unsigned int, 2>
4238{
4239 // set up a table that for each
4240 // line describes a) from which
4241 // quad to take it, b) which line
4242 // therein it is if the face is
4243 // oriented correctly
4244 static const unsigned int lookup_table[GeometryInfo<3>::lines_per_cell][2] = {
4245 {4, 0}, // take first four lines from bottom face
4246 {4, 1},
4247 {4, 2},
4248 {4, 3},
4249
4250 {5, 0}, // second four lines from top face
4251 {5, 1},
4252 {5, 2},
4253 {5, 3},
4254
4255 {0, 0}, // the rest randomly
4256 {1, 0},
4257 {0, 1},
4258 {1, 1}};
4259
4260 return {{lookup_table[i][0], lookup_table[i][1]}};
4261}
4262
4263
4264
4265template <int dim>
4266inline std::array<unsigned int, 2>
4268{
4270 (void)line;
4271 return {{0, 0}};
4272}
4273
4274
4275
4276template <>
4277inline std::array<unsigned int, 2>
4279 const unsigned int vertex)
4280{
4281 // get the corner indices by asking either the bottom or the top face for its
4282 // vertices. handle non-standard faces by calling the vertex reordering
4283 // function from GeometryInfo
4284
4285 // bottom face (4) for first four vertices, top face (5) for the rest
4286 return {{4 + vertex / 4, vertex % 4}};
4287}
4288
4289
4290
4291template <int dim>
4292inline std::array<unsigned int, 2>
4294 const unsigned int vertex)
4295{
4297 (void)vertex;
4298 return {{0, 0}};
4299}
4300
4301
4302
4303template <>
4304inline unsigned int
4305GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
4306 const bool face_orientation,
4307 const bool face_flip,
4308 const bool face_rotation)
4309{
4311
4312
4313 // make sure we handle
4314 // non-standard faces correctly
4315 //
4316 // so set up a table that for each line (of a
4317 // quad) describes which line to take
4318 //
4319 // first index: four lines 0...3
4320 //
4321 // second index: face_orientation; 0:
4322 // opposite normal, 1: standard
4323 //
4324 // third index: face_flip; 0: standard, 1:
4325 // face rotated by 180 degrees
4326 //
4327 // forth index: face_rotation: 0: standard,
4328 // 1: face rotated by 90 degrees
4329
4330 const unsigned int line_translation[4][2][2][2] = {
4331 {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4332 // face_rotation=false and true
4333 {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4334 // face_rotation=false and true
4335 {{0, 2}, // line 0, face_orientation=true, face_flip=false,
4336 // face_rotation=false and true
4337 {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
4338 // face_rotation=false and true
4339
4340 {{{3, 1}, // line 1 ...
4341 {2, 0}},
4342 {{1, 3}, {0, 2}}},
4343
4344 {{{0, 3}, // line 2 ...
4345 {1, 2}},
4346 {{2, 1}, {3, 0}}},
4347
4348 {{{1, 2}, // line 3 ...
4349 {0, 3}},
4350 {{3, 0}, {2, 1}}}};
4351
4352 return line_translation[line][face_orientation][face_flip][face_rotation];
4353}
4354
4355
4356
4357template <int dim>
4358inline unsigned int
4359GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
4360 const bool,
4361 const bool,
4362 const bool)
4363{
4365 return line;
4366}
4367
4368
4369
4370template <>
4371inline unsigned int
4373 const unsigned int face,
4374 const unsigned int subface,
4375 const bool,
4376 const bool,
4377 const bool,
4378 const RefinementCase<0> &)
4379{
4380 (void)subface;
4381 AssertIndexRange(face, faces_per_cell);
4382 AssertIndexRange(subface, max_children_per_face);
4383
4384 return face;
4385}
4386
4387
4388
4389template <>
4390inline unsigned int
4392 const unsigned int face,
4393 const unsigned int subface,
4394 const bool face_orientation,
4395 const bool /*face_flip*/,
4396 const bool /*face_rotation*/,
4397 const RefinementCase<1> &)
4398{
4399 AssertIndexRange(face, faces_per_cell);
4400 AssertIndexRange(subface, max_children_per_face);
4401
4402 // always return the child adjacent to the specified
4403 // subface. if the face of a cell is not refined, don't
4404 // throw an assertion but deliver the child adjacent to
4405 // the face nevertheless, i.e. deliver the child of
4406 // this cell adjacent to the subface of a possibly
4407 // refined neighbor. this simplifies setting neighbor
4408 // information in execute_refinement.
4409 constexpr unsigned int
4410 subcells[/* possible face orientation */ 2]
4411 [/* number of different ways to refine a cell */ 4]
4412 [/* faces_per_cell */ 4][/* max_children_per_face */ 2] = {
4413 {
4414 // Flipped orientation (face_orientation = false)
4415 {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
4416 {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
4417 {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
4418 },
4419 {
4420 // Normal orientation (face_orientation = true)
4421 {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
4422 {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
4423 {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
4424 }};
4425
4426 return subcells[face_orientation][ref_case - 1][face][subface];
4427}
4428
4429
4430
4431template <>
4432inline unsigned int
4434 const unsigned int face,
4435 const unsigned int subface,
4436 const bool face_orientation,
4437 const bool face_flip,
4438 const bool face_rotation,
4439 const RefinementCase<2> &face_ref_case)
4440{
4441 const unsigned int dim = 3;
4442
4444 ExcMessage("Cell has no children."));
4445 AssertIndexRange(face, faces_per_cell);
4446 if (!(subface == 0 &&
4448 {
4449 AssertIndexRange(subface,
4450 GeometryInfo<dim - 1>::n_children(face_ref_case));
4451 }
4452
4453 // invalid number used for invalid cases,
4454 // e.g. when the children are more refined at
4455 // a given face than the face itself
4456 const unsigned int e = numbers::invalid_unsigned_int;
4457
4458 // the whole process of finding a child cell
4459 // at a given subface considering the
4460 // possibly anisotropic refinement cases of
4461 // the cell and the face as well as
4462 // orientation, flip and rotation of the face
4463 // is quite complicated. thus, we break it
4464 // down into several steps.
4465
4466 // first step: convert the given face refine
4467 // case to a face refine case concerning the
4468 // face in standard orientation (, flip and
4469 // rotation). This only affects cut_x and
4470 // cut_y
4471 const RefinementCase<dim - 1> flip[4] = {
4472 RefinementCase<dim - 1>::no_refinement,
4473 RefinementCase<dim - 1>::cut_y,
4474 RefinementCase<dim - 1>::cut_x,
4475 RefinementCase<dim - 1>::cut_xy};
4476 // for face_orientation, 'true' is the
4477 // default value whereas for face_rotation,
4478 // 'false' is standard. If
4479 // <tt>face_rotation==face_orientation</tt>,
4480 // then one of them is non-standard and we
4481 // have to swap cut_x and cut_y, otherwise no
4482 // change is necessary.
4483 const RefinementCase<dim - 1> std_face_ref =
4484 (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4485
4486 // second step: convert the given subface
4487 // index to the one for a standard face
4488 // respecting face_orientation, face_flip and
4489 // face_rotation
4490
4491 // first index: face_ref_case
4492 // second index: face_orientation
4493 // third index: face_flip
4494 // forth index: face_rotation
4495 // fifth index: subface index
4496 const unsigned int subface_exchange[4][2][2][2][4] = {
4497 // no_refinement (subface 0 stays 0,
4498 // all others are invalid)
4499 {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4500 {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4501 // cut_x (here, if the face is only
4502 // rotated OR only falsely oriented,
4503 // then subface 0 of the non-standard
4504 // face does NOT correspond to one of
4505 // the subfaces of a standard
4506 // face. Thus we indicate the subface
4507 // which is located at the lower left
4508 // corner (the origin of the face's
4509 // local coordinate system) with
4510 // '0'. The rest of this issue is
4511 // taken care of using the above
4512 // conversion to a 'standard face
4513 // refine case')
4514 {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4515 {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4516 // cut_y (the same applies as for
4517 // cut_x)
4518 {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4519 {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4520 // cut_xyz: this information is
4521 // identical to the information
4522 // returned by
4523 // GeometryInfo<3>::real_to_standard_face_vertex()
4524 {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4525 // face_rotation=false, subfaces 0,1,2,3
4526 {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4527 // face_rotation=true, subfaces 0,1,2,3
4528 {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4529 // face_rotation=false, subfaces 0,1,2,3
4530 {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4531 // face_rotation=true, subfaces 0,1,2,3
4532 {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4533 // face_rotation=false, subfaces 0,1,2,3
4534 {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4535 // face_rotation=true, subfaces 0,1,2,3
4536 {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4537 // face_rotation=false, subfaces 0,1,2,3
4538 {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4539 // face_rotation=true, subfaces 0,1,2,3
4540
4541 const unsigned int std_subface =
4542 subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4543 [subface];
4544 Assert(std_subface != e, ExcInternalError());
4545
4546 // third step: these are the children, which
4547 // can be found at the given subfaces of an
4548 // isotropically refined (standard) face
4549 //
4550 // first index: (refinement_case-1)
4551 // second index: face_index
4552 // third index: subface_index (isotropic refinement)
4553 const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4554 [max_children_per_face] = {
4555 // cut_x
4556 {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4557 {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4558 {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4559 {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4560 {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4561 {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4562 // cut_y
4563 {{0, 1, 0, 1},
4564 {0, 1, 0, 1},
4565 {0, 0, 0, 0},
4566 {1, 1, 1, 1},
4567 {0, 0, 1, 1},
4568 {0, 0, 1, 1}},
4569 // cut_xy
4570 {{0, 2, 0, 2},
4571 {1, 3, 1, 3},
4572 {0, 0, 1, 1},
4573 {2, 2, 3, 3},
4574 {0, 1, 2, 3},
4575 {0, 1, 2, 3}},
4576 // cut_z
4577 {{0, 0, 1, 1},
4578 {0, 0, 1, 1},
4579 {0, 1, 0, 1},
4580 {0, 1, 0, 1},
4581 {0, 0, 0, 0},
4582 {1, 1, 1, 1}},
4583 // cut_xz
4584 {{0, 0, 1, 1},
4585 {2, 2, 3, 3},
4586 {0, 1, 2, 3},
4587 {0, 1, 2, 3},
4588 {0, 2, 0, 2},
4589 {1, 3, 1, 3}},
4590 // cut_yz
4591 {{0, 1, 2, 3},
4592 {0, 1, 2, 3},
4593 {0, 2, 0, 2},
4594 {1, 3, 1, 3},
4595 {0, 0, 1, 1},
4596 {2, 2, 3, 3}},
4597 // cut_xyz
4598 {{0, 2, 4, 6},
4599 {1, 3, 5, 7},
4600 {0, 4, 1, 5},
4601 {2, 6, 3, 7},
4602 {0, 1, 2, 3},
4603 {4, 5, 6, 7}}};
4604
4605 // forth step: check, whether the given face
4606 // refine case is valid for the given cell
4607 // refine case. this is the case, if the
4608 // given face refine case is at least as
4609 // refined as the face is for the given cell
4610 // refine case
4611
4612 // note, that we are considering standard
4613 // face refinement cases here and thus must
4614 // not pass the given orientation, flip and
4615 // rotation flags
4616 if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4617 face_refinement_case(ref_case, face))
4618 {
4619 // all is fine. for anisotropic face
4620 // refine cases, select one of the
4621 // isotropic subfaces which neighbors the
4622 // same child
4623
4624 // first index: (standard) face refine case
4625 // second index: subface index
4626 const unsigned int equivalent_iso_subface[4][4] = {
4627 {0, e, e, e}, // no_refinement
4628 {0, 3, e, e}, // cut_x
4629 {0, 3, e, e}, // cut_y
4630 {0, 1, 2, 3}}; // cut_xy
4631
4632 const unsigned int equ_std_subface =
4633 equivalent_iso_subface[std_face_ref][std_subface];
4634 Assert(equ_std_subface != e, ExcInternalError());
4635
4636 return iso_children[ref_case - 1][face][equ_std_subface];
4637 }
4638 else
4639 {
4640 // the face_ref_case was too coarse,
4641 // throw an error
4642 Assert(false,
4643 ExcMessage("The face RefineCase is too coarse "
4644 "for the given cell RefineCase."));
4645 }
4646 // we only get here in case of an error
4647 return e;
4648}
4649
4650
4651
4652template <>
4653inline unsigned int
4655 const unsigned int,
4656 const unsigned int,
4657 const bool,
4658 const bool,
4659 const bool,
4660 const RefinementCase<3> &)
4661{
4664}
4665
4666
4667
4668template <>
4669inline unsigned int
4670GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4671 const unsigned int line,
4672 const bool,
4673 const bool,
4674 const bool)
4675{
4676 (void)line;
4677 AssertIndexRange(face, faces_per_cell);
4678 AssertIndexRange(line, lines_per_face);
4679
4680 // The face is a line itself.
4681 return face;
4682}
4683
4684
4685
4686template <>
4687inline unsigned int
4688GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4689 const unsigned int line,
4690 const bool face_orientation,
4691 const bool face_flip,
4692 const bool face_rotation)
4693{
4694 AssertIndexRange(face, faces_per_cell);
4695 AssertIndexRange(line, lines_per_face);
4696
4697 const unsigned lines[faces_per_cell][lines_per_face] = {
4698 {8, 10, 0, 4}, // left face
4699 {9, 11, 1, 5}, // right face
4700 {2, 6, 8, 9}, // front face
4701 {3, 7, 10, 11}, // back face
4702 {0, 1, 2, 3}, // bottom face
4703 {4, 5, 6, 7}}; // top face
4704 return lines[face][real_to_standard_face_line(
4705 line, face_orientation, face_flip, face_rotation)];
4706}
4707
4708
4709
4710inline unsigned int
4711GeometryInfo<0>::face_to_cell_lines(const unsigned int,
4712 const unsigned int,
4713 const bool,
4714 const bool,
4715 const bool)
4716{
4719}
4720
4721
4722
4723template <int dim>
4724inline unsigned int
4725GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4726 const unsigned int,
4727 const bool,
4728 const bool,
4729 const bool)
4730{
4733}
4734
4735
4736
4737template <int dim>
4738inline unsigned int
4739GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4740 const unsigned int vertex,
4741 const bool face_orientation,
4742 const bool face_flip,
4743 const bool face_rotation)
4744{
4745 return child_cell_on_face(RefinementCase<dim>::isotropic_refinement,
4746 face,
4747 vertex,
4748 face_orientation,
4749 face_flip,
4750 face_rotation);
4751}
4752
4753
4754
4755inline unsigned int
4757 const unsigned int,
4758 const bool,
4759 const bool,
4760 const bool)
4761{
4764}
4765
4766
4767
4768template <int dim>
4769template <typename Number>
4770inline Point<dim, Number>
4772{
4774 for (unsigned int i = 0; i < dim; ++i)
4775 p[i] = std::min(std::max(q[i], Number(0.)), Number(1.));
4776
4777 return p;
4778}
4779
4780
4781
4782template <int dim>
4783inline double
4785{
4786 double result = 0.0;
4787
4788 for (unsigned int i = 0; i < dim; ++i)
4789 {
4790 result = std::max(result, -p[i]);
4791 result = std::max(result, p[i] - 1.);
4792 }
4793
4794 return result;
4795}
4796
4797
4798
4799template <int dim>
4800inline double
4802 const unsigned int i)
4803{
4805
4806 switch (dim)
4807 {
4808 case 1:
4809 {
4810 const double x = xi[0];
4811 switch (i)
4812 {
4813 case 0:
4814 return 1 - x;
4815 case 1:
4816 return x;
4817 }
4818 break;
4819 }
4820
4821 case 2:
4822 {
4823 const double x = xi[0];
4824 const double y = xi[1];
4825 switch (i)
4826 {
4827 case 0:
4828 return (1 - x) * (1 - y);
4829 case 1:
4830 return x * (1 - y);
4831 case 2:
4832 return (1 - x) * y;
4833 case 3:
4834 return x * y;
4835 }
4836 break;
4837 }
4838
4839 case 3:
4840 {
4841 const double x = xi[0];
4842 const double y = xi[1];
4843 const double z = xi[2];
4844 switch (i)
4845 {
4846 case 0:
4847 return (1 - x) * (1 - y) * (1 - z);
4848 case 1:
4849 return x * (1 - y) * (1 - z);
4850 case 2:
4851 return (1 - x) * y * (1 - z);
4852 case 3:
4853 return x * y * (1 - z);
4854 case 4:
4855 return (1 - x) * (1 - y) * z;
4856 case 5:
4857 return x * (1 - y) * z;
4858 case 6:
4859 return (1 - x) * y * z;
4860 case 7:
4861 return x * y * z;
4862 }
4863 break;
4864 }
4865
4866 default:
4868 }
4869 return -1e9;
4870}
4871
4872
4873
4874template <>
4876 const Point<1> &,
4877 const unsigned int i)
4878{
4880
4881 switch (i)
4882 {
4883 case 0:
4884 return Point<1>(-1.);
4885 case 1:
4886 return Point<1>(1.);
4887 }
4888
4889 return Point<1>(-1e9);
4890}
4891
4892
4893
4894template <>
4896 const Point<2> &xi,
4897 const unsigned int i)
4898{
4900
4901 const double x = xi[0];
4902 const double y = xi[1];
4903 switch (i)
4904 {
4905 case 0:
4906 return Point<2>(-(1 - y), -(1 - x));
4907 case 1:
4908 return Point<2>(1 - y, -x);
4909 case 2:
4910 return Point<2>(-y, 1 - x);
4911 case 3:
4912 return Point<2>(y, x);
4913 }
4914 return Point<2>(-1e9, -1e9);
4915}
4916
4917
4918
4919template <>
4921 const Point<3> &xi,
4922 const unsigned int i)
4923{
4925
4926 const double x = xi[0];
4927 const double y = xi[1];
4928 const double z = xi[2];
4929 switch (i)
4930 {
4931 case 0:
4932 return Point<3>(-(1 - y) * (1 - z),
4933 -(1 - x) * (1 - z),
4934 -(1 - x) * (1 - y));
4935 case 1:
4936 return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4937 case 2:
4938 return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4939 case 3:
4940 return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4941 case 4:
4942 return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4943 case 5:
4944 return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4945 case 6:
4946 return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4947 case 7:
4948 return Point<3>(y * z, x * z, x * y);
4949 }
4950
4951 return Point<3>(-1e9, -1e9, -1e9);
4952}
4953
4954
4955
4956template <int dim>
4957inline Tensor<1, dim>
4959 const unsigned int)
4960{
4962 return Tensor<1, dim>();
4963}
4964
4965
4966unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4967{
4968 return 0;
4969}
4970
4971
4972namespace internal
4973{
4974 namespace GeometryInfoHelper
4975 {
4976 // wedge product of a single
4977 // vector in 2d: we just have to
4978 // rotate it by 90 degrees to the
4979 // right
4980 inline Tensor<1, 2>
4981 wedge_product(const Tensor<1, 2> (&derivative)[1])
4982 {
4983 Tensor<1, 2> result;
4984 result[0] = derivative[0][1];
4985 result[1] = -derivative[0][0];
4986
4987 return result;
4988 }
4989
4990
4991 // wedge product of 2 vectors in
4992 // 3d is the cross product
4993 inline Tensor<1, 3>
4994 wedge_product(const Tensor<1, 3> (&derivative)[2])
4995 {
4996 return cross_product_3d(derivative[0], derivative[1]);
4997 }
4998
4999
5000 // wedge product of dim vectors
5001 // in dim-d: that's the
5002 // determinant of the matrix
5003 template <int dim>
5004 inline Tensor<0, dim>
5005 wedge_product(const Tensor<1, dim> (&derivative)[dim])
5006 {
5007 Tensor<2, dim> jacobian;
5008 for (unsigned int i = 0; i < dim; ++i)
5009 jacobian[i] = derivative[i];
5010
5011 return determinant(jacobian);
5012 }
5013 } // namespace GeometryInfoHelper
5014} // namespace internal
5015
5016
5017template <int dim>
5018template <int spacedim>
5019inline void
5021# ifndef DEAL_II_CXX14_CONSTEXPR_BUG
5022 (const Point<spacedim> (&vertices)[vertices_per_cell],
5023 Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell])
5024# else
5025 (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms)
5026# endif
5027{
5028 // for each of the vertices,
5029 // compute the alternating form
5030 // of the mapped unit
5031 // vectors. consider for
5032 // example the case of a quad
5033 // in spacedim==3: to do so, we
5034 // need to see how the
5035 // infinitesimal vectors
5036 // (d\xi_1,0) and (0,d\xi_2)
5037 // are transformed into
5038 // spacedim-dimensional space
5039 // and then form their cross
5040 // product (i.e. the wedge product
5041 // of two vectors). to this end, note
5042 // that
5043 // \vec x = sum_i \vec v_i phi_i(\vec xi)
5044 // so the transformed vectors are
5045 // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
5046 // and
5047 // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
5048 // which boils down to the columns
5049 // of the 3x2 matrix \grad_\xi x(\xi)
5050 //
5051 // a similar reasoning would
5052 // hold for all dim,spacedim
5053 // pairs -- we only have to
5054 // compute the wedge product of
5055 // the columns of the
5056 // derivatives
5057 for (unsigned int i = 0; i < vertices_per_cell; ++i)
5058 {
5059 Tensor<1, spacedim> derivatives[dim];
5060
5061 for (unsigned int j = 0; j < vertices_per_cell; ++j)
5062 {
5063 const Tensor<1, dim> grad_phi_j =
5064 d_linear_shape_function_gradient(unit_cell_vertex(i), j);
5065 for (unsigned int l = 0; l < dim; ++l)
5066 derivatives[l] += vertices[j] * grad_phi_j[l];
5067 }
5068
5069 forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
5070 }
5071}
5072
5073#endif // DOXYGEN
5075
5076#endif
GeometryPrimitive(const unsigned int object_dimension)
GeometryPrimitive(const Object object)
Definition point.h:113
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:40
#define DEAL_II_HOST_DEVICE
Definition config.h:169
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static ::ExceptionBase & ExcInvalidCoordinate(double arg1)
static ::ExceptionBase & ExcInvalidSubface(int arg1, int arg2, int arg3)
static ::ExceptionBase & ExcInvalidRefinementCase(int arg1)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
#define DeclException3(Exception3, type1, type2, type3, outsequence)
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcInvalidSubfaceCase(int arg1)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:193
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)
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
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
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 > &)