deal.II version GIT relicensing-3535-g02565b63bb 2025-06-19 13:50:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
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
1963template <int dim>
1965{
1973 static constexpr unsigned int max_children_per_cell = 1 << dim;
1974
1978 static constexpr unsigned int faces_per_cell = 2 * dim;
1979
1998
2006 static constexpr unsigned int max_children_per_face =
2008
2012 static constexpr unsigned int vertices_per_cell = 1 << dim;
2013
2031
2035 static constexpr unsigned int vertices_per_face =
2037
2041 static constexpr unsigned int lines_per_face =
2043
2047 static constexpr unsigned int quads_per_face =
2049
2059 static constexpr unsigned int lines_per_cell =
2060 (2 * GeometryInfo<dim - 1>::lines_per_cell +
2062
2070 static constexpr unsigned int quads_per_cell =
2071 (2 * GeometryInfo<dim - 1>::quads_per_cell +
2073
2077 static constexpr unsigned int hexes_per_cell =
2078 (2 * GeometryInfo<dim - 1>::hexes_per_cell +
2080
2098 static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
2099 internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal();
2100
2114 static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
2115 internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal();
2116
2129 internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face();
2130
2135 static unsigned int
2136 n_children(const RefinementCase<dim> &refinement_case);
2137
2142 static unsigned int
2144
2154 static double
2156 const unsigned int subface_no);
2157
2163 static RefinementCase<dim - 1>
2164 face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2165 const unsigned int face_no,
2166 const bool face_orientation = true,
2167 const bool face_flip = false,
2168 const bool face_rotation = false);
2169
2175 static RefinementCase<dim>
2178 const unsigned int face_no,
2179 const bool face_orientation = true,
2180 const bool face_flip = false,
2181 const bool face_rotation = false);
2182
2187 static RefinementCase<1>
2188 line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2189 const unsigned int line_no);
2190
2195 static RefinementCase<dim>
2197
2244 static unsigned int
2246 const unsigned int face,
2247 const unsigned int subface,
2248 const bool face_orientation = true,
2249 const bool face_flip = false,
2250 const bool face_rotation = false,
2253
2267 static unsigned int
2268 line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
2269
2290 static unsigned int
2291 face_to_cell_vertices(const unsigned int face,
2292 const unsigned int vertex,
2293 const bool face_orientation = true,
2294 const bool face_flip = false,
2295 const bool face_rotation = false);
2296
2308 static unsigned int
2309 face_to_cell_lines(const unsigned int face,
2310 const unsigned int line,
2311 const bool face_orientation = true,
2312 const bool face_flip = false,
2313 const bool face_rotation = false);
2314
2324 static unsigned int
2325 standard_to_real_face_vertex(const unsigned int vertex,
2326 const bool face_orientation = true,
2327 const bool face_flip = false,
2328 const bool face_rotation = false);
2329
2339 static unsigned int
2340 real_to_standard_face_vertex(const unsigned int vertex,
2341 const bool face_orientation = true,
2342 const bool face_flip = false,
2343 const bool face_rotation = false);
2344
2354 static unsigned int
2355 standard_to_real_face_line(const unsigned int line,
2356 const bool face_orientation = true,
2357 const bool face_flip = false,
2358 const bool face_rotation = false);
2359
2365 static unsigned int
2366 standard_to_real_line_vertex(const unsigned int vertex,
2367 const bool line_orientation = true);
2368
2376 static std::array<unsigned int, 2>
2378
2386 static std::array<unsigned int, 2>
2388
2396 static std::array<unsigned int, 2>
2398
2408 static unsigned int
2409 real_to_standard_face_line(const unsigned int line,
2410 const bool face_orientation = true,
2411 const bool face_flip = false,
2412 const bool face_rotation = false);
2413
2419 static Point<dim>
2420 unit_cell_vertex(const unsigned int vertex);
2421
2431 static unsigned int
2433
2441 static Point<dim>
2443 const unsigned int child_index,
2444 const RefinementCase<dim> refine_case =
2446
2452 static Point<dim>
2454 const unsigned int child_index,
2455 const RefinementCase<dim> refine_case =
2457
2462 static bool
2464
2476 static bool
2477 is_inside_unit_cell(const Point<dim> &p, const double eps);
2478
2483 template <typename Number = double>
2484 static Point<dim, Number>
2486
2492 static double
2494
2499 static double
2500 d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2501
2506 static Tensor<1, dim>
2507 d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2508
2560 template <int spacedim>
2561 static void
2563#ifndef DEAL_II_CXX14_CONSTEXPR_BUG
2564 (const Point<spacedim> (&vertices)[vertices_per_cell],
2565 Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2566#else
2567 (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2568#endif
2569
2579 static constexpr std::array<unsigned int, faces_per_cell>
2581 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction();
2582
2599 static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2600 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation();
2601
2612 static constexpr std::array<Tensor<1, dim>, faces_per_cell>
2614 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_vector();
2615
2629 static constexpr ndarray<Tensor<1, dim>, faces_per_cell, dim - 1>
2630 unit_tangential_vectors = internal::GeometryInfoHelper::Initializers<
2632
2638 static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2639 internal::GeometryInfoHelper::Initializers<dim>::opposite_face();
2640
2641
2646 double,
2647 << "The coordinates must satisfy 0 <= x_i <= 1, "
2648 << "but here we have x_i=" << arg1);
2649
2654 int,
2655 int,
2656 int,
2657 << "RefinementCase<dim> " << arg1 << ": face " << arg2
2658 << " has no subface " << arg3);
2659};
2660
2661
2662
2663#ifndef DOXYGEN
2664
2665
2666/* -------------- declaration of explicit specializations ------------- */
2667
2668
2669template <>
2672 const unsigned int i);
2673template <>
2676 const unsigned int i);
2677template <>
2680 const unsigned int i);
2681
2682
2683
2684/* -------------- inline functions ------------- */
2685
2686
2687inline GeometryPrimitive::GeometryPrimitive(const Object object)
2688 : object(object)
2689{}
2690
2691
2692
2693inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2694 : object(static_cast<Object>(object_dimension))
2695{}
2696
2697
2698inline GeometryPrimitive::operator unsigned int() const
2699{
2700 return static_cast<unsigned int>(object);
2701}
2702
2703
2704
2705namespace internal
2706{
2707 template <int dim>
2709 const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2710 : value(subface_possibility)
2711 {}
2712
2713
2714 template <int dim>
2715 inline SubfaceCase<dim>::operator std::uint8_t() const
2716 {
2717 return value;
2718 }
2719
2720
2721} // namespace internal
2722
2723
2724template <int dim>
2726RefinementCase<dim>::cut_axis(const unsigned int)
2727{
2729 return static_cast<std::uint8_t>(-1);
2730}
2731
2732
2733template <>
2734inline RefinementCase<1>
2735RefinementCase<1>::cut_axis(const unsigned int i)
2736{
2737 AssertIndexRange(i, 1);
2738
2740 return options[i];
2741}
2742
2743
2744
2745template <>
2746inline RefinementCase<2>
2747RefinementCase<2>::cut_axis(const unsigned int i)
2748{
2749 AssertIndexRange(i, 2);
2750
2753 return options[i];
2754}
2755
2756
2757
2758template <>
2759inline RefinementCase<3>
2760RefinementCase<3>::cut_axis(const unsigned int i)
2761{
2762 AssertIndexRange(i, 3);
2763
2767 return options[i];
2768}
2769
2770
2771
2772template <>
2773inline std::array<RefinementCase<1>, 2>
2775{
2778}
2779
2780
2781
2782template <>
2783inline std::array<RefinementCase<2>, 4>
2785{
2790}
2791
2792
2793
2794template <>
2795inline std::array<RefinementCase<3>, 8>
2797{
2806}
2807
2808
2809
2810template <int dim>
2812 : value(RefinementPossibilities<dim>::no_refinement)
2813{}
2814
2815
2816
2817template <int dim>
2819 const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2820 : value(refinement_case)
2821{
2822 // check that only those bits of
2823 // the given argument are set that
2824 // make sense for a given space
2825 // dimension
2826 Assert((refinement_case &
2828 refinement_case,
2829 ExcInvalidRefinementCase(refinement_case));
2830}
2831
2832
2833
2834template <int dim>
2835inline RefinementCase<dim>::RefinementCase(const std::uint8_t refinement_case)
2836 : value(refinement_case)
2837{
2838 // check that only those bits of
2839 // the given argument are set that
2840 // make sense for a given space
2841 // dimension
2842 Assert((refinement_case &
2844 refinement_case,
2845 ExcInvalidRefinementCase(refinement_case));
2846}
2847
2848
2849
2850template <int dim>
2851inline DEAL_II_HOST_DEVICE RefinementCase<dim>::operator std::uint8_t() const
2852{
2853 return value;
2854}
2855
2856
2857
2858template <int dim>
2861{
2862 return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2863}
2864
2865
2866
2867template <int dim>
2870{
2871 return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2872}
2873
2874
2875
2876template <int dim>
2879{
2880 return RefinementCase<dim>(static_cast<std::uint8_t>(
2882}
2883
2884
2885
2886template <int dim>
2887inline std::size_t
2889{
2890 return sizeof(RefinementCase<dim>);
2891}
2892
2893
2894
2895template <int dim>
2896template <class Archive>
2897inline void
2898RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2899{
2900 // serialization can't deal with bitfields, so copy from/to a full sized
2901 // std::uint8_t
2902 std::uint8_t uchar_value = value;
2903 ar &uchar_value;
2904 value = uchar_value;
2905}
2906
2907
2908
2909template <>
2910inline Point<1>
2911GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2912{
2913 AssertIndexRange(vertex, vertices_per_cell);
2914
2915 return Point<1>(static_cast<double>(vertex));
2916}
2917
2918
2919
2920template <>
2921inline Point<2>
2922GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2923{
2924 AssertIndexRange(vertex, vertices_per_cell);
2925
2926 return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2927}
2928
2929
2930
2931template <>
2932inline Point<3>
2933GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2934{
2935 AssertIndexRange(vertex, vertices_per_cell);
2936
2937 return {static_cast<double>(vertex % 2),
2938 static_cast<double>(vertex / 2 % 2),
2939 static_cast<double>(vertex / 4)};
2940}
2941
2942
2943
2944inline std::array<unsigned int, 0>
2946{
2947 return {{}};
2948}
2949
2950
2951
2952inline std::array<unsigned int, 1>
2954{
2955 return {{0}};
2956}
2957
2958
2959
2960template <int dim>
2963{
2965 0U, faces_per_cell);
2966}
2967
2968
2969
2970template <int dim>
2973{
2975 0U, vertices_per_cell);
2976}
2977
2978
2979
2980template <int dim>
2981inline Point<dim>
2982GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2983{
2985
2986 return {};
2987}
2988
2989
2990
2991template <>
2992inline unsigned int
2994{
2995 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2996
2997 return (p[0] <= 0.5 ? 0 : 1);
2998}
2999
3000
3001
3002template <>
3003inline unsigned int
3005{
3006 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
3007 Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
3008
3009 return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
3010}
3011
3012
3013
3014template <>
3015inline unsigned int
3017{
3018 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
3019 Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
3020 Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
3021
3022 return (p[0] <= 0.5 ?
3023 (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
3024 (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
3025}
3026
3027
3028template <int dim>
3029inline unsigned int
3031{
3033
3034 return 0;
3035}
3036
3037
3038
3039template <>
3040inline Point<1>
3042 const unsigned int child_index,
3043 const RefinementCase<1> refine_case)
3044
3045{
3046 AssertIndexRange(child_index, 2);
3048 (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3049
3050 return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
3051}
3052
3053
3054
3055template <>
3056inline Point<2>
3058 const unsigned int child_index,
3059 const RefinementCase<2> refine_case)
3060
3061{
3062 AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3063
3064 Point<2> point = p;
3065 switch (refine_case)
3066 {
3068 point[0] *= 2.0;
3069 if (child_index == 1)
3070 point[0] -= 1.0;
3071 break;
3073 point[1] *= 2.0;
3074 if (child_index == 1)
3075 point[1] -= 1.0;
3076 break;
3078 point *= 2.0;
3079 point -= unit_cell_vertex(child_index);
3080 break;
3081 default:
3083 }
3084
3085 return point;
3086}
3087
3088
3089
3090template <>
3091inline Point<3>
3093 const unsigned int child_index,
3094 const RefinementCase<3> refine_case)
3095
3096{
3097 AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3098
3099 Point<3> point = p;
3100 // there might be a cleverer way to do
3101 // this, but since this function is called
3102 // in very few places for initialization
3103 // purposes only, I don't care at the
3104 // moment
3105 switch (refine_case)
3106 {
3108 point[0] *= 2.0;
3109 if (child_index == 1)
3110 point[0] -= 1.0;
3111 break;
3113 point[1] *= 2.0;
3114 if (child_index == 1)
3115 point[1] -= 1.0;
3116 break;
3118 point[2] *= 2.0;
3119 if (child_index == 1)
3120 point[2] -= 1.0;
3121 break;
3123 point[0] *= 2.0;
3124 point[1] *= 2.0;
3125 if (child_index % 2 == 1)
3126 point[0] -= 1.0;
3127 if (child_index / 2 == 1)
3128 point[1] -= 1.0;
3129 break;
3131 // careful, this is slightly
3132 // different from xy and yz due to
3133 // different internal numbering of
3134 // children!
3135 point[0] *= 2.0;
3136 point[2] *= 2.0;
3137 if (child_index / 2 == 1)
3138 point[0] -= 1.0;
3139 if (child_index % 2 == 1)
3140 point[2] -= 1.0;
3141 break;
3143 point[1] *= 2.0;
3144 point[2] *= 2.0;
3145 if (child_index % 2 == 1)
3146 point[1] -= 1.0;
3147 if (child_index / 2 == 1)
3148 point[2] -= 1.0;
3149 break;
3151 point *= 2.0;
3152 point -= unit_cell_vertex(child_index);
3153 break;
3154 default:
3156 }
3157
3158 return point;
3159}
3160
3161
3162
3163template <int dim>
3164inline Point<dim>
3166 const Point<dim> & /*p*/,
3167 const unsigned int /*child_index*/,
3168 const RefinementCase<dim> /*refine_case*/)
3169
3170{
3172 return {};
3173}
3174
3175
3176
3177template <>
3178inline Point<1>
3180 const unsigned int child_index,
3181 const RefinementCase<1> refine_case)
3182
3183{
3184 AssertIndexRange(child_index, 2);
3186 (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3187
3188 return (p + unit_cell_vertex(child_index)) * 0.5;
3189}
3190
3191
3192
3193template <>
3194inline Point<3>
3196 const unsigned int child_index,
3197 const RefinementCase<3> refine_case)
3198
3199{
3200 AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3201
3202 Point<3> point = p;
3203 // there might be a cleverer way to do
3204 // this, but since this function is called
3205 // in very few places for initialization
3206 // purposes only, I don't care at the
3207 // moment
3208 switch (refine_case)
3209 {
3211 if (child_index == 1)
3212 point[0] += 1.0;
3213 point[0] *= 0.5;
3214 break;
3216 if (child_index == 1)
3217 point[1] += 1.0;
3218 point[1] *= 0.5;
3219 break;
3221 if (child_index == 1)
3222 point[2] += 1.0;
3223 point[2] *= 0.5;
3224 break;
3226 if (child_index % 2 == 1)
3227 point[0] += 1.0;
3228 if (child_index / 2 == 1)
3229 point[1] += 1.0;
3230 point[0] *= 0.5;
3231 point[1] *= 0.5;
3232 break;
3234 // careful, this is slightly
3235 // different from xy and yz due to
3236 // different internal numbering of
3237 // children!
3238 if (child_index / 2 == 1)
3239 point[0] += 1.0;
3240 if (child_index % 2 == 1)
3241 point[2] += 1.0;
3242 point[0] *= 0.5;
3243 point[2] *= 0.5;
3244 break;
3246 if (child_index % 2 == 1)
3247 point[1] += 1.0;
3248 if (child_index / 2 == 1)
3249 point[2] += 1.0;
3250 point[1] *= 0.5;
3251 point[2] *= 0.5;
3252 break;
3254 point += unit_cell_vertex(child_index);
3255 point *= 0.5;
3256 break;
3257 default:
3259 }
3260
3261 return point;
3262}
3263
3264
3265
3266template <>
3267inline Point<2>
3269 const unsigned int child_index,
3270 const RefinementCase<2> refine_case)
3271{
3272 AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3273
3274 Point<2> point = p;
3275 switch (refine_case)
3276 {
3278 if (child_index == 1)
3279 point[0] += 1.0;
3280 point[0] *= 0.5;
3281 break;
3283 if (child_index == 1)
3284 point[1] += 1.0;
3285 point[1] *= 0.5;
3286 break;
3288 point += unit_cell_vertex(child_index);
3289 point *= 0.5;
3290 break;
3291 default:
3293 }
3294
3295 return point;
3296}
3297
3298
3299
3300template <int dim>
3301inline Point<dim>
3303 const Point<dim> & /*p*/,
3304 const unsigned int /*child_index*/,
3305 const RefinementCase<dim> /*refine_case*/)
3306{
3308 return {};
3309}
3310
3311
3312
3313template <int dim>
3314inline bool
3316{
3318 return false;
3319}
3320
3321template <>
3322inline bool
3324{
3325 return (p[0] >= 0.) && (p[0] <= 1.);
3326}
3327
3328
3329
3330template <>
3331inline bool
3333{
3334 return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
3335}
3336
3337
3338
3339template <>
3340inline bool
3342{
3343 return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
3344 (p[2] >= 0.) && (p[2] <= 1.);
3345}
3346
3347
3348
3349template <int dim>
3350inline bool
3352{
3354 return false;
3355}
3356
3357template <>
3358inline bool
3359GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
3360{
3361 return (p[0] >= -eps) && (p[0] <= 1. + eps);
3362}
3363
3364
3365
3366template <>
3367inline bool
3368GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
3369{
3370 const double l = -eps, u = 1 + eps;
3371 return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
3372}
3373
3374
3375
3376template <>
3377inline bool
3378GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
3379{
3380 const double l = -eps, u = 1.0 + eps;
3381 return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
3382 (p[2] >= l) && (p[2] <= u);
3383}
3384
3385
3386
3387template <>
3388inline unsigned int
3389GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
3390 const unsigned int vertex)
3391{
3392 (void)line;
3393 AssertIndexRange(line, lines_per_cell);
3394 AssertIndexRange(vertex, 2);
3395
3396 return vertex;
3397}
3398
3399
3400template <>
3401inline unsigned int
3402GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
3403 const unsigned int vertex)
3404{
3405 constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
3406 return cell_vertices[line][vertex];
3407}
3408
3409
3410
3411template <>
3412inline unsigned int
3413GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
3414 const unsigned int vertex)
3415{
3416 AssertIndexRange(line, lines_per_cell);
3417 AssertIndexRange(vertex, 2);
3418
3419 constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3420 {1, 3},
3421 {0, 1},
3422 {2, 3},
3423 {4, 6}, // top face
3424 {5, 7},
3425 {4, 5},
3426 {6, 7},
3427 {0,
3428 4}, // connects of bottom
3429 {1, 5}, // top face
3430 {2, 6},
3431 {3, 7}};
3432 return vertices[line][vertex];
3433}
3434
3435
3436template <>
3437inline unsigned int
3438GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3439{
3442}
3443
3444template <>
3445inline unsigned int
3446GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3447 const bool face_orientation,
3448 const bool face_flip,
3449 const bool face_rotation)
3450{
3452
3453 // set up a table to make sure that
3454 // we handle non-standard faces correctly
3455 //
3456 // so set up a table that for each vertex (of
3457 // a quad in standard position) describes
3458 // which vertex to take
3459 //
3460 // first index: four vertices 0...3
3461 //
3462 // second index: face_orientation; 0:
3463 // opposite normal, 1: standard
3464 //
3465 // third index: face_flip; 0: standard, 1:
3466 // face rotated by 180 degrees
3467 //
3468 // forth index: face_rotation: 0: standard,
3469 // 1: face rotated by 90 degrees
3470
3471 constexpr unsigned int vertex_translation[4][2][2][2] = {
3472 {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3473 // face_rotation=false and true
3474 {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3475 // face_rotation=false and true
3476 {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3477 // face_rotation=false and true
3478 {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3479 // face_rotation=false and true
3480
3481 {{{2, 3}, // vertex 1 ...
3482 {1, 0}},
3483 {{1, 0}, {2, 3}}},
3484
3485 {{{1, 0}, // vertex 2 ...
3486 {2, 3}},
3487 {{2, 3}, {1, 0}}},
3488
3489 {{{3, 1}, // vertex 3 ...
3490 {0, 2}},
3491 {{3, 1}, {0, 2}}}};
3492
3493 return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3494}
3495
3496
3497
3498template <int dim>
3499inline unsigned int
3500GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3501 const bool,
3502 const bool,
3503 const bool)
3504{
3505 Assert(dim > 1, ExcImpossibleInDim(dim));
3507 return vertex;
3508}
3509
3510template <int dim>
3511inline unsigned int
3513{
3514 constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3515 0, 2, 2, 4, 2, 4, 4, 8};
3516
3517 return n_children[ref_case];
3518}
3519
3520
3521
3522template <int dim>
3523inline unsigned int
3525{
3527 return 0;
3528}
3529
3530template <>
3531inline unsigned int
3533{
3534 Assert(false, ExcImpossibleInDim(1));
3535 return 0;
3536}
3537
3538template <>
3539inline unsigned int
3541{
3542 return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3543}
3544
3545
3546
3547template <>
3548inline unsigned int
3550{
3551 const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3552 0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3553 return nsubs[subface_case];
3554}
3555
3556
3557
3558template <int dim>
3559inline double
3561 const unsigned int)
3562{
3564 return 0.;
3565}
3566
3567template <>
3568inline double
3570 const unsigned int)
3571{
3572 return 1;
3573}
3574
3575
3576template <>
3577inline double
3579 const unsigned int)
3580{
3581 double ratio = 1;
3582 switch (subface_case)
3583 {
3585 // Here, an
3586 // Assert(false,ExcInternalError())
3587 // would be the right
3588 // choice, but
3589 // unfortunately the
3590 // current function is
3591 // also called for faces
3592 // without children (see
3593 // tests/fe/mapping.cc).
3594 // Assert(false, ExcMessage("Face has no subfaces."));
3595 // Furthermore, assign
3596 // following value as
3597 // otherwise the
3598 // bits/volume_x tests
3599 // break
3601 break;
3603 ratio = 0.5;
3604 break;
3605 default:
3606 // there should be no
3607 // cases left
3609 break;
3610 }
3611
3612 return ratio;
3613}
3614
3615
3616template <>
3617inline double
3619 const unsigned int subface_no)
3620{
3621 double ratio = 1;
3622 switch (subface_case)
3623 {
3625 // Here, an
3626 // Assert(false,ExcInternalError())
3627 // would be the right
3628 // choice, but
3629 // unfortunately the
3630 // current function is
3631 // also called for faces
3632 // without children (see
3633 // tests/bits/mesh_3d_16.cc). Add
3634 // following switch to
3635 // avoid diffs in
3636 // tests/bits/mesh_3d_16
3638 break;
3641 ratio = 0.5;
3642 break;
3646 ratio = 0.25;
3647 break;
3650 if (subface_no < 2)
3651 ratio = 0.25;
3652 else
3653 ratio = 0.5;
3654 break;
3657 if (subface_no == 0)
3658 ratio = 0.5;
3659 else
3660 ratio = 0.25;
3661 break;
3662 default:
3663 // there should be no
3664 // cases left
3666 break;
3667 }
3668
3669 return ratio;
3670}
3671
3672
3673
3674template <int dim>
3676 const RefinementCase<dim> &,
3677 const unsigned int,
3678 const bool,
3679 const bool,
3680 const bool)
3681{
3683 return RefinementCase<dim - 1>::no_refinement;
3684}
3685
3686template <>
3688 const RefinementCase<1> &,
3689 const unsigned int,
3690 const bool,
3691 const bool,
3692 const bool)
3693{
3694 Assert(false, ExcImpossibleInDim(1));
3695
3697}
3698
3699
3700template <>
3701inline RefinementCase<1>
3703 const RefinementCase<2> &cell_refinement_case,
3704 const unsigned int face_no,
3705 const bool,
3706 const bool,
3707 const bool)
3708{
3709 const unsigned int dim = 2;
3710 AssertIndexRange(cell_refinement_case,
3713
3714 // simple special case
3715 if (cell_refinement_case == RefinementCase<dim>::cut_xy)
3717
3718 const RefinementCase<dim - 1>
3721 {RefinementCase<dim - 1>::no_refinement, // no_refinement
3722 RefinementCase<dim - 1>::no_refinement},
3723
3724 {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3725
3726 {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3727
3728 {RefinementCase<dim - 1>::cut_x, // cut_xy
3729 RefinementCase<dim - 1>::cut_x}};
3730
3731 return ref_cases[cell_refinement_case][face_no / 2];
3732}
3733
3734
3735template <>
3736inline RefinementCase<2>
3738 const RefinementCase<3> &cell_refinement_case,
3739 const unsigned int face_no,
3740 const bool face_orientation,
3741 const bool /*face_flip*/,
3742 const bool face_rotation)
3743{
3744 const unsigned int dim = 3;
3745 AssertIndexRange(cell_refinement_case,
3748
3749 // simple special case
3750 if (cell_refinement_case == RefinementCase<dim>::cut_xyz)
3751 return RefinementCase<dim - 1>::cut_xy;
3752
3753 const RefinementCase<dim - 1>
3756 {RefinementCase<dim - 1>::no_refinement, // no_refinement
3757 RefinementCase<dim - 1>::no_refinement,
3758 RefinementCase<dim - 1>::no_refinement},
3759
3760 {RefinementCase<dim - 1>::no_refinement, // cut_x
3761 RefinementCase<dim - 1>::cut_y,
3762 RefinementCase<dim - 1>::cut_x},
3763
3764 {RefinementCase<dim - 1>::cut_x, // cut_y
3765 RefinementCase<dim - 1>::no_refinement,
3766 RefinementCase<dim - 1>::cut_y},
3767
3768 {RefinementCase<dim - 1>::cut_x, // cut_xy
3769 RefinementCase<dim - 1>::cut_y,
3770 RefinementCase<dim - 1>::cut_xy},
3771
3772 {RefinementCase<dim - 1>::cut_y, // cut_z
3773 RefinementCase<dim - 1>::cut_x,
3774 RefinementCase<dim - 1>::no_refinement},
3775
3776 {RefinementCase<dim - 1>::cut_y, // cut_xz
3777 RefinementCase<dim - 1>::cut_xy,
3778 RefinementCase<dim - 1>::cut_x},
3779
3780 {RefinementCase<dim - 1>::cut_xy, // cut_yz
3781 RefinementCase<dim - 1>::cut_x,
3782 RefinementCase<dim - 1>::cut_y},
3783
3784 {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3785 RefinementCase<dim - 1>::cut_xy,
3786 RefinementCase<dim - 1>::cut_xy},
3787 };
3788
3789 const RefinementCase<dim - 1> ref_case =
3790 ref_cases[cell_refinement_case][face_no / 2];
3791
3792 const RefinementCase<dim - 1> flip[4] = {
3793 RefinementCase<dim - 1>::no_refinement,
3794 RefinementCase<dim - 1>::cut_y,
3795 RefinementCase<dim - 1>::cut_x,
3796 RefinementCase<dim - 1>::cut_xy};
3797
3798 // correct the ref_case for face_orientation
3799 // and face_rotation. for face_orientation,
3800 // 'true' is the default value whereas for
3801 // face_rotation, 'false' is standard. If
3802 // <tt>face_rotation==face_orientation</tt>,
3803 // then one of them is non-standard and we
3804 // have to swap cut_x and cut_y, otherwise no
3805 // change is necessary. face_flip has no
3806 // influence. however, in order to keep the
3807 // interface consistent with other functions,
3808 // we still include it as an argument to this
3809 // function
3810 return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3811}
3812
3813
3814
3815template <int dim>
3816inline RefinementCase<1>
3818 const unsigned int)
3819{
3822}
3823
3824template <>
3825inline RefinementCase<1>
3827 const RefinementCase<1> &cell_refinement_case,
3828 const unsigned int line_no)
3829{
3830 (void)line_no;
3831 const unsigned int dim = 1;
3832 (void)dim;
3833 AssertIndexRange(cell_refinement_case,
3836
3837 return cell_refinement_case;
3838}
3839
3840
3841template <>
3842inline RefinementCase<1>
3844 const RefinementCase<2> &cell_refinement_case,
3845 const unsigned int line_no)
3846{
3847 // Assertions are in face_refinement_case()
3848 return face_refinement_case(cell_refinement_case, line_no);
3849}
3850
3851
3852template <>
3853inline RefinementCase<1>
3855 const RefinementCase<3> &cell_refinement_case,
3856 const unsigned int line_no)
3857{
3858 const unsigned int dim = 3;
3859 AssertIndexRange(cell_refinement_case,
3862
3863 // simple special case
3864 if (cell_refinement_case == RefinementCase<dim>::cut_xyz)
3866
3867 // array indicating, which simple refine
3868 // case cuts a line in direction x, y or
3869 // z. For example, cut_y and everything
3870 // containing cut_y (cut_xy, cut_yz,
3871 // cut_xyz) cuts lines, which are in y
3872 // direction.
3873 const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3876
3877 // order the direction of lines
3878 // 0->x, 1->y, 2->z
3879 const unsigned int direction[lines_per_cell] = {
3880 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3881
3882 return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3884 RefinementCase<1>::no_refinement);
3885}
3886
3887
3888
3889template <int dim>
3893 const unsigned int,
3894 const bool,
3895 const bool,
3896 const bool)
3897{
3899
3901}
3902
3903template <>
3904inline RefinementCase<1>
3906 const RefinementCase<0> &,
3907 const unsigned int,
3908 const bool,
3909 const bool,
3910 const bool)
3911{
3912 const unsigned int dim = 1;
3913 Assert(false, ExcImpossibleInDim(dim));
3914
3916}
3917
3918
3919template <>
3920inline RefinementCase<2>
3922 const RefinementCase<1> &face_refinement_case,
3923 const unsigned int face_no,
3924 const bool,
3925 const bool,
3926 const bool)
3927{
3928 const unsigned int dim = 2;
3929 AssertIndexRange(face_refinement_case,
3932
3933 if (face_refinement_case == RefinementCase<dim>::cut_x)
3934 return (face_no / 2) != 0u ? RefinementCase<dim>::cut_x :
3936 else
3938}
3939
3940
3941template <>
3942inline RefinementCase<3>
3944 const RefinementCase<2> &face_refinement_case,
3945 const unsigned int face_no,
3946 const bool face_orientation,
3947 const bool /*face_flip*/,
3948 const bool face_rotation)
3949{
3950 const unsigned int dim = 3;
3951 AssertIndexRange(face_refinement_case,
3954
3959
3960 // correct the face_refinement_case for
3961 // face_orientation and face_rotation. for
3962 // face_orientation, 'true' is the default
3963 // value whereas for face_rotation, 'false'
3964 // is standard. If
3965 // <tt>face_rotation==face_orientation</tt>,
3966 // then one of them is non-standard and we
3967 // have to swap cut_x and cut_y, otherwise no
3968 // change is necessary. face_flip has no
3969 // influence. however, in order to keep the
3970 // interface consistent with other functions,
3971 // we still include it as an argument to this
3972 // function
3973 const RefinementCase<dim - 1> std_face_ref =
3974 (face_orientation == face_rotation) ? flip[face_refinement_case] :
3975 face_refinement_case;
3976
3977 const RefinementCase<dim> face_to_cell[3][4] = {
3978 {RefinementCase<dim>::no_refinement, // faces 0 and 1
3979 RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3982
3983 {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3984 // "exchanged on faces 2 and 3")
3988
3989 {RefinementCase<dim>::no_refinement, // faces 4 and 5
3993
3994 return face_to_cell[face_no / 2][std_face_ref];
3995}
3996
3997
3998
3999template <int dim>
4002 const unsigned int)
4003{
4005
4007}
4008
4009template <>
4010inline RefinementCase<1>
4012 const unsigned int line_no)
4013{
4014 (void)line_no;
4015 AssertIndexRange(line_no, 1);
4016
4018}
4019
4020
4021template <>
4022inline RefinementCase<2>
4024 const unsigned int line_no)
4025{
4026 const unsigned int dim = 2;
4027 (void)dim;
4029
4030 return (line_no / 2) != 0u ? RefinementCase<2>::cut_x :
4032}
4033
4034
4035template <>
4036inline RefinementCase<3>
4038 const unsigned int line_no)
4039{
4040 const unsigned int dim = 3;
4042
4043 const RefinementCase<dim> ref_cases[6] = {
4044 RefinementCase<dim>::cut_y, // lines 0 and 1
4045 RefinementCase<dim>::cut_x, // lines 2 and 3
4046 RefinementCase<dim>::cut_y, // lines 4 and 5
4047 RefinementCase<dim>::cut_x, // lines 6 and 7
4048 RefinementCase<dim>::cut_z, // lines 8 and 9
4049 RefinementCase<dim>::cut_z}; // lines 10 and 11
4050
4051 return ref_cases[line_no / 2];
4052}
4053
4054
4055
4056template <>
4057inline unsigned int
4058GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
4059 const bool face_orientation,
4060 const bool face_flip,
4061 const bool face_rotation)
4062{
4064
4065 // set up a table to make sure that
4066 // we handle non-standard faces correctly
4067 //
4068 // so set up a table that for each vertex (of
4069 // a quad in standard position) describes
4070 // which vertex to take
4071 //
4072 // first index: four vertices 0...3
4073 //
4074 // second index: face_orientation; 0:
4075 // opposite normal, 1: standard
4076 //
4077 // third index: face_flip; 0: standard, 1:
4078 // face rotated by 180 degrees
4079 //
4080 // forth index: face_rotation: 0: standard,
4081 // 1: face rotated by 90 degrees
4082
4083 const unsigned int vertex_translation[4][2][2][2] = {
4084 {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
4085 // face_rotation=false and true
4086 {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
4087 // face_rotation=false and true
4088 {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
4089 // face_rotation=false and true
4090 {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
4091 // face_rotation=false and true
4092
4093 {{{2, 3}, // vertex 1 ...
4094 {1, 0}},
4095 {{1, 3}, {2, 0}}},
4096
4097 {{{1, 0}, // vertex 2 ...
4098 {2, 3}},
4099 {{2, 0}, {1, 3}}},
4100
4101 {{{3, 1}, // vertex 3 ...
4102 {0, 2}},
4103 {{3, 2}, {0, 1}}}};
4104
4105 return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
4106}
4107
4108
4109
4110template <int dim>
4111inline unsigned int
4112GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
4113 const bool,
4114 const bool,
4115 const bool)
4116{
4117 Assert(dim > 1, ExcImpossibleInDim(dim));
4119 return vertex;
4120}
4121
4122
4123
4124template <>
4125inline unsigned int
4126GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
4127 const bool face_orientation,
4128 const bool face_flip,
4129 const bool face_rotation)
4130{
4132
4133
4134 // make sure we handle
4135 // non-standard faces correctly
4136 //
4137 // so set up a table that for each line (of a
4138 // quad) describes which line to take
4139 //
4140 // first index: four lines 0...3
4141 //
4142 // second index: face_orientation; 0:
4143 // opposite normal, 1: standard
4144 //
4145 // third index: face_flip; 0: standard, 1:
4146 // face rotated by 180 degrees
4147 //
4148 // forth index: face_rotation: 0: standard,
4149 // 1: face rotated by 90 degrees
4150
4151 const unsigned int line_translation[4][2][2][2] = {
4152 {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4153 // face_rotation=false and true
4154 {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4155 // face_rotation=false and true
4156 {{0, 3}, // line 0, face_orientation=true, face_flip=false,
4157 // face_rotation=false and true
4158 {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
4159 // face_rotation=false and true
4160
4161 {{{3, 1}, // line 1 ...
4162 {2, 0}},
4163 {{1, 2}, {0, 3}}},
4164
4165 {{{0, 3}, // line 2 ...
4166 {1, 2}},
4167 {{2, 0}, {3, 1}}},
4168
4169 {{{1, 2}, // line 3 ...
4170 {0, 3}},
4171 {{3, 1}, {2, 0}}}};
4172
4173 return line_translation[line][face_orientation][face_flip][face_rotation];
4174}
4175
4176
4177
4178template <int dim>
4179inline unsigned int
4180GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
4181 const bool,
4182 const bool,
4183 const bool)
4184{
4186 return line;
4187}
4188
4189
4190
4191template <>
4192inline unsigned int
4193GeometryInfo<2>::standard_to_real_line_vertex(const unsigned int vertex,
4194 const bool line_orientation)
4195{
4196 return line_orientation ? vertex : (1 - vertex);
4197}
4198
4199
4200
4201template <int dim>
4202inline unsigned int
4203GeometryInfo<dim>::standard_to_real_line_vertex(const unsigned int vertex,
4204 const bool)
4205{
4207 return vertex;
4208}
4209
4210
4211
4212template <>
4213inline std::array<unsigned int, 2>
4215 const unsigned int vertex)
4216{
4217 return {{vertex % 2, vertex / 2}};
4218}
4219
4220
4221
4222template <int dim>
4223inline std::array<unsigned int, 2>
4225 const unsigned int vertex)
4226{
4228 (void)vertex;
4229 return {{0, 0}};
4230}
4231
4232
4233
4234template <>
4235inline std::array<unsigned int, 2>
4237{
4238 // set up a table that for each
4239 // line describes a) from which
4240 // quad to take it, b) which line
4241 // therein it is if the face is
4242 // oriented correctly
4243 static const unsigned int lookup_table[GeometryInfo<3>::lines_per_cell][2] = {
4244 {4, 0}, // take first four lines from bottom face
4245 {4, 1},
4246 {4, 2},
4247 {4, 3},
4248
4249 {5, 0}, // second four lines from top face
4250 {5, 1},
4251 {5, 2},
4252 {5, 3},
4253
4254 {0, 0}, // the rest randomly
4255 {1, 0},
4256 {0, 1},
4257 {1, 1}};
4258
4259 return {{lookup_table[i][0], lookup_table[i][1]}};
4260}
4261
4262
4263
4264template <int dim>
4265inline std::array<unsigned int, 2>
4267{
4269 (void)line;
4270 return {{0, 0}};
4271}
4272
4273
4274
4275template <>
4276inline std::array<unsigned int, 2>
4278 const unsigned int vertex)
4279{
4280 // get the corner indices by asking either the bottom or the top face for its
4281 // vertices. handle non-standard faces by calling the vertex reordering
4282 // function from GeometryInfo
4283
4284 // bottom face (4) for first four vertices, top face (5) for the rest
4285 return {{4 + vertex / 4, vertex % 4}};
4286}
4287
4288
4289
4290template <int dim>
4291inline std::array<unsigned int, 2>
4293 const unsigned int vertex)
4294{
4296 (void)vertex;
4297 return {{0, 0}};
4298}
4299
4300
4301
4302template <>
4303inline unsigned int
4304GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
4305 const bool face_orientation,
4306 const bool face_flip,
4307 const bool face_rotation)
4308{
4310
4311
4312 // make sure we handle
4313 // non-standard faces correctly
4314 //
4315 // so set up a table that for each line (of a
4316 // quad) describes which line to take
4317 //
4318 // first index: four lines 0...3
4319 //
4320 // second index: face_orientation; 0:
4321 // opposite normal, 1: standard
4322 //
4323 // third index: face_flip; 0: standard, 1:
4324 // face rotated by 180 degrees
4325 //
4326 // forth index: face_rotation: 0: standard,
4327 // 1: face rotated by 90 degrees
4328
4329 const unsigned int line_translation[4][2][2][2] = {
4330 {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4331 // face_rotation=false and true
4332 {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4333 // face_rotation=false and true
4334 {{0, 2}, // line 0, face_orientation=true, face_flip=false,
4335 // face_rotation=false and true
4336 {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
4337 // face_rotation=false and true
4338
4339 {{{3, 1}, // line 1 ...
4340 {2, 0}},
4341 {{1, 3}, {0, 2}}},
4342
4343 {{{0, 3}, // line 2 ...
4344 {1, 2}},
4345 {{2, 1}, {3, 0}}},
4346
4347 {{{1, 2}, // line 3 ...
4348 {0, 3}},
4349 {{3, 0}, {2, 1}}}};
4350
4351 return line_translation[line][face_orientation][face_flip][face_rotation];
4352}
4353
4354
4355
4356template <int dim>
4357inline unsigned int
4358GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
4359 const bool,
4360 const bool,
4361 const bool)
4362{
4364 return line;
4365}
4366
4367
4368
4369template <>
4370inline unsigned int
4372 const unsigned int face,
4373 const unsigned int subface,
4374 const bool,
4375 const bool,
4376 const bool,
4377 const RefinementCase<0> &)
4378{
4379 (void)subface;
4380 AssertIndexRange(face, faces_per_cell);
4381 AssertIndexRange(subface, max_children_per_face);
4382
4383 return face;
4384}
4385
4386
4387
4388template <>
4389inline unsigned int
4391 const unsigned int face,
4392 const unsigned int subface,
4393 const bool face_orientation,
4394 const bool /*face_flip*/,
4395 const bool /*face_rotation*/,
4396 const RefinementCase<1> &)
4397{
4398 AssertIndexRange(face, faces_per_cell);
4399 AssertIndexRange(subface, max_children_per_face);
4400
4401 // always return the child adjacent to the specified
4402 // subface. if the face of a cell is not refined, don't
4403 // throw an assertion but deliver the child adjacent to
4404 // the face nevertheless, i.e. deliver the child of
4405 // this cell adjacent to the subface of a possibly
4406 // refined neighbor. this simplifies setting neighbor
4407 // information in execute_refinement.
4408 constexpr unsigned int
4409 subcells[/* possible face orientation */ 2]
4410 [/* number of different ways to refine a cell */ 4]
4411 [/* faces_per_cell */ 4][/* max_children_per_face */ 2] = {
4412 {
4413 // Flipped orientation (face_orientation = false)
4414 {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
4415 {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
4416 {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
4417 },
4418 {
4419 // Normal orientation (face_orientation = true)
4420 {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
4421 {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
4422 {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
4423 }};
4424
4425 return subcells[face_orientation][ref_case - 1][face][subface];
4426}
4427
4428
4429
4430template <>
4431inline unsigned int
4433 const unsigned int face,
4434 const unsigned int subface,
4435 const bool face_orientation,
4436 const bool face_flip,
4437 const bool face_rotation,
4438 const RefinementCase<2> &face_ref_case)
4439{
4440 const unsigned int dim = 3;
4441
4443 ExcMessage("Cell has no children."));
4444 AssertIndexRange(face, faces_per_cell);
4445 if (!(subface == 0 &&
4447 {
4448 AssertIndexRange(subface,
4449 GeometryInfo<dim - 1>::n_children(face_ref_case));
4450 }
4451
4452 // invalid number used for invalid cases,
4453 // e.g. when the children are more refined at
4454 // a given face than the face itself
4455 const unsigned int e = numbers::invalid_unsigned_int;
4456
4457 // the whole process of finding a child cell
4458 // at a given subface considering the
4459 // possibly anisotropic refinement cases of
4460 // the cell and the face as well as
4461 // orientation, flip and rotation of the face
4462 // is quite complicated. thus, we break it
4463 // down into several steps.
4464
4465 // first step: convert the given face refine
4466 // case to a face refine case concerning the
4467 // face in standard orientation (, flip and
4468 // rotation). This only affects cut_x and
4469 // cut_y
4470 const RefinementCase<dim - 1> flip[4] = {
4471 RefinementCase<dim - 1>::no_refinement,
4472 RefinementCase<dim - 1>::cut_y,
4473 RefinementCase<dim - 1>::cut_x,
4474 RefinementCase<dim - 1>::cut_xy};
4475 // for face_orientation, 'true' is the
4476 // default value whereas for face_rotation,
4477 // 'false' is standard. If
4478 // <tt>face_rotation==face_orientation</tt>,
4479 // then one of them is non-standard and we
4480 // have to swap cut_x and cut_y, otherwise no
4481 // change is necessary.
4482 const RefinementCase<dim - 1> std_face_ref =
4483 (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4484
4485 // second step: convert the given subface
4486 // index to the one for a standard face
4487 // respecting face_orientation, face_flip and
4488 // face_rotation
4489
4490 // first index: face_ref_case
4491 // second index: face_orientation
4492 // third index: face_flip
4493 // forth index: face_rotation
4494 // fifth index: subface index
4495 const unsigned int subface_exchange[4][2][2][2][4] = {
4496 // no_refinement (subface 0 stays 0,
4497 // all others are invalid)
4498 {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4499 {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4500 // cut_x (here, if the face is only
4501 // rotated OR only falsely oriented,
4502 // then subface 0 of the non-standard
4503 // face does NOT correspond to one of
4504 // the subfaces of a standard
4505 // face. Thus we indicate the subface
4506 // which is located at the lower left
4507 // corner (the origin of the face's
4508 // local coordinate system) with
4509 // '0'. The rest of this issue is
4510 // taken care of using the above
4511 // conversion to a 'standard face
4512 // refine case')
4513 {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4514 {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4515 // cut_y (the same applies as for
4516 // cut_x)
4517 {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4518 {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4519 // cut_xyz: this information is
4520 // identical to the information
4521 // returned by
4522 // GeometryInfo<3>::real_to_standard_face_vertex()
4523 {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4524 // face_rotation=false, subfaces 0,1,2,3
4525 {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4526 // face_rotation=true, subfaces 0,1,2,3
4527 {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4528 // face_rotation=false, subfaces 0,1,2,3
4529 {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4530 // face_rotation=true, subfaces 0,1,2,3
4531 {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4532 // face_rotation=false, subfaces 0,1,2,3
4533 {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4534 // face_rotation=true, subfaces 0,1,2,3
4535 {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4536 // face_rotation=false, subfaces 0,1,2,3
4537 {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4538 // face_rotation=true, subfaces 0,1,2,3
4539
4540 const unsigned int std_subface =
4541 subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4542 [subface];
4543 Assert(std_subface != e, ExcInternalError());
4544
4545 // third step: these are the children, which
4546 // can be found at the given subfaces of an
4547 // isotropically refined (standard) face
4548 //
4549 // first index: (refinement_case-1)
4550 // second index: face_index
4551 // third index: subface_index (isotropic refinement)
4552 const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4553 [max_children_per_face] = {
4554 // cut_x
4555 {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4556 {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4557 {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4558 {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4559 {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4560 {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4561 // cut_y
4562 {{0, 1, 0, 1},
4563 {0, 1, 0, 1},
4564 {0, 0, 0, 0},
4565 {1, 1, 1, 1},
4566 {0, 0, 1, 1},
4567 {0, 0, 1, 1}},
4568 // cut_xy
4569 {{0, 2, 0, 2},
4570 {1, 3, 1, 3},
4571 {0, 0, 1, 1},
4572 {2, 2, 3, 3},
4573 {0, 1, 2, 3},
4574 {0, 1, 2, 3}},
4575 // cut_z
4576 {{0, 0, 1, 1},
4577 {0, 0, 1, 1},
4578 {0, 1, 0, 1},
4579 {0, 1, 0, 1},
4580 {0, 0, 0, 0},
4581 {1, 1, 1, 1}},
4582 // cut_xz
4583 {{0, 0, 1, 1},
4584 {2, 2, 3, 3},
4585 {0, 1, 2, 3},
4586 {0, 1, 2, 3},
4587 {0, 2, 0, 2},
4588 {1, 3, 1, 3}},
4589 // cut_yz
4590 {{0, 1, 2, 3},
4591 {0, 1, 2, 3},
4592 {0, 2, 0, 2},
4593 {1, 3, 1, 3},
4594 {0, 0, 1, 1},
4595 {2, 2, 3, 3}},
4596 // cut_xyz
4597 {{0, 2, 4, 6},
4598 {1, 3, 5, 7},
4599 {0, 4, 1, 5},
4600 {2, 6, 3, 7},
4601 {0, 1, 2, 3},
4602 {4, 5, 6, 7}}};
4603
4604 // forth step: check, whether the given face
4605 // refine case is valid for the given cell
4606 // refine case. this is the case, if the
4607 // given face refine case is at least as
4608 // refined as the face is for the given cell
4609 // refine case
4610
4611 // note, that we are considering standard
4612 // face refinement cases here and thus must
4613 // not pass the given orientation, flip and
4614 // rotation flags
4615 if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4616 face_refinement_case(ref_case, face))
4617 {
4618 // all is fine. for anisotropic face
4619 // refine cases, select one of the
4620 // isotropic subfaces which neighbors the
4621 // same child
4622
4623 // first index: (standard) face refine case
4624 // second index: subface index
4625 const unsigned int equivalent_iso_subface[4][4] = {
4626 {0, e, e, e}, // no_refinement
4627 {0, 3, e, e}, // cut_x
4628 {0, 3, e, e}, // cut_y
4629 {0, 1, 2, 3}}; // cut_xy
4630
4631 const unsigned int equ_std_subface =
4632 equivalent_iso_subface[std_face_ref][std_subface];
4633 Assert(equ_std_subface != e, ExcInternalError());
4634
4635 return iso_children[ref_case - 1][face][equ_std_subface];
4636 }
4637 else
4638 {
4639 // the face_ref_case was too coarse,
4640 // throw an error
4641 Assert(false,
4642 ExcMessage("The face RefineCase is too coarse "
4643 "for the given cell RefineCase."));
4644 }
4645 // we only get here in case of an error
4646 return e;
4647}
4648
4649
4650
4651template <>
4652inline unsigned int
4654 const unsigned int,
4655 const unsigned int,
4656 const bool,
4657 const bool,
4658 const bool,
4659 const RefinementCase<3> &)
4660{
4663}
4664
4665
4666
4667template <>
4668inline unsigned int
4669GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4670 const unsigned int line,
4671 const bool,
4672 const bool,
4673 const bool)
4674{
4675 (void)line;
4676 AssertIndexRange(face, faces_per_cell);
4677 AssertIndexRange(line, lines_per_face);
4678
4679 // The face is a line itself.
4680 return face;
4681}
4682
4683
4684
4685template <>
4686inline unsigned int
4687GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4688 const unsigned int line,
4689 const bool face_orientation,
4690 const bool face_flip,
4691 const bool face_rotation)
4692{
4693 AssertIndexRange(face, faces_per_cell);
4694 AssertIndexRange(line, lines_per_face);
4695
4696 const unsigned lines[faces_per_cell][lines_per_face] = {
4697 {8, 10, 0, 4}, // left face
4698 {9, 11, 1, 5}, // right face
4699 {2, 6, 8, 9}, // front face
4700 {3, 7, 10, 11}, // back face
4701 {0, 1, 2, 3}, // bottom face
4702 {4, 5, 6, 7}}; // top face
4703 return lines[face][real_to_standard_face_line(
4704 line, face_orientation, face_flip, face_rotation)];
4705}
4706
4707
4708
4709inline unsigned int
4710GeometryInfo<0>::face_to_cell_lines(const unsigned int,
4711 const unsigned int,
4712 const bool,
4713 const bool,
4714 const bool)
4715{
4718}
4719
4720
4721
4722template <int dim>
4723inline unsigned int
4724GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4725 const unsigned int,
4726 const bool,
4727 const bool,
4728 const bool)
4729{
4732}
4733
4734
4735
4736template <int dim>
4737inline unsigned int
4738GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4739 const unsigned int vertex,
4740 const bool face_orientation,
4741 const bool face_flip,
4742 const bool face_rotation)
4743{
4744 return child_cell_on_face(RefinementCase<dim>::isotropic_refinement,
4745 face,
4746 vertex,
4747 face_orientation,
4748 face_flip,
4749 face_rotation);
4750}
4751
4752
4753
4754inline unsigned int
4756 const unsigned int,
4757 const bool,
4758 const bool,
4759 const bool)
4760{
4763}
4764
4765
4766
4767template <int dim>
4768template <typename Number>
4769inline Point<dim, Number>
4771{
4773 for (unsigned int i = 0; i < dim; ++i)
4774 p[i] = std::min(std::max(q[i], Number(0.)), Number(1.));
4775
4776 return p;
4777}
4778
4779
4780
4781template <int dim>
4782inline double
4784{
4785 double result = 0.0;
4786
4787 for (unsigned int i = 0; i < dim; ++i)
4788 {
4789 result = std::max(result, -p[i]);
4790 result = std::max(result, p[i] - 1.);
4791 }
4792
4793 return result;
4794}
4795
4796
4797
4798template <int dim>
4799inline double
4801 const unsigned int i)
4802{
4804
4805 switch (dim)
4806 {
4807 case 1:
4808 {
4809 const double x = xi[0];
4810 switch (i)
4811 {
4812 case 0:
4813 return 1 - x;
4814 case 1:
4815 return x;
4816 }
4817 break;
4818 }
4819
4820 case 2:
4821 {
4822 const double x = xi[0];
4823 const double y = xi[1];
4824 switch (i)
4825 {
4826 case 0:
4827 return (1 - x) * (1 - y);
4828 case 1:
4829 return x * (1 - y);
4830 case 2:
4831 return (1 - x) * y;
4832 case 3:
4833 return x * y;
4834 }
4835 break;
4836 }
4837
4838 case 3:
4839 {
4840 const double x = xi[0];
4841 const double y = xi[1];
4842 const double z = xi[2];
4843 switch (i)
4844 {
4845 case 0:
4846 return (1 - x) * (1 - y) * (1 - z);
4847 case 1:
4848 return x * (1 - y) * (1 - z);
4849 case 2:
4850 return (1 - x) * y * (1 - z);
4851 case 3:
4852 return x * y * (1 - z);
4853 case 4:
4854 return (1 - x) * (1 - y) * z;
4855 case 5:
4856 return x * (1 - y) * z;
4857 case 6:
4858 return (1 - x) * y * z;
4859 case 7:
4860 return x * y * z;
4861 }
4862 break;
4863 }
4864
4865 default:
4867 }
4868 return -1e9;
4869}
4870
4871
4872
4873template <>
4875 const Point<1> &,
4876 const unsigned int i)
4877{
4879
4880 switch (i)
4881 {
4882 case 0:
4883 return Point<1>(-1.);
4884 case 1:
4885 return Point<1>(1.);
4886 }
4887
4888 return Point<1>(-1e9);
4889}
4890
4891
4892
4893template <>
4895 const Point<2> &xi,
4896 const unsigned int i)
4897{
4899
4900 const double x = xi[0];
4901 const double y = xi[1];
4902 switch (i)
4903 {
4904 case 0:
4905 return Point<2>(-(1 - y), -(1 - x));
4906 case 1:
4907 return Point<2>(1 - y, -x);
4908 case 2:
4909 return Point<2>(-y, 1 - x);
4910 case 3:
4911 return Point<2>(y, x);
4912 }
4913 return Point<2>(-1e9, -1e9);
4914}
4915
4916
4917
4918template <>
4920 const Point<3> &xi,
4921 const unsigned int i)
4922{
4924
4925 const double x = xi[0];
4926 const double y = xi[1];
4927 const double z = xi[2];
4928 switch (i)
4929 {
4930 case 0:
4931 return Point<3>(-(1 - y) * (1 - z),
4932 -(1 - x) * (1 - z),
4933 -(1 - x) * (1 - y));
4934 case 1:
4935 return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4936 case 2:
4937 return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4938 case 3:
4939 return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4940 case 4:
4941 return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4942 case 5:
4943 return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4944 case 6:
4945 return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4946 case 7:
4947 return Point<3>(y * z, x * z, x * y);
4948 }
4949
4950 return Point<3>(-1e9, -1e9, -1e9);
4951}
4952
4953
4954
4955template <int dim>
4956inline Tensor<1, dim>
4958 const unsigned int)
4959{
4961 return Tensor<1, dim>();
4962}
4963
4964
4965unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4966{
4967 return 0;
4968}
4969
4970
4971namespace internal
4972{
4973 namespace GeometryInfoHelper
4974 {
4975 // wedge product of a single
4976 // vector in 2d: we just have to
4977 // rotate it by 90 degrees to the
4978 // right
4979 inline Tensor<1, 2>
4980 wedge_product(const Tensor<1, 2> (&derivative)[1])
4981 {
4982 Tensor<1, 2> result;
4983 result[0] = derivative[0][1];
4984 result[1] = -derivative[0][0];
4985
4986 return result;
4987 }
4988
4989
4990 // wedge product of 2 vectors in
4991 // 3d is the cross product
4992 inline Tensor<1, 3>
4993 wedge_product(const Tensor<1, 3> (&derivative)[2])
4994 {
4995 return cross_product_3d(derivative[0], derivative[1]);
4996 }
4997
4998
4999 // wedge product of dim vectors
5000 // in dim-d: that's the
5001 // determinant of the matrix
5002 template <int dim>
5003 inline Tensor<0, dim>
5004 wedge_product(const Tensor<1, dim> (&derivative)[dim])
5005 {
5006 Tensor<2, dim> jacobian;
5007 for (unsigned int i = 0; i < dim; ++i)
5008 jacobian[i] = derivative[i];
5009
5010 return determinant(jacobian);
5011 }
5012 } // namespace GeometryInfoHelper
5013} // namespace internal
5014
5015
5016template <int dim>
5017template <int spacedim>
5018inline void
5020# ifndef DEAL_II_CXX14_CONSTEXPR_BUG
5021 (const Point<spacedim> (&vertices)[vertices_per_cell],
5022 Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell])
5023# else
5024 (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms)
5025# endif
5026{
5027 // for each of the vertices,
5028 // compute the alternating form
5029 // of the mapped unit
5030 // vectors. consider for
5031 // example the case of a quad
5032 // in spacedim==3: to do so, we
5033 // need to see how the
5034 // infinitesimal vectors
5035 // (d\xi_1,0) and (0,d\xi_2)
5036 // are transformed into
5037 // spacedim-dimensional space
5038 // and then form their cross
5039 // product (i.e. the wedge product
5040 // of two vectors). to this end, note
5041 // that
5042 // \vec x = sum_i \vec v_i phi_i(\vec xi)
5043 // so the transformed vectors are
5044 // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
5045 // and
5046 // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
5047 // which boils down to the columns
5048 // of the 3x2 matrix \grad_\xi x(\xi)
5049 //
5050 // a similar reasoning would
5051 // hold for all dim,spacedim
5052 // pairs -- we only have to
5053 // compute the wedge product of
5054 // the columns of the
5055 // derivatives
5056 for (unsigned int i = 0; i < vertices_per_cell; ++i)
5057 {
5058 Tensor<1, spacedim> derivatives[dim];
5059
5060 for (unsigned int j = 0; j < vertices_per_cell; ++j)
5061 {
5062 const Tensor<1, dim> grad_phi_j =
5063 d_linear_shape_function_gradient(unit_cell_vertex(i), j);
5064 for (unsigned int l = 0; l < dim; ++l)
5065 derivatives[l] += vertices[j] * grad_phi_j[l];
5066 }
5067
5068 forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
5069 }
5070}
5071
5072#endif // DOXYGEN
5074
5075#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:35
#define DEAL_II_HOST_DEVICE
Definition config.h:166
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#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:238
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 > &)