Loading [MathJax]/extensions/TeX/AMSsymbols.js
 Reference documentation for deal.II version 9.4.1
\(\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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
geometry_info.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_geometry_info_h
17#define dealii_geometry_info_h
18
19
20#include <deal.II/base/config.h>
21
24#include <deal.II/base/point.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 }
76
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, 1, 5, 4, 2, 3, 7, 6}};
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
486template <int dim>
488{
524 {
529
542 };
543};
544
545
546
556template <>
558{
594 {
602 cut_x = 1,
607 };
608};
609
610
611
622template <>
624{
660 {
668 cut_x = 1,
672 cut_y = 2,
676 cut_xy = cut_x | cut_y,
677
681 isotropic_refinement = cut_xy
682 };
683};
684
685
686
697template <>
699{
735 {
743 cut_x = 1,
747 cut_y = 2,
751 cut_xy = cut_x | cut_y,
755 cut_z = 4,
759 cut_xz = cut_x | cut_z,
763 cut_yz = cut_y | cut_z,
767 cut_xyz = cut_x | cut_y | cut_z,
768
772 isotropic_refinement = cut_xyz
773 };
774};
775
776
777
788template <int dim>
790{
791public:
796
802 const typename RefinementPossibilities<dim>::Possibilities refinement_case);
803
809 explicit RefinementCase(const std::uint8_t refinement_case);
810
822 operator std::uint8_t() const;
823
829 operator|(const RefinementCase &r) const;
830
836 operator&(const RefinementCase &r) const;
837
846 operator~() const;
847
848
854 static RefinementCase
855 cut_axis(const unsigned int i);
856
860 static std::size_t
862
868 template <class Archive>
869 void
870 serialize(Archive &ar, const unsigned int version);
871
877 int,
878 << "The refinement flags given (" << arg1
879 << ") contain set bits that do not "
880 << "make sense for the space dimension of the object to which they are applied.");
881
882private:
887 std::uint8_t value : (dim > 0 ? dim : 1);
888};
889
890
891namespace internal
892{
912 template <int dim>
914 {
919 {
924
928 case_isotropic = static_cast<std::uint8_t>(-1)
929 };
930 };
931
932
941 template <>
943 {
950 {
955
960 };
961 };
962
963
964
974 template <>
976 {
983 {
988
993 };
994 };
995
996
997
1008 template <>
1010 {
1018 {
1026 case_x = 1,
1030 case_isotropic = case_x
1032 };
1033
1034
1035
1127 template <>
1129 {
1137 {
1139 case_x = 1,
1140 case_x1y = 2,
1141 case_x2y = 3,
1142 case_x1y2y = 4,
1143 case_y = 5,
1144 case_y1x = 6,
1145 case_y2x = 7,
1146 case_y1x2x = 8,
1147 case_xy = 9,
1148
1149 case_isotropic = case_xy
1151 };
1152
1153
1154
1161 template <int dim>
1163 {
1164 public:
1171 subface_possibility);
1172
1184 operator std::uint8_t() const;
1185
1189 static constexpr std::size_t
1191
1197 int,
1198 << "The subface case given (" << arg1 << ") does not make sense "
1199 << "for the space dimension of the object to which they are applied.");
1200
1201 private:
1206 std::uint8_t value : (dim == 3 ? 4 : 1);
1207 };
1208
1209} // namespace internal
1210
1211
1212
1213template <int dim>
1214struct GeometryInfo;
1215
1216
1217
1237template <>
1239{
1247 static constexpr unsigned int max_children_per_cell = 1;
1248
1252 static constexpr unsigned int faces_per_cell = 0;
1253
1270 static std::array<unsigned int, 0>
1272
1280 static constexpr unsigned int max_children_per_face = 0;
1281
1287 static unsigned int
1288 n_children(const RefinementCase<0> &refinement_case);
1289
1293 static constexpr unsigned int vertices_per_cell = 1;
1294
1313 static std::array<unsigned int, vertices_per_cell>
1315
1339 static unsigned int
1340 face_to_cell_vertices(const unsigned int face,
1341 const unsigned int vertex,
1342 const bool face_orientation = true,
1343 const bool face_flip = false,
1344 const bool face_rotation = false);
1345
1360 static unsigned int
1361 face_to_cell_lines(const unsigned int face,
1362 const unsigned int line,
1363 const bool face_orientation = true,
1364 const bool face_flip = false,
1365 const bool face_rotation = false);
1366
1373 static constexpr unsigned int vertices_per_face = 0;
1374
1378 static constexpr unsigned int lines_per_face = 0;
1379
1383 static constexpr unsigned int quads_per_face = 0;
1384
1388 static constexpr unsigned int lines_per_cell = 0;
1389
1393 static constexpr unsigned int quads_per_cell = 0;
1394
1398 static constexpr unsigned int hexes_per_cell = 0;
1399
1417 static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1418
1432 static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1433};
1434
1435
1436
1967template <int dim>
1969{
1977 static constexpr unsigned int max_children_per_cell = 1 << dim;
1978
1982 static constexpr unsigned int faces_per_cell = 2 * dim;
1983
2002
2010 static constexpr unsigned int max_children_per_face =
2012
2016 static constexpr unsigned int vertices_per_cell = 1 << dim;
2017
2035
2039 static constexpr unsigned int vertices_per_face =
2041
2045 static constexpr unsigned int lines_per_face =
2047
2051 static constexpr unsigned int quads_per_face =
2053
2063 static constexpr unsigned int lines_per_cell =
2064 (2 * GeometryInfo<dim - 1>::lines_per_cell +
2066
2074 static constexpr unsigned int quads_per_cell =
2075 (2 * GeometryInfo<dim - 1>::quads_per_cell +
2077
2081 static constexpr unsigned int hexes_per_cell =
2082 (2 * GeometryInfo<dim - 1>::hexes_per_cell +
2084
2102 static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
2103 internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal();
2104
2118 static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
2119 internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal();
2120
2133 internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face();
2134
2139 static unsigned int
2140 n_children(const RefinementCase<dim> &refinement_case);
2141
2146 static unsigned int
2148
2158 static double
2160 const unsigned int subface_no);
2161
2167 static RefinementCase<dim - 1>
2168 face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2169 const unsigned int face_no,
2170 const bool face_orientation = true,
2171 const bool face_flip = false,
2172 const bool face_rotation = false);
2173
2179 static RefinementCase<dim>
2182 const unsigned int face_no,
2183 const bool face_orientation = true,
2184 const bool face_flip = false,
2185 const bool face_rotation = false);
2186
2191 static RefinementCase<1>
2192 line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2193 const unsigned int line_no);
2194
2199 static RefinementCase<dim>
2201
2248 static unsigned int
2250 const unsigned int face,
2251 const unsigned int subface,
2252 const bool face_orientation = true,
2253 const bool face_flip = false,
2254 const bool face_rotation = false,
2257
2271 static unsigned int
2272 line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
2273
2294 static unsigned int
2295 face_to_cell_vertices(const unsigned int face,
2296 const unsigned int vertex,
2297 const bool face_orientation = true,
2298 const bool face_flip = false,
2299 const bool face_rotation = false);
2300
2312 static unsigned int
2313 face_to_cell_lines(const unsigned int face,
2314 const unsigned int line,
2315 const bool face_orientation = true,
2316 const bool face_flip = false,
2317 const bool face_rotation = false);
2318
2328 static unsigned int
2329 standard_to_real_face_vertex(const unsigned int vertex,
2330 const bool face_orientation = true,
2331 const bool face_flip = false,
2332 const bool face_rotation = false);
2333
2343 static unsigned int
2344 real_to_standard_face_vertex(const unsigned int vertex,
2345 const bool face_orientation = true,
2346 const bool face_flip = false,
2347 const bool face_rotation = false);
2348
2358 static unsigned int
2359 standard_to_real_face_line(const unsigned int line,
2360 const bool face_orientation = true,
2361 const bool face_flip = false,
2362 const bool face_rotation = false);
2363
2369 static unsigned int
2370 standard_to_real_line_vertex(const unsigned int vertex,
2371 const bool line_orientation = true);
2372
2380 static std::array<unsigned int, 2>
2382
2390 static std::array<unsigned int, 2>
2392
2400 static std::array<unsigned int, 2>
2402
2412 static unsigned int
2413 real_to_standard_face_line(const unsigned int line,
2414 const bool face_orientation = true,
2415 const bool face_flip = false,
2416 const bool face_rotation = false);
2417
2423 static Point<dim>
2424 unit_cell_vertex(const unsigned int vertex);
2425
2435 static unsigned int
2437
2445 static Point<dim>
2447 const unsigned int child_index,
2448 const RefinementCase<dim> refine_case =
2450
2456 static Point<dim>
2458 const unsigned int child_index,
2459 const RefinementCase<dim> refine_case =
2461
2466 static bool
2468
2480 static bool
2481 is_inside_unit_cell(const Point<dim> &p, const double eps);
2482
2487 template <typename Number = double>
2488 static Point<dim, Number>
2490
2496 static double
2498
2503 static double
2504 d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2505
2510 static Tensor<1, dim>
2511 d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2512
2564 template <int spacedim>
2565 static void
2567#ifndef DEAL_II_CXX14_CONSTEXPR_BUG
2569 Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2570#else
2571 (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2572#endif
2573
2583 static constexpr std::array<unsigned int, faces_per_cell>
2585 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction();
2586
2603 static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2604 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation();
2605
2616 static constexpr std::array<Tensor<1, dim>, faces_per_cell>
2618 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_vector();
2619
2633 static constexpr ndarray<Tensor<1, dim>, faces_per_cell, dim - 1>
2634 unit_tangential_vectors = internal::GeometryInfoHelper::Initializers<
2636
2642 static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2643 internal::GeometryInfoHelper::Initializers<dim>::opposite_face();
2644
2645
2650 double,
2651 << "The coordinates must satisfy 0 <= x_i <= 1, "
2652 << "but here we have x_i=" << arg1);
2653
2658 int,
2659 int,
2660 int,
2661 << "RefinementCase<dim> " << arg1 << ": face " << arg2
2662 << " has no subface " << arg3);
2663};
2664
2665
2666
2667#ifndef DOXYGEN
2668
2669
2670/* -------------- declaration of explicit specializations ------------- */
2671
2672
2673template <>
2676 const unsigned int i);
2677template <>
2680 const unsigned int i);
2681template <>
2684 const unsigned int i);
2685
2686
2687
2688/* -------------- inline functions ------------- */
2689
2690
2691inline GeometryPrimitive::GeometryPrimitive(const Object object)
2692 : object(object)
2693{}
2694
2695
2696
2697inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2698 : object(static_cast<Object>(object_dimension))
2699{}
2700
2701
2702inline GeometryPrimitive::operator unsigned int() const
2703{
2704 return static_cast<unsigned int>(object);
2705}
2706
2707
2708
2709namespace internal
2710{
2711 template <int dim>
2713 const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2714 : value(subface_possibility)
2715 {}
2716
2717
2718 template <int dim>
2719 inline SubfaceCase<dim>::operator std::uint8_t() const
2720 {
2721 return value;
2722 }
2723
2724
2725} // namespace internal
2726
2727
2728template <int dim>
2730RefinementCase<dim>::cut_axis(const unsigned int)
2731{
2732 Assert(false, ExcInternalError());
2733 return static_cast<std::uint8_t>(-1);
2734}
2735
2736
2737template <>
2738inline RefinementCase<1>
2739RefinementCase<1>::cut_axis(const unsigned int i)
2740{
2741 AssertIndexRange(i, 1);
2742
2744 return options[i];
2745}
2746
2747
2748
2749template <>
2750inline RefinementCase<2>
2751RefinementCase<2>::cut_axis(const unsigned int i)
2752{
2753 AssertIndexRange(i, 2);
2754
2757 return options[i];
2758}
2759
2760
2761
2762template <>
2763inline RefinementCase<3>
2764RefinementCase<3>::cut_axis(const unsigned int i)
2765{
2766 AssertIndexRange(i, 3);
2767
2771 return options[i];
2772}
2773
2774
2775
2776template <int dim>
2778 : value(RefinementPossibilities<dim>::no_refinement)
2779{}
2780
2781
2782
2783template <int dim>
2785 const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2786 : value(refinement_case)
2787{
2788 // check that only those bits of
2789 // the given argument are set that
2790 // make sense for a given space
2791 // dimension
2792 Assert((refinement_case &
2794 refinement_case,
2795 ExcInvalidRefinementCase(refinement_case));
2796}
2797
2798
2799
2800template <int dim>
2801inline RefinementCase<dim>::RefinementCase(const std::uint8_t refinement_case)
2802 : value(refinement_case)
2803{
2804 // check that only those bits of
2805 // the given argument are set that
2806 // make sense for a given space
2807 // dimension
2808 Assert((refinement_case &
2810 refinement_case,
2811 ExcInvalidRefinementCase(refinement_case));
2812}
2813
2814
2815
2816template <int dim>
2817inline RefinementCase<dim>::operator std::uint8_t() const
2818{
2819 return value;
2820}
2821
2822
2823
2824template <int dim>
2827{
2828 return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2829}
2830
2831
2832
2833template <int dim>
2836{
2837 return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2838}
2839
2840
2841
2842template <int dim>
2845{
2846 return RefinementCase<dim>(static_cast<std::uint8_t>(
2848}
2849
2850
2851
2852template <int dim>
2853inline std::size_t
2855{
2856 return sizeof(RefinementCase<dim>);
2857}
2858
2859
2860
2861template <int dim>
2862template <class Archive>
2863inline void
2864RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2865{
2866 // serialization can't deal with bitfields, so copy from/to a full sized
2867 // std::uint8_t
2868 std::uint8_t uchar_value = value;
2869 ar & uchar_value;
2870 value = uchar_value;
2871}
2872
2873
2874
2875template <>
2876inline Point<1>
2877GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2878{
2879 AssertIndexRange(vertex, vertices_per_cell);
2880
2881 return Point<1>(static_cast<double>(vertex));
2882}
2883
2884
2885
2886template <>
2887inline Point<2>
2888GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2889{
2890 AssertIndexRange(vertex, vertices_per_cell);
2891
2892 return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2893}
2894
2895
2896
2897template <>
2898inline Point<3>
2899GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2900{
2901 AssertIndexRange(vertex, vertices_per_cell);
2902
2903 return {static_cast<double>(vertex % 2),
2904 static_cast<double>(vertex / 2 % 2),
2905 static_cast<double>(vertex / 4)};
2906}
2907
2908
2909
2910inline std::array<unsigned int, 0>
2912{
2913 return {{}};
2914}
2915
2916
2917
2918inline std::array<unsigned int, 1>
2920{
2921 return {{0}};
2922}
2923
2924
2925
2926template <int dim>
2929{
2930 return {0U, faces_per_cell};
2931}
2932
2933
2934
2935template <int dim>
2938{
2939 return {0U, vertices_per_cell};
2940}
2941
2942
2943
2944template <int dim>
2945inline Point<dim>
2946GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2947{
2948 Assert(false, ExcNotImplemented());
2949
2950 return {};
2951}
2952
2953
2954
2955template <>
2956inline unsigned int
2958{
2959 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2960
2961 return (p[0] <= 0.5 ? 0 : 1);
2962}
2963
2964
2965
2966template <>
2967inline unsigned int
2969{
2970 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2971 Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2972
2973 return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
2974}
2975
2976
2977
2978template <>
2979inline unsigned int
2981{
2982 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2983 Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2984 Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
2985
2986 return (p[0] <= 0.5 ?
2987 (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
2988 (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
2989}
2990
2991
2992template <int dim>
2993inline unsigned int
2995{
2996 Assert(false, ExcNotImplemented());
2997
2998 return 0;
2999}
3000
3001
3002
3003template <>
3004inline Point<1>
3006 const unsigned int child_index,
3007 const RefinementCase<1> refine_case)
3008
3009{
3010 AssertIndexRange(child_index, 2);
3012 (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3013
3014 return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
3015}
3016
3017
3018
3019template <>
3020inline Point<2>
3022 const unsigned int child_index,
3023 const RefinementCase<2> refine_case)
3024
3025{
3026 AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3027
3028 Point<2> point = p;
3029 switch (refine_case)
3030 {
3032 point[0] *= 2.0;
3033 if (child_index == 1)
3034 point[0] -= 1.0;
3035 break;
3037 point[1] *= 2.0;
3038 if (child_index == 1)
3039 point[1] -= 1.0;
3040 break;
3042 point *= 2.0;
3043 point -= unit_cell_vertex(child_index);
3044 break;
3045 default:
3046 Assert(false, ExcInternalError());
3047 }
3048
3049 return point;
3050}
3051
3052
3053
3054template <>
3055inline Point<3>
3057 const unsigned int child_index,
3058 const RefinementCase<3> refine_case)
3059
3060{
3061 AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3062
3063 Point<3> point = p;
3064 // there might be a cleverer way to do
3065 // this, but since this function is called
3066 // in very few places for initialization
3067 // purposes only, I don't care at the
3068 // moment
3069 switch (refine_case)
3070 {
3072 point[0] *= 2.0;
3073 if (child_index == 1)
3074 point[0] -= 1.0;
3075 break;
3077 point[1] *= 2.0;
3078 if (child_index == 1)
3079 point[1] -= 1.0;
3080 break;
3082 point[2] *= 2.0;
3083 if (child_index == 1)
3084 point[2] -= 1.0;
3085 break;
3087 point[0] *= 2.0;
3088 point[1] *= 2.0;
3089 if (child_index % 2 == 1)
3090 point[0] -= 1.0;
3091 if (child_index / 2 == 1)
3092 point[1] -= 1.0;
3093 break;
3095 // careful, this is slightly
3096 // different from xy and yz due to
3097 // different internal numbering of
3098 // children!
3099 point[0] *= 2.0;
3100 point[2] *= 2.0;
3101 if (child_index / 2 == 1)
3102 point[0] -= 1.0;
3103 if (child_index % 2 == 1)
3104 point[2] -= 1.0;
3105 break;
3107 point[1] *= 2.0;
3108 point[2] *= 2.0;
3109 if (child_index % 2 == 1)
3110 point[1] -= 1.0;
3111 if (child_index / 2 == 1)
3112 point[2] -= 1.0;
3113 break;
3115 point *= 2.0;
3116 point -= unit_cell_vertex(child_index);
3117 break;
3118 default:
3119 Assert(false, ExcInternalError());
3120 }
3121
3122 return point;
3123}
3124
3125
3126
3127template <int dim>
3128inline Point<dim>
3130 const Point<dim> & /*p*/,
3131 const unsigned int /*child_index*/,
3132 const RefinementCase<dim> /*refine_case*/)
3133
3134{
3135 Assert(false, ExcNotImplemented());
3136 return {};
3137}
3138
3139
3140
3141template <>
3142inline Point<1>
3144 const unsigned int child_index,
3145 const RefinementCase<1> refine_case)
3146
3147{
3148 AssertIndexRange(child_index, 2);
3150 (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3151
3152 return (p + unit_cell_vertex(child_index)) * 0.5;
3153}
3154
3155
3156
3157template <>
3158inline Point<3>
3160 const unsigned int child_index,
3161 const RefinementCase<3> refine_case)
3162
3163{
3164 AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3165
3166 Point<3> point = p;
3167 // there might be a cleverer way to do
3168 // this, but since this function is called
3169 // in very few places for initialization
3170 // purposes only, I don't care at the
3171 // moment
3172 switch (refine_case)
3173 {
3175 if (child_index == 1)
3176 point[0] += 1.0;
3177 point[0] *= 0.5;
3178 break;
3180 if (child_index == 1)
3181 point[1] += 1.0;
3182 point[1] *= 0.5;
3183 break;
3185 if (child_index == 1)
3186 point[2] += 1.0;
3187 point[2] *= 0.5;
3188 break;
3190 if (child_index % 2 == 1)
3191 point[0] += 1.0;
3192 if (child_index / 2 == 1)
3193 point[1] += 1.0;
3194 point[0] *= 0.5;
3195 point[1] *= 0.5;
3196 break;
3198 // careful, this is slightly
3199 // different from xy and yz due to
3200 // different internal numbering of
3201 // children!
3202 if (child_index / 2 == 1)
3203 point[0] += 1.0;
3204 if (child_index % 2 == 1)
3205 point[2] += 1.0;
3206 point[0] *= 0.5;
3207 point[2] *= 0.5;
3208 break;
3210 if (child_index % 2 == 1)
3211 point[1] += 1.0;
3212 if (child_index / 2 == 1)
3213 point[2] += 1.0;
3214 point[1] *= 0.5;
3215 point[2] *= 0.5;
3216 break;
3218 point += unit_cell_vertex(child_index);
3219 point *= 0.5;
3220 break;
3221 default:
3222 Assert(false, ExcInternalError());
3223 }
3224
3225 return point;
3226}
3227
3228
3229
3230template <>
3231inline Point<2>
3233 const unsigned int child_index,
3234 const RefinementCase<2> refine_case)
3235{
3236 AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3237
3238 Point<2> point = p;
3239 switch (refine_case)
3240 {
3242 if (child_index == 1)
3243 point[0] += 1.0;
3244 point[0] *= 0.5;
3245 break;
3247 if (child_index == 1)
3248 point[1] += 1.0;
3249 point[1] *= 0.5;
3250 break;
3252 point += unit_cell_vertex(child_index);
3253 point *= 0.5;
3254 break;
3255 default:
3256 Assert(false, ExcInternalError());
3257 }
3258
3259 return point;
3260}
3261
3262
3263
3264template <int dim>
3265inline Point<dim>
3267 const Point<dim> & /*p*/,
3268 const unsigned int /*child_index*/,
3269 const RefinementCase<dim> /*refine_case*/)
3270{
3271 Assert(false, ExcNotImplemented());
3272 return {};
3273}
3274
3275
3276
3277template <int dim>
3278inline bool
3280{
3281 Assert(false, ExcNotImplemented());
3282 return false;
3283}
3284
3285template <>
3286inline bool
3288{
3289 return (p[0] >= 0.) && (p[0] <= 1.);
3290}
3291
3292
3293
3294template <>
3295inline bool
3297{
3298 return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
3299}
3300
3301
3302
3303template <>
3304inline bool
3306{
3307 return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
3308 (p[2] >= 0.) && (p[2] <= 1.);
3309}
3310
3311
3312
3313template <int dim>
3314inline bool
3316{
3317 Assert(false, ExcNotImplemented());
3318 return false;
3319}
3320
3321template <>
3322inline bool
3323GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
3324{
3325 return (p[0] >= -eps) && (p[0] <= 1. + eps);
3326}
3327
3328
3329
3330template <>
3331inline bool
3332GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
3333{
3334 const double l = -eps, u = 1 + eps;
3335 return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
3336}
3337
3338
3339
3340template <>
3341inline bool
3342GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
3343{
3344 const double l = -eps, u = 1.0 + eps;
3345 return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
3346 (p[2] >= l) && (p[2] <= u);
3347}
3348
3349
3350
3351template <>
3352inline unsigned int
3353GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
3354 const unsigned int vertex)
3355{
3356 (void)line;
3357 AssertIndexRange(line, lines_per_cell);
3358 AssertIndexRange(vertex, 2);
3359
3360 return vertex;
3361}
3362
3363
3364template <>
3365inline unsigned int
3366GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
3367 const unsigned int vertex)
3368{
3369 constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
3370 return cell_vertices[line][vertex];
3371}
3372
3373
3374
3375template <>
3376inline unsigned int
3377GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
3378 const unsigned int vertex)
3379{
3380 AssertIndexRange(line, lines_per_cell);
3381 AssertIndexRange(vertex, 2);
3382
3383 constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3384 {1, 3},
3385 {0, 1},
3386 {2, 3},
3387 {4, 6}, // top face
3388 {5, 7},
3389 {4, 5},
3390 {6, 7},
3391 {0,
3392 4}, // connects of bottom
3393 {1, 5}, // top face
3394 {2, 6},
3395 {3, 7}};
3396 return vertices[line][vertex];
3397}
3398
3399
3400template <>
3401inline unsigned int
3402GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3403{
3404 Assert(false, ExcNotImplemented());
3406}
3407
3408template <>
3409inline unsigned int
3410GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3411 const bool face_orientation,
3412 const bool face_flip,
3413 const bool face_rotation)
3414{
3416
3417 // set up a table to make sure that
3418 // we handle non-standard faces correctly
3419 //
3420 // so set up a table that for each vertex (of
3421 // a quad in standard position) describes
3422 // which vertex to take
3423 //
3424 // first index: four vertices 0...3
3425 //
3426 // second index: face_orientation; 0:
3427 // opposite normal, 1: standard
3428 //
3429 // third index: face_flip; 0: standard, 1:
3430 // face rotated by 180 degrees
3431 //
3432 // forth index: face_rotation: 0: standard,
3433 // 1: face rotated by 90 degrees
3434
3435 constexpr unsigned int vertex_translation[4][2][2][2] = {
3436 {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3437 // face_rotation=false and true
3438 {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3439 // face_rotation=false and true
3440 {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3441 // face_rotation=false and true
3442 {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3443 // face_rotation=false and true
3444
3445 {{{2, 3}, // vertex 1 ...
3446 {1, 0}},
3447 {{1, 0}, {2, 3}}},
3448
3449 {{{1, 0}, // vertex 2 ...
3450 {2, 3}},
3451 {{2, 3}, {1, 0}}},
3452
3453 {{{3, 1}, // vertex 3 ...
3454 {0, 2}},
3455 {{3, 1}, {0, 2}}}};
3456
3457 return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3458}
3459
3460
3461
3462template <int dim>
3463inline unsigned int
3464GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3465 const bool,
3466 const bool,
3467 const bool)
3468{
3469 Assert(dim > 1, ExcImpossibleInDim(dim));
3471 return vertex;
3472}
3473
3474template <int dim>
3475inline unsigned int
3477{
3478 constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3479 0, 2, 2, 4, 2, 4, 4, 8};
3480
3481 return n_children[ref_case];
3482}
3483
3484
3485
3486template <int dim>
3487inline unsigned int
3489{
3490 Assert(false, ExcNotImplemented());
3491 return 0;
3492}
3493
3494template <>
3495inline unsigned int
3497{
3498 Assert(false, ExcImpossibleInDim(1));
3499 return 0;
3500}
3501
3502template <>
3503inline unsigned int
3505{
3506 return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3507}
3508
3509
3510
3511template <>
3512inline unsigned int
3514{
3515 const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3516 0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3517 return nsubs[subface_case];
3518}
3519
3520
3521
3522template <int dim>
3523inline double
3525 const unsigned int)
3526{
3527 Assert(false, ExcNotImplemented());
3528 return 0.;
3529}
3530
3531template <>
3532inline double
3534 const unsigned int)
3535{
3536 return 1;
3537}
3538
3539
3540template <>
3541inline double
3543 const unsigned int)
3544{
3545 double ratio = 1;
3546 switch (subface_case)
3547 {
3549 // Here, an
3550 // Assert(false,ExcInternalError())
3551 // would be the right
3552 // choice, but
3553 // unfortunately the
3554 // current function is
3555 // also called for faces
3556 // without children (see
3557 // tests/fe/mapping.cc).
3558 // Assert(false, ExcMessage("Face has no subfaces."));
3559 // Furthermore, assign
3560 // following value as
3561 // otherwise the
3562 // bits/volume_x tests
3563 // break
3565 break;
3567 ratio = 0.5;
3568 break;
3569 default:
3570 // there should be no
3571 // cases left
3572 Assert(false, ExcInternalError());
3573 break;
3574 }
3575
3576 return ratio;
3577}
3578
3579
3580template <>
3581inline double
3583 const unsigned int subface_no)
3584{
3585 double ratio = 1;
3586 switch (subface_case)
3587 {
3589 // Here, an
3590 // Assert(false,ExcInternalError())
3591 // would be the right
3592 // choice, but
3593 // unfortunately the
3594 // current function is
3595 // also called for faces
3596 // without children (see
3597 // tests/bits/mesh_3d_16.cc). Add
3598 // following switch to
3599 // avoid diffs in
3600 // tests/bits/mesh_3d_16
3602 break;
3605 ratio = 0.5;
3606 break;
3610 ratio = 0.25;
3611 break;
3614 if (subface_no < 2)
3615 ratio = 0.25;
3616 else
3617 ratio = 0.5;
3618 break;
3621 if (subface_no == 0)
3622 ratio = 0.5;
3623 else
3624 ratio = 0.25;
3625 break;
3626 default:
3627 // there should be no
3628 // cases left
3629 Assert(false, ExcInternalError());
3630 break;
3631 }
3632
3633 return ratio;
3634}
3635
3636
3637
3638template <int dim>
3640 const RefinementCase<dim> &,
3641 const unsigned int,
3642 const bool,
3643 const bool,
3644 const bool)
3645{
3646 Assert(false, ExcNotImplemented());
3647 return RefinementCase<dim - 1>::no_refinement;
3648}
3649
3650template <>
3652 const RefinementCase<1> &,
3653 const unsigned int,
3654 const bool,
3655 const bool,
3656 const bool)
3657{
3658 Assert(false, ExcImpossibleInDim(1));
3659
3661}
3662
3663
3664template <>
3665inline RefinementCase<1>
3667 const RefinementCase<2> &cell_refinement_case,
3668 const unsigned int face_no,
3669 const bool,
3670 const bool,
3671 const bool)
3672{
3673 const unsigned int dim = 2;
3674 AssertIndexRange(cell_refinement_case,
3677
3678 const RefinementCase<dim - 1>
3681 {RefinementCase<dim - 1>::no_refinement, // no_refinement
3682 RefinementCase<dim - 1>::no_refinement},
3683
3684 {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3685
3686 {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3687
3688 {RefinementCase<dim - 1>::cut_x, // cut_xy
3689 RefinementCase<dim - 1>::cut_x}};
3690
3691 return ref_cases[cell_refinement_case][face_no / 2];
3692}
3693
3694
3695template <>
3696inline RefinementCase<2>
3698 const RefinementCase<3> &cell_refinement_case,
3699 const unsigned int face_no,
3700 const bool face_orientation,
3701 const bool /*face_flip*/,
3702 const bool face_rotation)
3703{
3704 const unsigned int dim = 3;
3705 AssertIndexRange(cell_refinement_case,
3708
3709 const RefinementCase<dim - 1>
3712 {RefinementCase<dim - 1>::no_refinement, // no_refinement
3713 RefinementCase<dim - 1>::no_refinement,
3714 RefinementCase<dim - 1>::no_refinement},
3715
3716 {RefinementCase<dim - 1>::no_refinement, // cut_x
3717 RefinementCase<dim - 1>::cut_y,
3718 RefinementCase<dim - 1>::cut_x},
3719
3720 {RefinementCase<dim - 1>::cut_x, // cut_y
3721 RefinementCase<dim - 1>::no_refinement,
3722 RefinementCase<dim - 1>::cut_y},
3723
3724 {RefinementCase<dim - 1>::cut_x, // cut_xy
3725 RefinementCase<dim - 1>::cut_y,
3726 RefinementCase<dim - 1>::cut_xy},
3727
3728 {RefinementCase<dim - 1>::cut_y, // cut_z
3729 RefinementCase<dim - 1>::cut_x,
3730 RefinementCase<dim - 1>::no_refinement},
3731
3732 {RefinementCase<dim - 1>::cut_y, // cut_xz
3733 RefinementCase<dim - 1>::cut_xy,
3734 RefinementCase<dim - 1>::cut_x},
3735
3736 {RefinementCase<dim - 1>::cut_xy, // cut_yz
3737 RefinementCase<dim - 1>::cut_x,
3738 RefinementCase<dim - 1>::cut_y},
3739
3740 {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3741 RefinementCase<dim - 1>::cut_xy,
3742 RefinementCase<dim - 1>::cut_xy},
3743 };
3744
3745 const RefinementCase<dim - 1> ref_case =
3746 ref_cases[cell_refinement_case][face_no / 2];
3747
3748 const RefinementCase<dim - 1> flip[4] = {
3749 RefinementCase<dim - 1>::no_refinement,
3750 RefinementCase<dim - 1>::cut_y,
3751 RefinementCase<dim - 1>::cut_x,
3752 RefinementCase<dim - 1>::cut_xy};
3753
3754 // correct the ref_case for face_orientation
3755 // and face_rotation. for face_orientation,
3756 // 'true' is the default value whereas for
3757 // face_rotation, 'false' is standard. If
3758 // <tt>face_rotation==face_orientation</tt>,
3759 // then one of them is non-standard and we
3760 // have to swap cut_x and cut_y, otherwise no
3761 // change is necessary. face_flip has no
3762 // influence. however, in order to keep the
3763 // interface consistent with other functions,
3764 // we still include it as an argument to this
3765 // function
3766 return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3767}
3768
3769
3770
3771template <int dim>
3772inline RefinementCase<1>
3774 const unsigned int)
3775{
3776 Assert(false, ExcNotImplemented());
3778}
3779
3780template <>
3781inline RefinementCase<1>
3783 const RefinementCase<1> &cell_refinement_case,
3784 const unsigned int line_no)
3785{
3786 (void)line_no;
3787 const unsigned int dim = 1;
3788 (void)dim;
3789 AssertIndexRange(cell_refinement_case,
3792
3793 return cell_refinement_case;
3794}
3795
3796
3797template <>
3798inline RefinementCase<1>
3800 const RefinementCase<2> &cell_refinement_case,
3801 const unsigned int line_no)
3802{
3803 // Assertions are in face_refinement_case()
3804 return face_refinement_case(cell_refinement_case, line_no);
3805}
3806
3807
3808template <>
3809inline RefinementCase<1>
3811 const RefinementCase<3> &cell_refinement_case,
3812 const unsigned int line_no)
3813{
3814 const unsigned int dim = 3;
3815 AssertIndexRange(cell_refinement_case,
3818
3819 // array indicating, which simple refine
3820 // case cuts a line in direction x, y or
3821 // z. For example, cut_y and everything
3822 // containing cut_y (cut_xy, cut_yz,
3823 // cut_xyz) cuts lines, which are in y
3824 // direction.
3825 const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3828
3829 // order the direction of lines
3830 // 0->x, 1->y, 2->z
3831 const unsigned int direction[lines_per_cell] = {
3832 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3833
3834 return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3837}
3838
3839
3840
3841template <int dim>
3845 const unsigned int,
3846 const bool,
3847 const bool,
3848 const bool)
3849{
3850 Assert(false, ExcNotImplemented());
3851
3853}
3854
3855template <>
3856inline RefinementCase<1>
3858 const RefinementCase<0> &,
3859 const unsigned int,
3860 const bool,
3861 const bool,
3862 const bool)
3863{
3864 const unsigned int dim = 1;
3865 Assert(false, ExcImpossibleInDim(dim));
3866
3868}
3869
3870
3871template <>
3872inline RefinementCase<2>
3874 const RefinementCase<1> &face_refinement_case,
3875 const unsigned int face_no,
3876 const bool,
3877 const bool,
3878 const bool)
3879{
3880 const unsigned int dim = 2;
3881 AssertIndexRange(face_refinement_case,
3884
3885 if (face_refinement_case == RefinementCase<dim>::cut_x)
3886 return (face_no / 2) != 0u ? RefinementCase<dim>::cut_x :
3888 else
3890}
3891
3892
3893template <>
3894inline RefinementCase<3>
3896 const RefinementCase<2> &face_refinement_case,
3897 const unsigned int face_no,
3898 const bool face_orientation,
3899 const bool /*face_flip*/,
3900 const bool face_rotation)
3901{
3902 const unsigned int dim = 3;
3903 AssertIndexRange(face_refinement_case,
3906
3911
3912 // correct the face_refinement_case for
3913 // face_orientation and face_rotation. for
3914 // face_orientation, 'true' is the default
3915 // value whereas for face_rotation, 'false'
3916 // is standard. If
3917 // <tt>face_rotation==face_orientation</tt>,
3918 // then one of them is non-standard and we
3919 // have to swap cut_x and cut_y, otherwise no
3920 // change is necessary. face_flip has no
3921 // influence. however, in order to keep the
3922 // interface consistent with other functions,
3923 // we still include it as an argument to this
3924 // function
3925 const RefinementCase<dim - 1> std_face_ref =
3926 (face_orientation == face_rotation) ? flip[face_refinement_case] :
3927 face_refinement_case;
3928
3929 const RefinementCase<dim> face_to_cell[3][4] = {
3930 {RefinementCase<dim>::no_refinement, // faces 0 and 1
3931 RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3934
3935 {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3936 // "exchanged on faces 2 and 3")
3940
3941 {RefinementCase<dim>::no_refinement, // faces 4 and 5
3945
3946 return face_to_cell[face_no / 2][std_face_ref];
3947}
3948
3949
3950
3951template <int dim>
3954 const unsigned int)
3955{
3956 Assert(false, ExcNotImplemented());
3957
3959}
3960
3961template <>
3962inline RefinementCase<1>
3964 const unsigned int line_no)
3965{
3966 (void)line_no;
3967 AssertIndexRange(line_no, 1);
3968
3970}
3971
3972
3973template <>
3974inline RefinementCase<2>
3976 const unsigned int line_no)
3977{
3978 const unsigned int dim = 2;
3979 (void)dim;
3981
3982 return (line_no / 2) != 0u ? RefinementCase<2>::cut_x :
3984}
3985
3986
3987template <>
3988inline RefinementCase<3>
3990 const unsigned int line_no)
3991{
3992 const unsigned int dim = 3;
3994
3995 const RefinementCase<dim> ref_cases[6] = {
3996 RefinementCase<dim>::cut_y, // lines 0 and 1
3997 RefinementCase<dim>::cut_x, // lines 2 and 3
3998 RefinementCase<dim>::cut_y, // lines 4 and 5
3999 RefinementCase<dim>::cut_x, // lines 6 and 7
4000 RefinementCase<dim>::cut_z, // lines 8 and 9
4001 RefinementCase<dim>::cut_z}; // lines 10 and 11
4002
4003 return ref_cases[line_no / 2];
4004}
4005
4006
4007
4008template <>
4009inline unsigned int
4010GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
4011 const bool face_orientation,
4012 const bool face_flip,
4013 const bool face_rotation)
4014{
4016
4017 // set up a table to make sure that
4018 // we handle non-standard faces correctly
4019 //
4020 // so set up a table that for each vertex (of
4021 // a quad in standard position) describes
4022 // which vertex to take
4023 //
4024 // first index: four vertices 0...3
4025 //
4026 // second index: face_orientation; 0:
4027 // opposite normal, 1: standard
4028 //
4029 // third index: face_flip; 0: standard, 1:
4030 // face rotated by 180 degrees
4031 //
4032 // forth index: face_rotation: 0: standard,
4033 // 1: face rotated by 90 degrees
4034
4035 const unsigned int vertex_translation[4][2][2][2] = {
4036 {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
4037 // face_rotation=false and true
4038 {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
4039 // face_rotation=false and true
4040 {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
4041 // face_rotation=false and true
4042 {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
4043 // face_rotation=false and true
4044
4045 {{{2, 3}, // vertex 1 ...
4046 {1, 0}},
4047 {{1, 3}, {2, 0}}},
4048
4049 {{{1, 0}, // vertex 2 ...
4050 {2, 3}},
4051 {{2, 0}, {1, 3}}},
4052
4053 {{{3, 1}, // vertex 3 ...
4054 {0, 2}},
4055 {{3, 2}, {0, 1}}}};
4056
4057 return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
4058}
4059
4060
4061
4062template <int dim>
4063inline unsigned int
4064GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
4065 const bool,
4066 const bool,
4067 const bool)
4068{
4069 Assert(dim > 1, ExcImpossibleInDim(dim));
4071 return vertex;
4072}
4073
4074
4075
4076template <>
4077inline unsigned int
4078GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
4079 const bool face_orientation,
4080 const bool face_flip,
4081 const bool face_rotation)
4082{
4084
4085
4086 // make sure we handle
4087 // non-standard faces correctly
4088 //
4089 // so set up a table that for each line (of a
4090 // quad) describes which line to take
4091 //
4092 // first index: four lines 0...3
4093 //
4094 // second index: face_orientation; 0:
4095 // opposite normal, 1: standard
4096 //
4097 // third index: face_flip; 0: standard, 1:
4098 // face rotated by 180 degrees
4099 //
4100 // forth index: face_rotation: 0: standard,
4101 // 1: face rotated by 90 degrees
4102
4103 const unsigned int line_translation[4][2][2][2] = {
4104 {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4105 // face_rotation=false and true
4106 {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4107 // face_rotation=false and true
4108 {{0, 3}, // line 0, face_orientation=true, face_flip=false,
4109 // face_rotation=false and true
4110 {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
4111 // face_rotation=false and true
4112
4113 {{{3, 1}, // line 1 ...
4114 {2, 0}},
4115 {{1, 2}, {0, 3}}},
4116
4117 {{{0, 3}, // line 2 ...
4118 {1, 2}},
4119 {{2, 0}, {3, 1}}},
4120
4121 {{{1, 2}, // line 3 ...
4122 {0, 3}},
4123 {{3, 1}, {2, 0}}}};
4124
4125 return line_translation[line][face_orientation][face_flip][face_rotation];
4126}
4127
4128
4129
4130template <int dim>
4131inline unsigned int
4132GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
4133 const bool,
4134 const bool,
4135 const bool)
4136{
4137 Assert(false, ExcNotImplemented());
4138 return line;
4139}
4140
4141
4142
4143template <>
4144inline unsigned int
4145GeometryInfo<2>::standard_to_real_line_vertex(const unsigned int vertex,
4146 const bool line_orientation)
4147{
4148 return line_orientation ? vertex : (1 - vertex);
4149}
4150
4151
4152
4153template <int dim>
4154inline unsigned int
4155GeometryInfo<dim>::standard_to_real_line_vertex(const unsigned int vertex,
4156 const bool)
4157{
4158 Assert(false, ExcNotImplemented());
4159 return vertex;
4160}
4161
4162
4163
4164template <>
4165inline std::array<unsigned int, 2>
4167 const unsigned int vertex)
4168{
4169 return {{vertex % 2, vertex / 2}};
4170}
4171
4172
4173
4174template <int dim>
4175inline std::array<unsigned int, 2>
4177 const unsigned int vertex)
4178{
4179 Assert(false, ExcNotImplemented());
4180 (void)vertex;
4181 return {{0, 0}};
4182}
4183
4184
4185
4186template <>
4187inline std::array<unsigned int, 2>
4189{
4190 // set up a table that for each
4191 // line describes a) from which
4192 // quad to take it, b) which line
4193 // therein it is if the face is
4194 // oriented correctly
4195 static const unsigned int lookup_table[GeometryInfo<3>::lines_per_cell][2] = {
4196 {4, 0}, // take first four lines from bottom face
4197 {4, 1},
4198 {4, 2},
4199 {4, 3},
4200
4201 {5, 0}, // second four lines from top face
4202 {5, 1},
4203 {5, 2},
4204 {5, 3},
4205
4206 {0, 0}, // the rest randomly
4207 {1, 0},
4208 {0, 1},
4209 {1, 1}};
4210
4211 return {{lookup_table[i][0], lookup_table[i][1]}};
4212}
4213
4214
4215
4216template <int dim>
4217inline std::array<unsigned int, 2>
4219{
4220 Assert(false, ExcNotImplemented());
4221 (void)line;
4222 return {{0, 0}};
4223}
4224
4225
4226
4227template <>
4228inline std::array<unsigned int, 2>
4230 const unsigned int vertex)
4231{
4232 // get the corner indices by asking either the bottom or the top face for its
4233 // vertices. handle non-standard faces by calling the vertex reordering
4234 // function from GeometryInfo
4235
4236 // bottom face (4) for first four vertices, top face (5) for the rest
4237 return {{4 + vertex / 4, vertex % 4}};
4238}
4239
4240
4241
4242template <int dim>
4243inline std::array<unsigned int, 2>
4245 const unsigned int vertex)
4246{
4247 Assert(false, ExcNotImplemented());
4248 (void)vertex;
4249 return {{0, 0}};
4250}
4251
4252
4253
4254template <>
4255inline unsigned int
4256GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
4257 const bool face_orientation,
4258 const bool face_flip,
4259 const bool face_rotation)
4260{
4262
4263
4264 // make sure we handle
4265 // non-standard faces correctly
4266 //
4267 // so set up a table that for each line (of a
4268 // quad) describes which line to take
4269 //
4270 // first index: four lines 0...3
4271 //
4272 // second index: face_orientation; 0:
4273 // opposite normal, 1: standard
4274 //
4275 // third index: face_flip; 0: standard, 1:
4276 // face rotated by 180 degrees
4277 //
4278 // forth index: face_rotation: 0: standard,
4279 // 1: face rotated by 90 degrees
4280
4281 const unsigned int line_translation[4][2][2][2] = {
4282 {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4283 // face_rotation=false and true
4284 {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4285 // face_rotation=false and true
4286 {{0, 2}, // line 0, face_orientation=true, face_flip=false,
4287 // face_rotation=false and true
4288 {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
4289 // face_rotation=false and true
4290
4291 {{{3, 1}, // line 1 ...
4292 {2, 0}},
4293 {{1, 3}, {0, 2}}},
4294
4295 {{{0, 3}, // line 2 ...
4296 {1, 2}},
4297 {{2, 1}, {3, 0}}},
4298
4299 {{{1, 2}, // line 3 ...
4300 {0, 3}},
4301 {{3, 0}, {2, 1}}}};
4302
4303 return line_translation[line][face_orientation][face_flip][face_rotation];
4304}
4305
4306
4307
4308template <int dim>
4309inline unsigned int
4310GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
4311 const bool,
4312 const bool,
4313 const bool)
4314{
4315 Assert(false, ExcNotImplemented());
4316 return line;
4317}
4318
4319
4320
4321template <>
4322inline unsigned int
4324 const unsigned int face,
4325 const unsigned int subface,
4326 const bool,
4327 const bool,
4328 const bool,
4329 const RefinementCase<0> &)
4330{
4331 (void)subface;
4332 AssertIndexRange(face, faces_per_cell);
4333 AssertIndexRange(subface, max_children_per_face);
4334
4335 return face;
4336}
4337
4338
4339
4340template <>
4341inline unsigned int
4343 const unsigned int face,
4344 const unsigned int subface,
4345 const bool /*face_orientation*/,
4346 const bool face_flip,
4347 const bool /*face_rotation*/,
4348 const RefinementCase<1> &)
4349{
4350 AssertIndexRange(face, faces_per_cell);
4351 AssertIndexRange(subface, max_children_per_face);
4352
4353 // always return the child adjacent to the specified
4354 // subface. if the face of a cell is not refined, don't
4355 // throw an assertion but deliver the child adjacent to
4356 // the face nevertheless, i.e. deliver the child of
4357 // this cell adjacent to the subface of a possibly
4358 // refined neighbor. this simplifies setting neighbor
4359 // information in execute_refinement.
4360 const unsigned int
4361 subcells[2][RefinementCase<2>::isotropic_refinement][faces_per_cell]
4362 [max_children_per_face] = {
4363 {
4364 // Normal orientation (face_flip = false)
4365 {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
4366 {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
4367 {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
4368 },
4369 {
4370 // Flipped orientation (face_flip = true)
4371 {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
4372 {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
4373 {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
4374 }};
4375
4376 return subcells[face_flip][ref_case - 1][face][subface];
4377}
4378
4379
4380
4381template <>
4382inline unsigned int
4384 const unsigned int face,
4385 const unsigned int subface,
4386 const bool face_orientation,
4387 const bool face_flip,
4388 const bool face_rotation,
4389 const RefinementCase<2> &face_ref_case)
4390{
4391 const unsigned int dim = 3;
4392
4394 ExcMessage("Cell has no children."));
4395 AssertIndexRange(face, faces_per_cell);
4396 if (!(subface == 0 &&
4398 {
4399 AssertIndexRange(subface,
4400 GeometryInfo<dim - 1>::n_children(face_ref_case));
4401 }
4402
4403 // invalid number used for invalid cases,
4404 // e.g. when the children are more refined at
4405 // a given face than the face itself
4406 const unsigned int e = numbers::invalid_unsigned_int;
4407
4408 // the whole process of finding a child cell
4409 // at a given subface considering the
4410 // possibly anisotropic refinement cases of
4411 // the cell and the face as well as
4412 // orientation, flip and rotation of the face
4413 // is quite complicated. thus, we break it
4414 // down into several steps.
4415
4416 // first step: convert the given face refine
4417 // case to a face refine case concerning the
4418 // face in standard orientation (, flip and
4419 // rotation). This only affects cut_x and
4420 // cut_y
4421 const RefinementCase<dim - 1> flip[4] = {
4422 RefinementCase<dim - 1>::no_refinement,
4423 RefinementCase<dim - 1>::cut_y,
4424 RefinementCase<dim - 1>::cut_x,
4425 RefinementCase<dim - 1>::cut_xy};
4426 // for face_orientation, 'true' is the
4427 // default value whereas for face_rotation,
4428 // 'false' is standard. If
4429 // <tt>face_rotation==face_orientation</tt>,
4430 // then one of them is non-standard and we
4431 // have to swap cut_x and cut_y, otherwise no
4432 // change is necessary.
4433 const RefinementCase<dim - 1> std_face_ref =
4434 (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4435
4436 // second step: convert the given subface
4437 // index to the one for a standard face
4438 // respecting face_orientation, face_flip and
4439 // face_rotation
4440
4441 // first index: face_ref_case
4442 // second index: face_orientation
4443 // third index: face_flip
4444 // forth index: face_rotation
4445 // fifth index: subface index
4446 const unsigned int subface_exchange[4][2][2][2][4] = {
4447 // no_refinement (subface 0 stays 0,
4448 // all others are invalid)
4449 {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4450 {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4451 // cut_x (here, if the face is only
4452 // rotated OR only falsely oriented,
4453 // then subface 0 of the non-standard
4454 // face does NOT correspond to one of
4455 // the subfaces of a standard
4456 // face. Thus we indicate the subface
4457 // which is located at the lower left
4458 // corner (the origin of the face's
4459 // local coordinate system) with
4460 // '0'. The rest of this issue is
4461 // taken care of using the above
4462 // conversion to a 'standard face
4463 // refine case')
4464 {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4465 {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4466 // cut_y (the same applies as for
4467 // cut_x)
4468 {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4469 {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4470 // cut_xyz: this information is
4471 // identical to the information
4472 // returned by
4473 // GeometryInfo<3>::real_to_standard_face_vertex()
4474 {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4475 // face_rotation=false, subfaces 0,1,2,3
4476 {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4477 // face_rotation=true, subfaces 0,1,2,3
4478 {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4479 // face_rotation=false, subfaces 0,1,2,3
4480 {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4481 // face_rotation=true, subfaces 0,1,2,3
4482 {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4483 // face_rotation=false, subfaces 0,1,2,3
4484 {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4485 // face_rotation=true, subfaces 0,1,2,3
4486 {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4487 // face_rotation=false, subfaces 0,1,2,3
4488 {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4489 // face_rotation=true, subfaces 0,1,2,3
4490
4491 const unsigned int std_subface =
4492 subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4493 [subface];
4494 Assert(std_subface != e, ExcInternalError());
4495
4496 // third step: these are the children, which
4497 // can be found at the given subfaces of an
4498 // isotropically refined (standard) face
4499 //
4500 // first index: (refinement_case-1)
4501 // second index: face_index
4502 // third index: subface_index (isotropic refinement)
4503 const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4504 [max_children_per_face] = {
4505 // cut_x
4506 {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4507 {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4508 {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4509 {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4510 {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4511 {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4512 // cut_y
4513 {{0, 1, 0, 1},
4514 {0, 1, 0, 1},
4515 {0, 0, 0, 0},
4516 {1, 1, 1, 1},
4517 {0, 0, 1, 1},
4518 {0, 0, 1, 1}},
4519 // cut_xy
4520 {{0, 2, 0, 2},
4521 {1, 3, 1, 3},
4522 {0, 0, 1, 1},
4523 {2, 2, 3, 3},
4524 {0, 1, 2, 3},
4525 {0, 1, 2, 3}},
4526 // cut_z
4527 {{0, 0, 1, 1},
4528 {0, 0, 1, 1},
4529 {0, 1, 0, 1},
4530 {0, 1, 0, 1},
4531 {0, 0, 0, 0},
4532 {1, 1, 1, 1}},
4533 // cut_xz
4534 {{0, 0, 1, 1},
4535 {2, 2, 3, 3},
4536 {0, 1, 2, 3},
4537 {0, 1, 2, 3},
4538 {0, 2, 0, 2},
4539 {1, 3, 1, 3}},
4540 // cut_yz
4541 {{0, 1, 2, 3},
4542 {0, 1, 2, 3},
4543 {0, 2, 0, 2},
4544 {1, 3, 1, 3},
4545 {0, 0, 1, 1},
4546 {2, 2, 3, 3}},
4547 // cut_xyz
4548 {{0, 2, 4, 6},
4549 {1, 3, 5, 7},
4550 {0, 4, 1, 5},
4551 {2, 6, 3, 7},
4552 {0, 1, 2, 3},
4553 {4, 5, 6, 7}}};
4554
4555 // forth step: check, whether the given face
4556 // refine case is valid for the given cell
4557 // refine case. this is the case, if the
4558 // given face refine case is at least as
4559 // refined as the face is for the given cell
4560 // refine case
4561
4562 // note, that we are considering standard
4563 // face refinement cases here and thus must
4564 // not pass the given orientation, flip and
4565 // rotation flags
4566 if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4567 face_refinement_case(ref_case, face))
4568 {
4569 // all is fine. for anisotropic face
4570 // refine cases, select one of the
4571 // isotropic subfaces which neighbors the
4572 // same child
4573
4574 // first index: (standard) face refine case
4575 // second index: subface index
4576 const unsigned int equivalent_iso_subface[4][4] = {
4577 {0, e, e, e}, // no_refinement
4578 {0, 3, e, e}, // cut_x
4579 {0, 3, e, e}, // cut_y
4580 {0, 1, 2, 3}}; // cut_xy
4581
4582 const unsigned int equ_std_subface =
4583 equivalent_iso_subface[std_face_ref][std_subface];
4584 Assert(equ_std_subface != e, ExcInternalError());
4585
4586 return iso_children[ref_case - 1][face][equ_std_subface];
4587 }
4588 else
4589 {
4590 // the face_ref_case was too coarse,
4591 // throw an error
4592 Assert(false,
4593 ExcMessage("The face RefineCase is too coarse "
4594 "for the given cell RefineCase."));
4595 }
4596 // we only get here in case of an error
4597 return e;
4598}
4599
4600
4601
4602template <>
4603inline unsigned int
4605 const unsigned int,
4606 const unsigned int,
4607 const bool,
4608 const bool,
4609 const bool,
4610 const RefinementCase<3> &)
4611{
4612 Assert(false, ExcNotImplemented());
4614}
4615
4616
4617
4618template <>
4619inline unsigned int
4620GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4621 const unsigned int line,
4622 const bool,
4623 const bool,
4624 const bool)
4625{
4626 (void)line;
4627 AssertIndexRange(face, faces_per_cell);
4628 AssertIndexRange(line, lines_per_face);
4629
4630 // The face is a line itself.
4631 return face;
4632}
4633
4634
4635
4636template <>
4637inline unsigned int
4638GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4639 const unsigned int line,
4640 const bool face_orientation,
4641 const bool face_flip,
4642 const bool face_rotation)
4643{
4644 AssertIndexRange(face, faces_per_cell);
4645 AssertIndexRange(line, lines_per_face);
4646
4647 const unsigned lines[faces_per_cell][lines_per_face] = {
4648 {8, 10, 0, 4}, // left face
4649 {9, 11, 1, 5}, // right face
4650 {2, 6, 8, 9}, // front face
4651 {3, 7, 10, 11}, // back face
4652 {0, 1, 2, 3}, // bottom face
4653 {4, 5, 6, 7}}; // top face
4654 return lines[face][real_to_standard_face_line(
4655 line, face_orientation, face_flip, face_rotation)];
4656}
4657
4658
4659
4660inline unsigned int
4661GeometryInfo<0>::face_to_cell_lines(const unsigned int,
4662 const unsigned int,
4663 const bool,
4664 const bool,
4665 const bool)
4666{
4667 Assert(false, ExcNotImplemented());
4669}
4670
4671
4672
4673template <int dim>
4674inline unsigned int
4675GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4676 const unsigned int,
4677 const bool,
4678 const bool,
4679 const bool)
4680{
4681 Assert(false, ExcNotImplemented());
4683}
4684
4685
4686
4687template <int dim>
4688inline unsigned int
4689GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4690 const unsigned int vertex,
4691 const bool face_orientation,
4692 const bool face_flip,
4693 const bool face_rotation)
4694{
4695 return child_cell_on_face(RefinementCase<dim>::isotropic_refinement,
4696 face,
4697 vertex,
4698 face_orientation,
4699 face_flip,
4700 face_rotation);
4701}
4702
4703
4704
4705inline unsigned int
4707 const unsigned int,
4708 const bool,
4709 const bool,
4710 const bool)
4711{
4712 Assert(false, ExcNotImplemented());
4714}
4715
4716
4717
4718template <int dim>
4719template <typename Number>
4720inline Point<dim, Number>
4722{
4724 for (unsigned int i = 0; i < dim; ++i)
4725 p[i] = std::min(std::max(q[i], Number(0.)), Number(1.));
4726
4727 return p;
4728}
4729
4730
4731
4732template <int dim>
4733inline double
4735{
4736 double result = 0.0;
4737
4738 for (unsigned int i = 0; i < dim; ++i)
4739 {
4740 result = std::max(result, -p[i]);
4741 result = std::max(result, p[i] - 1.);
4742 }
4743
4744 return result;
4745}
4746
4747
4748
4749template <int dim>
4750inline double
4752 const unsigned int i)
4753{
4755
4756 switch (dim)
4757 {
4758 case 1:
4759 {
4760 const double x = xi[0];
4761 switch (i)
4762 {
4763 case 0:
4764 return 1 - x;
4765 case 1:
4766 return x;
4767 }
4768 break;
4769 }
4770
4771 case 2:
4772 {
4773 const double x = xi[0];
4774 const double y = xi[1];
4775 switch (i)
4776 {
4777 case 0:
4778 return (1 - x) * (1 - y);
4779 case 1:
4780 return x * (1 - y);
4781 case 2:
4782 return (1 - x) * y;
4783 case 3:
4784 return x * y;
4785 }
4786 break;
4787 }
4788
4789 case 3:
4790 {
4791 const double x = xi[0];
4792 const double y = xi[1];
4793 const double z = xi[2];
4794 switch (i)
4795 {
4796 case 0:
4797 return (1 - x) * (1 - y) * (1 - z);
4798 case 1:
4799 return x * (1 - y) * (1 - z);
4800 case 2:
4801 return (1 - x) * y * (1 - z);
4802 case 3:
4803 return x * y * (1 - z);
4804 case 4:
4805 return (1 - x) * (1 - y) * z;
4806 case 5:
4807 return x * (1 - y) * z;
4808 case 6:
4809 return (1 - x) * y * z;
4810 case 7:
4811 return x * y * z;
4812 }
4813 break;
4814 }
4815
4816 default:
4817 Assert(false, ExcNotImplemented());
4818 }
4819 return -1e9;
4820}
4821
4822
4823
4824template <>
4826 const Point<1> &,
4827 const unsigned int i)
4828{
4830
4831 switch (i)
4832 {
4833 case 0:
4834 return Point<1>(-1.);
4835 case 1:
4836 return Point<1>(1.);
4837 }
4838
4839 return Point<1>(-1e9);
4840}
4841
4842
4843
4844template <>
4846 const Point<2> & xi,
4847 const unsigned int i)
4848{
4850
4851 const double x = xi[0];
4852 const double y = xi[1];
4853 switch (i)
4854 {
4855 case 0:
4856 return Point<2>(-(1 - y), -(1 - x));
4857 case 1:
4858 return Point<2>(1 - y, -x);
4859 case 2:
4860 return Point<2>(-y, 1 - x);
4861 case 3:
4862 return Point<2>(y, x);
4863 }
4864 return Point<2>(-1e9, -1e9);
4865}
4866
4867
4868
4869template <>
4871 const Point<3> & xi,
4872 const unsigned int i)
4873{
4875
4876 const double x = xi[0];
4877 const double y = xi[1];
4878 const double z = xi[2];
4879 switch (i)
4880 {
4881 case 0:
4882 return Point<3>(-(1 - y) * (1 - z),
4883 -(1 - x) * (1 - z),
4884 -(1 - x) * (1 - y));
4885 case 1:
4886 return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4887 case 2:
4888 return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4889 case 3:
4890 return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4891 case 4:
4892 return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4893 case 5:
4894 return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4895 case 6:
4896 return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4897 case 7:
4898 return Point<3>(y * z, x * z, x * y);
4899 }
4900
4901 return Point<3>(-1e9, -1e9, -1e9);
4902}
4903
4904
4905
4906template <int dim>
4907inline Tensor<1, dim>
4909 const unsigned int)
4910{
4911 Assert(false, ExcNotImplemented());
4912 return Tensor<1, dim>();
4913}
4914
4915
4916unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4917{
4918 return 0;
4919}
4920
4921
4922namespace internal
4923{
4924 namespace GeometryInfoHelper
4925 {
4926 // wedge product of a single
4927 // vector in 2d: we just have to
4928 // rotate it by 90 degrees to the
4929 // right
4930 inline Tensor<1, 2>
4931 wedge_product(const Tensor<1, 2> (&derivative)[1])
4932 {
4933 Tensor<1, 2> result;
4934 result[0] = derivative[0][1];
4935 result[1] = -derivative[0][0];
4936
4937 return result;
4938 }
4939
4940
4941 // wedge product of 2 vectors in
4942 // 3d is the cross product
4943 inline Tensor<1, 3>
4944 wedge_product(const Tensor<1, 3> (&derivative)[2])
4945 {
4946 return cross_product_3d(derivative[0], derivative[1]);
4947 }
4948
4949
4950 // wedge product of dim vectors
4951 // in dim-d: that's the
4952 // determinant of the matrix
4953 template <int dim>
4954 inline Tensor<0, dim>
4955 wedge_product(const Tensor<1, dim> (&derivative)[dim])
4956 {
4957 Tensor<2, dim> jacobian;
4958 for (unsigned int i = 0; i < dim; ++i)
4959 jacobian[i] = derivative[i];
4960
4961 return determinant(jacobian);
4962 }
4963 } // namespace GeometryInfoHelper
4964} // namespace internal
4965
4966
4967template <int dim>
4968template <int spacedim>
4969inline void
4971# ifndef DEAL_II_CXX14_CONSTEXPR_BUG
4972 (const Point<spacedim> (&vertices)[vertices_per_cell],
4973 Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell])
4974# else
4975 (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms)
4976# endif
4977{
4978 // for each of the vertices,
4979 // compute the alternating form
4980 // of the mapped unit
4981 // vectors. consider for
4982 // example the case of a quad
4983 // in spacedim==3: to do so, we
4984 // need to see how the
4985 // infinitesimal vectors
4986 // (d\xi_1,0) and (0,d\xi_2)
4987 // are transformed into
4988 // spacedim-dimensional space
4989 // and then form their cross
4990 // product (i.e. the wedge product
4991 // of two vectors). to this end, note
4992 // that
4993 // \vec x = sum_i \vec v_i phi_i(\vec xi)
4994 // so the transformed vectors are
4995 // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
4996 // and
4997 // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
4998 // which boils down to the columns
4999 // of the 3x2 matrix \grad_\xi x(\xi)
5000 //
5001 // a similar reasoning would
5002 // hold for all dim,spacedim
5003 // pairs -- we only have to
5004 // compute the wedge product of
5005 // the columns of the
5006 // derivatives
5007 for (unsigned int i = 0; i < vertices_per_cell; ++i)
5008 {
5009 Tensor<1, spacedim> derivatives[dim];
5010
5011 for (unsigned int j = 0; j < vertices_per_cell; ++j)
5012 {
5013 const Tensor<1, dim> grad_phi_j =
5014 d_linear_shape_function_gradient(unit_cell_vertex(i), j);
5015 for (unsigned int l = 0; l < dim; ++l)
5016 derivatives[l] += vertices[j] * grad_phi_j[l];
5017 }
5018
5019 forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
5020 }
5021}
5022
5023#endif // DOXYGEN
5025
5026#endif
GeometryPrimitive(const unsigned int object_dimension)
GeometryPrimitive(const Object object)
Definition: point.h:111
RefinementCase operator|(const RefinementCase &r) const
static std::size_t memory_consumption()
RefinementCase operator~() const
RefinementCase operator&(const RefinementCase &r) const
RefinementCase(const typename RefinementPossibilities< dim >::Possibilities refinement_case)
static RefinementCase cut_axis(const unsigned int i)
std::uint8_t value
void serialize(Archive &ar, const unsigned int version)
RefinementCase(const std::uint8_t refinement_case)
Definition: tensor.h:503
static constexpr std::size_t memory_consumption()
SubfaceCase(const typename SubfacePossibilities< dim >::Possibilities subface_possibility)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Point< 3 > vertices[4]
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
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)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcInternalError()
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:555
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
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:190
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static const unsigned int invalid_unsigned_int
Definition: types.h:201
boost::integer_range< IncrementableType > iota_view
Definition: iota_view.h:46
::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:108
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 const std::array< unsigned int, vertices_per_cell > dx_to_deal
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 const std::array< unsigned int, vertices_per_cell > ucd_to_deal
static unsigned int child_cell_from_point(const Point< dim > &p)
static unsigned int standard_to_real_line_vertex(const unsigned int vertex, const bool line_orientation=true)
static constexpr std::array< unsigned int, vertices_per_cell > ucd_to_deal
static constexpr ndarray< Tensor< 1, dim >, faces_per_cell, dim - 1 > unit_tangential_vectors
static bool is_inside_unit_cell(const Point< dim > &p, const double eps)
static double distance_to_unit_cell(const Point< dim > &p)
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)
static constexpr unsigned int quads_per_face
static RefinementCase< dim - 1 > face_refinement_case(const RefinementCase< dim > &cell_refinement_case, const unsigned int face_no, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static unsigned int standard_to_real_face_vertex(const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static RefinementCase< dim > min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
static constexpr std::array< unsigned int, faces_per_cell > opposite_face
static constexpr unsigned int max_children_per_face
static unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static constexpr unsigned int vertices_per_cell
static constexpr std::array< unsigned int, vertices_per_cell > dx_to_deal
static constexpr unsigned int lines_per_cell
static constexpr std::array< unsigned int, faces_per_cell > unit_normal_direction
static std::array< unsigned int, 2 > standard_hex_line_to_quad_line_index(const unsigned int line)
static constexpr unsigned int faces_per_cell
static unsigned int standard_to_real_face_line(const unsigned int line, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
static constexpr unsigned int hexes_per_cell
static Point< dim, Number > project_to_unit_cell(const Point< dim, Number > &p)
static RefinementCase< 1 > line_refinement_case(const RefinementCase< dim > &cell_refinement_case, const unsigned int line_no)
static std::array< unsigned int, 2 > standard_hex_vertex_to_quad_vertex_index(const unsigned int vertex)
static std::array< unsigned int, 2 > standard_quad_vertex_to_line_vertex_index(const unsigned int vertex)
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
static constexpr std::array< Tensor< 1, dim >, faces_per_cell > unit_normal_vector
static constexpr unsigned int vertices_per_face
static unsigned int real_to_standard_face_line(const unsigned int line, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
static bool is_inside_unit_cell(const Point< dim > &p)
static RefinementCase< dim > min_cell_refinement_case_for_face_refinement(const RefinementCase< dim - 1 > &face_refinement_case, const unsigned int face_no, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
static unsigned int real_to_standard_face_vertex(const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static unsigned int face_to_cell_lines(const unsigned int face, const unsigned int line, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static Point< dim > cell_to_child_coordinates(const Point< dim > &p, const unsigned int child_index, const RefinementCase< dim > refine_case=RefinementCase< dim >::isotropic_refinement)
static constexpr ndarray< unsigned int, vertices_per_cell, dim > vertex_to_face
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
static void alternating_form_at_vertices(const Point< spacedim >(&vertices)[vertices_per_cell], Tensor< spacedim - dim, spacedim >(&forms)[vertices_per_cell])
static Point< dim > unit_cell_vertex(const unsigned int vertex)
static constexpr std::array< int, faces_per_cell > unit_normal_orientation
static constexpr unsigned int quads_per_cell
static Point< dim > child_to_cell_coordinates(const Point< dim > &p, const unsigned int child_index, const RefinementCase< dim > refine_case=RefinementCase< dim >::isotropic_refinement)
static constexpr unsigned int max_children_per_cell
static unsigned int n_subfaces(const internal::SubfaceCase< dim > &subface_case)
static constexpr unsigned int lines_per_face
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2660