deal.II version GIT relicensing-2289-g1e5549a87a 2024-12-21 21:30:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
geometry_info.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1998 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_geometry_info_h
16#define dealii_geometry_info_h
17
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/base/exceptions.h>
22#include <deal.II/base/ndarray.h>
23#include <deal.II/base/point.h>
24#include <deal.II/base/std_cxx20/iota_view.h>
25
26#include <array>
27#include <cstdint>
28
29
30
31DEAL_II_NAMESPACE_OPEN
32
33#ifndef DOXYGEN
34namespace internal
35{
36 namespace GeometryInfoHelper
37 {
38 // A struct that holds the values for all the arrays we want to initialize
39 // in GeometryInfo
40 template <int dim>
41 struct Initializers;
42
43 template <>
44 struct Initializers<1>
45 {
46 static constexpr std::array<unsigned int, 2>
47 ucd_to_deal()
48 {
49 return {{0, 1}};
50 }
51
52 static constexpr std::array<unsigned int, 2>
53 unit_normal_direction()
54 {
55 return {{0, 0}};
56 }
57
58 static constexpr std::array<int, 2>
59 unit_normal_orientation()
60 {
61 return {{-1, 1}};
62 }
63
64 static constexpr std::array<Tensor<1, 1>, 2>
65 unit_normal_vector()
66 {
67 return {{Tensor<1, 1>{{-1}}, Tensor<1, 1>{{1}}}};
68 }
69
70 static constexpr ::ndarray<Tensor<1, 1>, 2, 0>
71 unit_tangential_vectors()
72 {
73 return {{{{}}, {{}}}};
74 }
75
76 static constexpr std::array<unsigned int, 2>
77 opposite_face()
78 {
79 return {{1, 0}};
80 }
81
82 static constexpr std::array<unsigned int, 2>
83 dx_to_deal()
84 {
85 return {{0, 1}};
86 }
87
88 static constexpr ::ndarray<unsigned int, 2, 1>
89 vertex_to_face()
90 {
91 return {{{{0}}, {{1}}}};
92 }
93 };
94
95 template <>
96 struct Initializers<2>
97 {
98 static constexpr std::array<unsigned int, 4>
99 ucd_to_deal()
100 {
101 return {{0, 1, 3, 2}};
102 }
103
104 static constexpr std::array<unsigned int, 4>
105 unit_normal_direction()
106 {
107 return {{0, 0, 1, 1}};
108 }
109
110 static constexpr std::array<int, 4>
111 unit_normal_orientation()
112 {
113 return {{-1, 1, -1, 1}};
114 }
115
116 static constexpr std::array<Tensor<1, 2>, 4>
117 unit_normal_vector()
118 {
119 return {{Tensor<1, 2>{{-1., 0.}},
120 Tensor<1, 2>{{1., 0.}},
121 Tensor<1, 2>{{0., -1.}},
122 Tensor<1, 2>{{0., 1.}}}};
123 }
124
125 static constexpr ::ndarray<Tensor<1, 2>, 4, 1>
126 unit_tangential_vectors()
127 {
128 return {{{{Tensor<1, 2>{{0, -1}}}},
129 {{Tensor<1, 2>{{0, 1}}}},
130 {{Tensor<1, 2>{{1, 0}}}},
131 {{Tensor<1, 2>{{-1, 0}}}}}};
132 }
133
134 static constexpr std::array<unsigned int, 4>
135 opposite_face()
136 {
137 return {{1, 0, 3, 2}};
138 }
139
140 static constexpr std::array<unsigned int, 4>
141 dx_to_deal()
142 {
143 return {{0, 2, 1, 3}};
144 }
145
146 static constexpr ::ndarray<unsigned int, 4, 2>
147 vertex_to_face()
148 {
149 return {{{{0, 2}}, {{1, 2}}, {{0, 3}}, {{1, 3}}}};
150 }
151 };
152
153 template <>
154 struct Initializers<3>
155 {
156 static constexpr std::array<unsigned int, 8>
157 ucd_to_deal()
158 {
159 return {{0, 4, 5, 1, 2, 6, 7, 3}};
160 }
161
162 static constexpr std::array<unsigned int, 6>
163 unit_normal_direction()
164 {
165 return {{0, 0, 1, 1, 2, 2}};
166 }
167
168 static constexpr std::array<int, 6>
169 unit_normal_orientation()
170 {
171 return {{-1, 1, -1, 1, -1, 1}};
172 }
173
174 static constexpr std::array<Tensor<1, 3>, 6>
175 unit_normal_vector()
176 {
177 return {{Tensor<1, 3>{{-1, 0, 0}},
178 Tensor<1, 3>{{1, 0, 0}},
179 Tensor<1, 3>{{0, -1, 0}},
180 Tensor<1, 3>{{0, 1, 0}},
181 Tensor<1, 3>{{0, 0, -1}},
182 Tensor<1, 3>{{0, 0, 1}}}};
183 }
184
185 static constexpr ::ndarray<Tensor<1, 3>, 6, 2>
186 unit_tangential_vectors()
187 {
188 return {{{{Tensor<1, 3>{{0, -1, 0}}, Tensor<1, 3>{{0, 0, 1}}}},
189 {{Tensor<1, 3>{{0, 1, 0}}, Tensor<1, 3>{{0, 0, 1}}}},
190 {{Tensor<1, 3>{{0, 0, -1}}, Tensor<1, 3>{{1, 0, 0}}}},
191 {{Tensor<1, 3>{{0, 0, 1}}, Tensor<1, 3>{{1, 0, 0}}}},
192 {{Tensor<1, 3>{{-1, 0, 0}}, Tensor<1, 3>{{0, 1, 0}}}},
193 {{Tensor<1, 3>{{1, 0, 0}}, Tensor<1, 3>{{0, 1, 0}}}}}};
194 }
195
196 static constexpr std::array<unsigned int, 6>
197 opposite_face()
198 {
199 return {{1, 0, 3, 2, 5, 4}};
200 }
201
202 static constexpr std::array<unsigned int, 8>
203 dx_to_deal()
204 {
205 return {{0, 4, 2, 6, 1, 5, 3, 7}};
206 }
207
208 static constexpr ::ndarray<unsigned int, 8, 3>
209 vertex_to_face()
210 {
211 return {{{{0, 2, 4}},
212 {{1, 2, 4}},
213 {{0, 3, 4}},
214 {{1, 3, 4}},
215 {{0, 2, 5}},
216 {{1, 2, 5}},
217 {{0, 3, 5}},
218 {{1, 3, 5}}}};
219 }
220 };
221
222 template <>
223 struct Initializers<4>
224 {
225 static constexpr std::array<unsigned int, 16>
226 ucd_to_deal()
227 {
228 return {{numbers::invalid_unsigned_int,
229 numbers::invalid_unsigned_int,
230 numbers::invalid_unsigned_int,
231 numbers::invalid_unsigned_int,
232 numbers::invalid_unsigned_int,
233 numbers::invalid_unsigned_int,
234 numbers::invalid_unsigned_int,
235 numbers::invalid_unsigned_int,
236 numbers::invalid_unsigned_int,
237 numbers::invalid_unsigned_int,
238 numbers::invalid_unsigned_int,
239 numbers::invalid_unsigned_int,
240 numbers::invalid_unsigned_int,
241 numbers::invalid_unsigned_int,
242 numbers::invalid_unsigned_int,
243 numbers::invalid_unsigned_int}};
244 }
245
246 static constexpr std::array<unsigned int, 8>
247 unit_normal_direction()
248 {
249 return {{0, 0, 1, 1, 2, 2, 3, 3}};
250 }
251
252 static constexpr std::array<int, 8>
253 unit_normal_orientation()
254 {
255 return {{-1, 1, -1, 1, -1, 1, -1, 1}};
256 }
257
258 static constexpr std::array<Tensor<1, 4>, 8>
259 unit_normal_vector()
260 {
261 return {{Tensor<1, 4>{{-1, 0, 0, 0}},
262 Tensor<1, 4>{{1, 0, 0, 0}},
263 Tensor<1, 4>{{0, -1, 0, 0}},
264 Tensor<1, 4>{{0, 1, 0, 0}},
265 Tensor<1, 4>{{0, 0, -1, 0}},
266 Tensor<1, 4>{{0, 0, 1, 0}},
267 Tensor<1, 4>{{0, 0, 0, -1}},
268 Tensor<1, 4>{{0, 0, 0, 1}}}};
269 }
270
271 static constexpr ::ndarray<Tensor<1, 4>, 8, 3>
272 unit_tangential_vectors()
273 {
274 return {{{{Tensor<1, 4>{{0, -1, 0, 0}},
275 Tensor<1, 4>{{0, 0, 1, 0}},
276 Tensor<1, 4>{{0, 0, 0, 1}}}},
277 {{Tensor<1, 4>{{0, 1, 0, 0}},
278 Tensor<1, 4>{{0, 0, 1, 0}},
279 Tensor<1, 4>{{0, 0, 0, 1}}}},
280 {{Tensor<1, 4>{{0, 0, -1, 0}},
281 Tensor<1, 4>{{0, 0, 0, 1}},
282 Tensor<1, 4>{{1, 0, 0, 0}}}},
283 {{Tensor<1, 4>{{0, 0, 1, 0}},
284 Tensor<1, 4>{{0, 0, 0, 1}},
285 Tensor<1, 4>{{1, 0, 0, 0}}}},
286 {{Tensor<1, 4>{{0, 0, 0, -1}},
287 Tensor<1, 4>{{1, 0, 0, 0}},
288 Tensor<1, 4>{{0, 1, 0, 0}}}},
289 {{Tensor<1, 4>{{0, 0, 0, 1}},
290 Tensor<1, 4>{{1, 0, 0, 0}},
291 Tensor<1, 4>{{0, 1, 0, 0}}}},
292 {{Tensor<1, 4>{{-1, 0, 0, 0}},
293 Tensor<1, 4>{{0, 1, 0, 0}},
294 Tensor<1, 4>{{0, 0, 1, 0}}}},
295 {{Tensor<1, 4>{{1, 0, 0, 0}},
296 Tensor<1, 4>{{0, 1, 0, 0}},
297 Tensor<1, 4>{{0, 0, 1, 0}}}}}};
298 }
299
300 static constexpr std::array<unsigned int, 8>
301 opposite_face()
302 {
303 return {{1, 0, 3, 2, 5, 4, 7, 6}};
304 }
305
306 static constexpr std::array<unsigned int, 16>
307 dx_to_deal()
308 {
309 return {{numbers::invalid_unsigned_int,
310 numbers::invalid_unsigned_int,
311 numbers::invalid_unsigned_int,
312 numbers::invalid_unsigned_int,
313 numbers::invalid_unsigned_int,
314 numbers::invalid_unsigned_int,
315 numbers::invalid_unsigned_int,
316 numbers::invalid_unsigned_int,
317 numbers::invalid_unsigned_int,
318 numbers::invalid_unsigned_int,
319 numbers::invalid_unsigned_int,
320 numbers::invalid_unsigned_int,
321 numbers::invalid_unsigned_int,
322 numbers::invalid_unsigned_int,
323 numbers::invalid_unsigned_int,
324 numbers::invalid_unsigned_int}};
325 }
326
327 static constexpr ::ndarray<unsigned int, 16, 4>
328 vertex_to_face()
329 {
330 return {{{{numbers::invalid_unsigned_int,
331 numbers::invalid_unsigned_int,
332 numbers::invalid_unsigned_int,
333 numbers::invalid_unsigned_int}},
334 {{numbers::invalid_unsigned_int,
335 numbers::invalid_unsigned_int,
336 numbers::invalid_unsigned_int,
337 numbers::invalid_unsigned_int}},
338 {{numbers::invalid_unsigned_int,
339 numbers::invalid_unsigned_int,
340 numbers::invalid_unsigned_int,
341 numbers::invalid_unsigned_int}},
342 {{numbers::invalid_unsigned_int,
343 numbers::invalid_unsigned_int,
344 numbers::invalid_unsigned_int,
345 numbers::invalid_unsigned_int}},
346 {{numbers::invalid_unsigned_int,
347 numbers::invalid_unsigned_int,
348 numbers::invalid_unsigned_int,
349 numbers::invalid_unsigned_int}},
350 {{numbers::invalid_unsigned_int,
351 numbers::invalid_unsigned_int,
352 numbers::invalid_unsigned_int,
353 numbers::invalid_unsigned_int}},
354 {{numbers::invalid_unsigned_int,
355 numbers::invalid_unsigned_int,
356 numbers::invalid_unsigned_int,
357 numbers::invalid_unsigned_int}},
358 {{numbers::invalid_unsigned_int,
359 numbers::invalid_unsigned_int,
360 numbers::invalid_unsigned_int,
361 numbers::invalid_unsigned_int}},
362 {{numbers::invalid_unsigned_int,
363 numbers::invalid_unsigned_int,
364 numbers::invalid_unsigned_int,
365 numbers::invalid_unsigned_int}},
366 {{numbers::invalid_unsigned_int,
367 numbers::invalid_unsigned_int,
368 numbers::invalid_unsigned_int,
369 numbers::invalid_unsigned_int}},
370 {{numbers::invalid_unsigned_int,
371 numbers::invalid_unsigned_int,
372 numbers::invalid_unsigned_int,
373 numbers::invalid_unsigned_int}},
374 {{numbers::invalid_unsigned_int,
375 numbers::invalid_unsigned_int,
376 numbers::invalid_unsigned_int,
377 numbers::invalid_unsigned_int}},
378 {{numbers::invalid_unsigned_int,
379 numbers::invalid_unsigned_int,
380 numbers::invalid_unsigned_int,
381 numbers::invalid_unsigned_int}},
382 {{numbers::invalid_unsigned_int,
383 numbers::invalid_unsigned_int,
384 numbers::invalid_unsigned_int,
385 numbers::invalid_unsigned_int}},
386 {{numbers::invalid_unsigned_int,
387 numbers::invalid_unsigned_int,
388 numbers::invalid_unsigned_int,
389 numbers::invalid_unsigned_int}},
390 {{numbers::invalid_unsigned_int,
391 numbers::invalid_unsigned_int,
392 numbers::invalid_unsigned_int,
393 numbers::invalid_unsigned_int}}}};
394 }
395 };
396 } // namespace GeometryInfoHelper
397} // namespace internal
398#endif // DOXYGEN
399
400
416class GeometryPrimitive
417{
418public:
425 enum Object
426 {
430 vertex = 0,
434 line = 1,
438 quad = 2,
442 hex = 3
443 };
444
449 GeometryPrimitive(const Object object);
450
456 GeometryPrimitive(const unsigned int object_dimension);
457
462 operator unsigned int() const;
463
464private:
468 Object object;
469};
470
471
472
483template <int dim>
484struct RefinementPossibilities
485{
520 enum Possibilities : std::uint8_t
521 {
525 no_refinement = 0,
526
538 isotropic_refinement = 0xFF
539 };
540};
541
542
543
551template <>
552struct RefinementPossibilities<1>
553{
588 enum Possibilities : std::uint8_t
589 {
593 no_refinement = 0,
597 cut_x = 1,
601 isotropic_refinement = cut_x
602 };
603};
604
605
606
615template <>
616struct RefinementPossibilities<2>
617{
652 enum Possibilities : std::uint8_t
653 {
657 no_refinement = 0,
661 cut_x = 1,
665 cut_y = 2,
669 cut_xy = cut_x | cut_y,
670
674 isotropic_refinement = cut_xy
675 };
676};
677
678
679
688template <>
689struct RefinementPossibilities<3>
690{
725 enum Possibilities : std::uint8_t
726 {
730 no_refinement = 0,
734 cut_x = 1,
738 cut_y = 2,
742 cut_xy = cut_x | cut_y,
746 cut_z = 4,
750 cut_xz = cut_x | cut_z,
754 cut_yz = cut_y | cut_z,
758 cut_xyz = cut_x | cut_y | cut_z,
759
763 isotropic_refinement = cut_xyz
764 };
765};
766
767
768
777template <int dim>
778class RefinementCase : public RefinementPossibilities<dim>
779{
780public:
784 static constexpr unsigned int n_refinement_cases = (1 << dim);
785
789 RefinementCase();
790
795 RefinementCase(
796 const typename RefinementPossibilities<dim>::Possibilities refinement_case);
797
803 explicit RefinementCase(const std::uint8_t refinement_case);
804
816 DEAL_II_HOST_DEVICE operator std::uint8_t() const;
817
822 RefinementCase
823 operator|(const RefinementCase &r) const;
824
829 RefinementCase
830 operator&(const RefinementCase &r) const;
831
839 RefinementCase
840 operator~() const;
841
842
848 static RefinementCase
849 cut_axis(const unsigned int i);
850
859 static std::array<RefinementCase<dim>, n_refinement_cases>
860 all_refinement_cases();
861
865 static std::size_t
866 memory_consumption();
867
873 template <class Archive>
874 void
875 serialize(Archive &ar, const unsigned int version);
876
880 DeclException1(
881 ExcInvalidRefinementCase,
882 int,
883 << "The refinement flags given (" << arg1
884 << ") contain set bits that do not "
885 << "make sense for the space dimension of the object to which they are applied.");
886
887private:
892 std::uint8_t value : (dim > 0 ? dim : 1);
893};
894
895
896
897namespace internal
898{
916 template <int dim>
917 struct SubfacePossibilities
918 {
922 enum Possibilities
923 {
927 case_none = 0,
928
932 case_isotropic = static_cast<std::uint8_t>(-1)
933 };
934 };
935
936
943 template <>
944 struct SubfacePossibilities<0>
945 {
963 };
964
965
966
974 template <>
975 struct SubfacePossibilities<1>
976 {
994 };
995
996
997
1006 template <>
1007 struct SubfacePossibilities<2>
1008 {
1015 enum Possibilities
1016 {
1020 case_none = 0,
1024 case_x = 1,
1028 case_isotropic = case_x
1029 };
1030 };
1031
1032
1033
1123 template <>
1124 struct SubfacePossibilities<3>
1125 {
1132 enum Possibilities
1133 {
1134 case_none = 0,
1135 case_x = 1,
1136 case_x1y = 2,
1137 case_x2y = 3,
1138 case_x1y2y = 4,
1139 case_y = 5,
1140 case_y1x = 6,
1141 case_y2x = 7,
1142 case_y1x2x = 8,
1143 case_xy = 9,
1144
1145 case_isotropic = case_xy
1146 };
1147 };
1148
1149
1150
1155 template <int dim>
1156 class SubfaceCase : public SubfacePossibilities<dim>
1157 {
1158 public:
1164 SubfaceCase(const typename SubfacePossibilities<dim>::Possibilities
1165 subface_possibility);
1166
1178 operator std::uint8_t() const;
1179
1183 static constexpr std::size_t
1184 memory_consumption();
1185
1189 DeclException1(
1190 ExcInvalidSubfaceCase,
1191 int,
1192 << "The subface case given (" << arg1 << ") does not make sense "
1193 << "for the space dimension of the object to which they are applied.");
1194
1195 private:
1200 std::uint8_t value : (dim == 3 ? 4 : 1);
1201 };
1202
1203} // namespace internal
1204
1205
1206
1207template <int dim>
1208struct GeometryInfo;
1209
1210
1211
1231template <>
1232struct GeometryInfo<0>
1233{
1241 static constexpr unsigned int max_children_per_cell = 1;
1242
1246 static constexpr unsigned int faces_per_cell = 0;
1247
1264 static std::array<unsigned int, 0>
1265 face_indices();
1266
1274 static constexpr unsigned int max_children_per_face = 0;
1275
1281 static unsigned int
1282 n_children(const RefinementCase<0> &refinement_case);
1283
1287 static constexpr unsigned int vertices_per_cell = 1;
1288
1307 static std::array<unsigned int, vertices_per_cell>
1308 vertex_indices();
1309
1333 static unsigned int
1334 face_to_cell_vertices(const unsigned int face,
1335 const unsigned int vertex,
1336 const bool face_orientation = true,
1337 const bool face_flip = false,
1338 const bool face_rotation = false);
1339
1354 static unsigned int
1355 face_to_cell_lines(const unsigned int face,
1356 const unsigned int line,
1357 const bool face_orientation = true,
1358 const bool face_flip = false,
1359 const bool face_rotation = false);
1360
1367 static constexpr unsigned int vertices_per_face = 0;
1368
1372 static constexpr unsigned int lines_per_face = 0;
1373
1377 static constexpr unsigned int quads_per_face = 0;
1378
1382 static constexpr unsigned int lines_per_cell = 0;
1383
1387 static constexpr unsigned int quads_per_cell = 0;
1388
1392 static constexpr unsigned int hexes_per_cell = 0;
1393
1411 static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1412
1426 static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1427};
1428
1429
1430
1962template <int dim>
1963struct GeometryInfo
1964{
1972 static constexpr unsigned int max_children_per_cell = 1 << dim;
1973
1977 static constexpr unsigned int faces_per_cell = 2 * dim;
1978
1995 static std_cxx20::ranges::iota_view<unsigned int, unsigned int>
1996 face_indices();
1997
2005 static constexpr unsigned int max_children_per_face =
2006 GeometryInfo<dim - 1>::max_children_per_cell;
2007
2011 static constexpr unsigned int vertices_per_cell = 1 << dim;
2012
2028 static std_cxx20::ranges::iota_view<unsigned int, unsigned int>
2029 vertex_indices();
2030
2034 static constexpr unsigned int vertices_per_face =
2035 GeometryInfo<dim - 1>::vertices_per_cell;
2036
2040 static constexpr unsigned int lines_per_face =
2041 GeometryInfo<dim - 1>::lines_per_cell;
2042
2046 static constexpr unsigned int quads_per_face =
2047 GeometryInfo<dim - 1>::quads_per_cell;
2048
2058 static constexpr unsigned int lines_per_cell =
2059 (2 * GeometryInfo<dim - 1>::lines_per_cell +
2060 GeometryInfo<dim - 1>::vertices_per_cell);
2061
2069 static constexpr unsigned int quads_per_cell =
2070 (2 * GeometryInfo<dim - 1>::quads_per_cell +
2071 GeometryInfo<dim - 1>::lines_per_cell);
2072
2076 static constexpr unsigned int hexes_per_cell =
2077 (2 * GeometryInfo<dim - 1>::hexes_per_cell +
2078 GeometryInfo<dim - 1>::quads_per_cell);
2079
2097 static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
2098 internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal();
2099
2113 static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
2114 internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal();
2115
2126 static constexpr ndarray<unsigned int, vertices_per_cell, dim>
2127 vertex_to_face =
2128 internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face();
2129
2134 static unsigned int
2135 n_children(const RefinementCase<dim> &refinement_case);
2136
2141 static unsigned int
2142 n_subfaces(const internal::SubfaceCase<dim> &subface_case);
2143
2153 static double
2154 subface_ratio(const internal::SubfaceCase<dim> &subface_case,
2155 const unsigned int subface_no);
2156
2162 static RefinementCase<dim - 1>
2163 face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2164 const unsigned int face_no,
2165 const bool face_orientation = true,
2166 const bool face_flip = false,
2167 const bool face_rotation = false);
2168
2174 static RefinementCase<dim>
2175 min_cell_refinement_case_for_face_refinement(
2176 const RefinementCase<dim - 1> &face_refinement_case,
2177 const unsigned int face_no,
2178 const bool face_orientation = true,
2179 const bool face_flip = false,
2180 const bool face_rotation = false);
2181
2186 static RefinementCase<1>
2187 line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2188 const unsigned int line_no);
2189
2194 static RefinementCase<dim>
2195 min_cell_refinement_case_for_line_refinement(const unsigned int line_no);
2196
2243 static unsigned int
2244 child_cell_on_face(const RefinementCase<dim> &ref_case,
2245 const unsigned int face,
2246 const unsigned int subface,
2247 const bool face_orientation = true,
2248 const bool face_flip = false,
2249 const bool face_rotation = false,
2250 const RefinementCase<dim - 1> &face_refinement_case =
2251 RefinementCase<dim - 1>::isotropic_refinement);
2252
2266 static unsigned int
2267 line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
2268
2289 static unsigned int
2290 face_to_cell_vertices(const unsigned int face,
2291 const unsigned int vertex,
2292 const bool face_orientation = true,
2293 const bool face_flip = false,
2294 const bool face_rotation = false);
2295
2307 static unsigned int
2308 face_to_cell_lines(const unsigned int face,
2309 const unsigned int line,
2310 const bool face_orientation = true,
2311 const bool face_flip = false,
2312 const bool face_rotation = false);
2313
2323 static unsigned int
2324 standard_to_real_face_vertex(const unsigned int vertex,
2325 const bool face_orientation = true,
2326 const bool face_flip = false,
2327 const bool face_rotation = false);
2328
2338 static unsigned int
2339 real_to_standard_face_vertex(const unsigned int vertex,
2340 const bool face_orientation = true,
2341 const bool face_flip = false,
2342 const bool face_rotation = false);
2343
2353 static unsigned int
2354 standard_to_real_face_line(const unsigned int line,
2355 const bool face_orientation = true,
2356 const bool face_flip = false,
2357 const bool face_rotation = false);
2358
2364 static unsigned int
2365 standard_to_real_line_vertex(const unsigned int vertex,
2366 const bool line_orientation = true);
2367
2375 static std::array<unsigned int, 2>
2376 standard_quad_vertex_to_line_vertex_index(const unsigned int vertex);
2377
2385 static std::array<unsigned int, 2>
2386 standard_hex_vertex_to_quad_vertex_index(const unsigned int vertex);
2387
2395 static std::array<unsigned int, 2>
2396 standard_hex_line_to_quad_line_index(const unsigned int line);
2397
2407 static unsigned int
2408 real_to_standard_face_line(const unsigned int line,
2409 const bool face_orientation = true,
2410 const bool face_flip = false,
2411 const bool face_rotation = false);
2412
2418 static Point<dim>
2419 unit_cell_vertex(const unsigned int vertex);
2420
2430 static unsigned int
2431 child_cell_from_point(const Point<dim> &p);
2432
2440 static Point<dim>
2441 cell_to_child_coordinates(const Point<dim> &p,
2442 const unsigned int child_index,
2443 const RefinementCase<dim> refine_case =
2444 RefinementCase<dim>::isotropic_refinement);
2445
2451 static Point<dim>
2452 child_to_cell_coordinates(const Point<dim> &p,
2453 const unsigned int child_index,
2454 const RefinementCase<dim> refine_case =
2455 RefinementCase<dim>::isotropic_refinement);
2456
2461 static bool
2462 is_inside_unit_cell(const Point<dim> &p);
2463
2475 static bool
2476 is_inside_unit_cell(const Point<dim> &p, const double eps);
2477
2482 template <typename Number = double>
2483 static Point<dim, Number>
2484 project_to_unit_cell(const Point<dim, Number> &p);
2485
2491 static double
2492 distance_to_unit_cell(const Point<dim> &p);
2493
2498 static double
2499 d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2500
2505 static Tensor<1, dim>
2506 d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2507
2559 template <int spacedim>
2560 static void
2561 alternating_form_at_vertices
2562#ifndef DEAL_II_CXX14_CONSTEXPR_BUG
2563 (const Point<spacedim> (&vertices)[vertices_per_cell],
2564 Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2565#else
2566 (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2567#endif
2568
2578 static constexpr std::array<unsigned int, faces_per_cell>
2579 unit_normal_direction =
2580 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction();
2581
2598 static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2599 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation();
2600
2611 static constexpr std::array<Tensor<1, dim>, faces_per_cell>
2612 unit_normal_vector =
2613 internal::GeometryInfoHelper::Initializers<dim>::unit_normal_vector();
2614
2628 static constexpr ndarray<Tensor<1, dim>, faces_per_cell, dim - 1>
2629 unit_tangential_vectors = internal::GeometryInfoHelper::Initializers<
2630 dim>::unit_tangential_vectors();
2631
2637 static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2638 internal::GeometryInfoHelper::Initializers<dim>::opposite_face();
2639
2640
2644 DeclException1(ExcInvalidCoordinate,
2645 double,
2646 << "The coordinates must satisfy 0 <= x_i <= 1, "
2647 << "but here we have x_i=" << arg1);
2648
2652 DeclException3(ExcInvalidSubface,
2653 int,
2654 int,
2655 int,
2656 << "RefinementCase<dim> " << arg1 << ": face " << arg2
2657 << " has no subface " << arg3);
2658};
2659
2660
2661
2662#ifndef DOXYGEN
2663
2664
2665/* -------------- declaration of explicit specializations ------------- */
2666
2667
2668template <>
2669Tensor<1, 1>
2670GeometryInfo<1>::d_linear_shape_function_gradient(const Point<1> &xi,
2671 const unsigned int i);
2672template <>
2673Tensor<1, 2>
2674GeometryInfo<2>::d_linear_shape_function_gradient(const Point<2> &xi,
2675 const unsigned int i);
2676template <>
2677Tensor<1, 3>
2678GeometryInfo<3>::d_linear_shape_function_gradient(const Point<3> &xi,
2679 const unsigned int i);
2680
2681
2682
2683/* -------------- inline functions ------------- */
2684
2685
2686inline GeometryPrimitive::GeometryPrimitive(const Object object)
2687 : object(object)
2688{}
2689
2690
2691
2692inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2693 : object(static_cast<Object>(object_dimension))
2694{}
2695
2696
2697inline GeometryPrimitive::operator unsigned int() const
2698{
2699 return static_cast<unsigned int>(object);
2700}
2701
2702
2703
2704namespace internal
2705{
2706 template <int dim>
2707 inline SubfaceCase<dim>::SubfaceCase(
2708 const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2709 : value(subface_possibility)
2710 {}
2711
2712
2713 template <int dim>
2714 inline SubfaceCase<dim>::operator std::uint8_t() const
2715 {
2716 return value;
2717 }
2718
2719
2720} // namespace internal
2721
2722
2723template <int dim>
2724inline RefinementCase<dim>
2725RefinementCase<dim>::cut_axis(const unsigned int)
2726{
2727 DEAL_II_ASSERT_UNREACHABLE();
2728 return static_cast<std::uint8_t>(-1);
2729}
2730
2731
2732template <>
2733inline RefinementCase<1>
2734RefinementCase<1>::cut_axis(const unsigned int i)
2735{
2736 AssertIndexRange(i, 1);
2737
2738 const RefinementCase options[1] = {RefinementPossibilities<1>::cut_x};
2739 return options[i];
2740}
2741
2742
2743
2744template <>
2745inline RefinementCase<2>
2746RefinementCase<2>::cut_axis(const unsigned int i)
2747{
2748 AssertIndexRange(i, 2);
2749
2750 const RefinementCase options[2] = {RefinementPossibilities<2>::cut_x,
2751 RefinementPossibilities<2>::cut_y};
2752 return options[i];
2753}
2754
2755
2756
2757template <>
2758inline RefinementCase<3>
2759RefinementCase<3>::cut_axis(const unsigned int i)
2760{
2761 AssertIndexRange(i, 3);
2762
2763 const RefinementCase options[3] = {RefinementPossibilities<3>::cut_x,
2764 RefinementPossibilities<3>::cut_y,
2765 RefinementPossibilities<3>::cut_z};
2766 return options[i];
2767}
2768
2769
2770
2771template <>
2772inline std::array<RefinementCase<1>, 2>
2773RefinementCase<1>::all_refinement_cases()
2774{
2775 return {{RefinementPossibilities<1>::no_refinement,
2776 RefinementPossibilities<1>::cut_x}};
2777}
2778
2779
2780
2781template <>
2782inline std::array<RefinementCase<2>, 4>
2783RefinementCase<2>::all_refinement_cases()
2784{
2785 return {{RefinementPossibilities<2>::no_refinement,
2786 RefinementPossibilities<2>::cut_x,
2787 RefinementPossibilities<2>::cut_y,
2788 RefinementPossibilities<2>::cut_xy}};
2789}
2790
2791
2792
2793template <>
2794inline std::array<RefinementCase<3>, 8>
2795RefinementCase<3>::all_refinement_cases()
2796{
2797 return {{RefinementPossibilities<3>::no_refinement,
2798 RefinementPossibilities<3>::cut_x,
2799 RefinementPossibilities<3>::cut_y,
2800 RefinementPossibilities<3>::cut_xy,
2801 RefinementPossibilities<3>::cut_z,
2802 RefinementPossibilities<3>::cut_xz,
2803 RefinementPossibilities<3>::cut_yz,
2804 RefinementPossibilities<3>::cut_xyz}};
2805}
2806
2807
2808
2809template <int dim>
2810inline RefinementCase<dim>::RefinementCase()
2811 : value(RefinementPossibilities<dim>::no_refinement)
2812{}
2813
2814
2815
2816template <int dim>
2817inline RefinementCase<dim>::RefinementCase(
2818 const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2819 : value(refinement_case)
2820{
2821 // check that only those bits of
2822 // the given argument are set that
2823 // make sense for a given space
2824 // dimension
2825 Assert((refinement_case &
2826 RefinementPossibilities<dim>::isotropic_refinement) ==
2827 refinement_case,
2828 ExcInvalidRefinementCase(refinement_case));
2829}
2830
2831
2832
2833template <int dim>
2834inline RefinementCase<dim>::RefinementCase(const std::uint8_t refinement_case)
2835 : value(refinement_case)
2836{
2837 // check that only those bits of
2838 // the given argument are set that
2839 // make sense for a given space
2840 // dimension
2841 Assert((refinement_case &
2842 RefinementPossibilities<dim>::isotropic_refinement) ==
2843 refinement_case,
2844 ExcInvalidRefinementCase(refinement_case));
2845}
2846
2847
2848
2849template <int dim>
2850inline DEAL_II_HOST_DEVICE RefinementCase<dim>::operator std::uint8_t() const
2851{
2852 return value;
2853}
2854
2855
2856
2857template <int dim>
2858inline RefinementCase<dim>
2859RefinementCase<dim>::operator|(const RefinementCase<dim> &r) const
2860{
2861 return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2862}
2863
2864
2865
2866template <int dim>
2867inline RefinementCase<dim>
2868RefinementCase<dim>::operator&(const RefinementCase<dim> &r) const
2869{
2870 return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2871}
2872
2873
2874
2875template <int dim>
2876inline RefinementCase<dim>
2877RefinementCase<dim>::operator~() const
2878{
2879 return RefinementCase<dim>(static_cast<std::uint8_t>(
2880 (~value) & RefinementPossibilities<dim>::isotropic_refinement));
2881}
2882
2883
2884
2885template <int dim>
2886inline std::size_t
2887RefinementCase<dim>::memory_consumption()
2888{
2889 return sizeof(RefinementCase<dim>);
2890}
2891
2892
2893
2894template <int dim>
2895template <class Archive>
2896inline void
2897RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2898{
2899 // serialization can't deal with bitfields, so copy from/to a full sized
2900 // std::uint8_t
2901 std::uint8_t uchar_value = value;
2902 ar &uchar_value;
2903 value = uchar_value;
2904}
2905
2906
2907
2908template <>
2909inline Point<1>
2910GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2911{
2912 AssertIndexRange(vertex, vertices_per_cell);
2913
2914 return Point<1>(static_cast<double>(vertex));
2915}
2916
2917
2918
2919template <>
2920inline Point<2>
2921GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2922{
2923 AssertIndexRange(vertex, vertices_per_cell);
2924
2925 return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2926}
2927
2928
2929
2930template <>
2931inline Point<3>
2932GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2933{
2934 AssertIndexRange(vertex, vertices_per_cell);
2935
2936 return {static_cast<double>(vertex % 2),
2937 static_cast<double>(vertex / 2 % 2),
2938 static_cast<double>(vertex / 4)};
2939}
2940
2941
2942
2943inline std::array<unsigned int, 0>
2944GeometryInfo<0>::face_indices()
2945{
2946 return {{}};
2947}
2948
2949
2950
2951inline std::array<unsigned int, 1>
2952GeometryInfo<0>::vertex_indices()
2953{
2954 return {{0}};
2955}
2956
2957
2958
2959template <int dim>
2960inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
2961GeometryInfo<dim>::face_indices()
2962{
2963 return {0U, faces_per_cell};
2964}
2965
2966
2967
2968template <int dim>
2969inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
2970GeometryInfo<dim>::vertex_indices()
2971{
2972 return {0U, vertices_per_cell};
2973}
2974
2975
2976
2977template <int dim>
2978inline Point<dim>
2979GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2980{
2981 DEAL_II_NOT_IMPLEMENTED();
2982
2983 return {};
2984}
2985
2986
2987
2988template <>
2989inline unsigned int
2990GeometryInfo<1>::child_cell_from_point(const Point<1> &p)
2991{
2992 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2993
2994 return (p[0] <= 0.5 ? 0 : 1);
2995}
2996
2997
2998
2999template <>
3000inline unsigned int
3001GeometryInfo<2>::child_cell_from_point(const Point<2> &p)
3002{
3003 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
3004 Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
3005
3006 return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
3007}
3008
3009
3010
3011template <>
3012inline unsigned int
3013GeometryInfo<3>::child_cell_from_point(const Point<3> &p)
3014{
3015 Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
3016 Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
3017 Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
3018
3019 return (p[0] <= 0.5 ?
3020 (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
3021 (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
3022}
3023
3024
3025template <int dim>
3026inline unsigned int
3027GeometryInfo<dim>::child_cell_from_point(const Point<dim> &)
3028{
3029 DEAL_II_NOT_IMPLEMENTED();
3030
3031 return 0;
3032}
3033
3034
3035
3036template <>
3037inline Point<1>
3038GeometryInfo<1>::cell_to_child_coordinates(const Point<1> &p,
3039 const unsigned int child_index,
3040 const RefinementCase<1> refine_case)
3041
3042{
3043 AssertIndexRange(child_index, 2);
3044 Assert(refine_case == RefinementCase<1>::cut_x, ExcInternalError());
3045 (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3046
3047 return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
3048}
3049
3050
3051
3052template <>
3053inline Point<2>
3054GeometryInfo<2>::cell_to_child_coordinates(const Point<2> &p,
3055 const unsigned int child_index,
3056 const RefinementCase<2> refine_case)
3057
3058{
3059 AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3060
3061 Point<2> point = p;
3062 switch (refine_case)
3063 {
3064 case RefinementCase<2>::cut_x:
3065 point[0] *= 2.0;
3066 if (child_index == 1)
3067 point[0] -= 1.0;
3068 break;
3069 case RefinementCase<2>::cut_y:
3070 point[1] *= 2.0;
3071 if (child_index == 1)
3072 point[1] -= 1.0;
3073 break;
3074 case RefinementCase<2>::cut_xy:
3075 point *= 2.0;
3076 point -= unit_cell_vertex(child_index);
3077 break;
3078 default:
3079 DEAL_II_ASSERT_UNREACHABLE();
3080 }
3081
3082 return point;
3083}
3084
3085
3086
3087template <>
3088inline Point<3>
3089GeometryInfo<3>::cell_to_child_coordinates(const Point<3> &p,
3090 const unsigned int child_index,
3091 const RefinementCase<3> refine_case)
3092
3093{
3094 AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3095
3096 Point<3> point = p;
3097 // there might be a cleverer way to do
3098 // this, but since this function is called
3099 // in very few places for initialization
3100 // purposes only, I don't care at the
3101 // moment
3102 switch (refine_case)
3103 {
3104 case RefinementCase<3>::cut_x:
3105 point[0] *= 2.0;
3106 if (child_index == 1)
3107 point[0] -= 1.0;
3108 break;
3109 case RefinementCase<3>::cut_y:
3110 point[1] *= 2.0;
3111 if (child_index == 1)
3112 point[1] -= 1.0;
3113 break;
3114 case RefinementCase<3>::cut_z:
3115 point[2] *= 2.0;
3116 if (child_index == 1)
3117 point[2] -= 1.0;
3118 break;
3119 case RefinementCase<3>::cut_xy:
3120 point[0] *= 2.0;
3121 point[1] *= 2.0;
3122 if (child_index % 2 == 1)
3123 point[0] -= 1.0;
3124 if (child_index / 2 == 1)
3125 point[1] -= 1.0;
3126 break;
3127 case RefinementCase<3>::cut_xz:
3128 // careful, this is slightly
3129 // different from xy and yz due to
3130 // different internal numbering of
3131 // children!
3132 point[0] *= 2.0;
3133 point[2] *= 2.0;
3134 if (child_index / 2 == 1)
3135 point[0] -= 1.0;
3136 if (child_index % 2 == 1)
3137 point[2] -= 1.0;
3138 break;
3139 case RefinementCase<3>::cut_yz:
3140 point[1] *= 2.0;
3141 point[2] *= 2.0;
3142 if (child_index % 2 == 1)
3143 point[1] -= 1.0;
3144 if (child_index / 2 == 1)
3145 point[2] -= 1.0;
3146 break;
3147 case RefinementCase<3>::cut_xyz:
3148 point *= 2.0;
3149 point -= unit_cell_vertex(child_index);
3150 break;
3151 default:
3152 DEAL_II_ASSERT_UNREACHABLE();
3153 }
3154
3155 return point;
3156}
3157
3158
3159
3160template <int dim>
3161inline Point<dim>
3162GeometryInfo<dim>::cell_to_child_coordinates(
3163 const Point<dim> & /*p*/,
3164 const unsigned int /*child_index*/,
3165 const RefinementCase<dim> /*refine_case*/)
3166
3167{
3168 DEAL_II_NOT_IMPLEMENTED();
3169 return {};
3170}
3171
3172
3173
3174template <>
3175inline Point<1>
3176GeometryInfo<1>::child_to_cell_coordinates(const Point<1> &p,
3177 const unsigned int child_index,
3178 const RefinementCase<1> refine_case)
3179
3180{
3181 AssertIndexRange(child_index, 2);
3182 Assert(refine_case == RefinementCase<1>::cut_x, ExcInternalError());
3183 (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3184
3185 return (p + unit_cell_vertex(child_index)) * 0.5;
3186}
3187
3188
3189
3190template <>
3191inline Point<3>
3192GeometryInfo<3>::child_to_cell_coordinates(const Point<3> &p,
3193 const unsigned int child_index,
3194 const RefinementCase<3> refine_case)
3195
3196{
3197 AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3198
3199 Point<3> point = p;
3200 // there might be a cleverer way to do
3201 // this, but since this function is called
3202 // in very few places for initialization
3203 // purposes only, I don't care at the
3204 // moment
3205 switch (refine_case)
3206 {
3207 case RefinementCase<3>::cut_x:
3208 if (child_index == 1)
3209 point[0] += 1.0;
3210 point[0] *= 0.5;
3211 break;
3212 case RefinementCase<3>::cut_y:
3213 if (child_index == 1)
3214 point[1] += 1.0;
3215 point[1] *= 0.5;
3216 break;
3217 case RefinementCase<3>::cut_z:
3218 if (child_index == 1)
3219 point[2] += 1.0;
3220 point[2] *= 0.5;
3221 break;
3222 case RefinementCase<3>::cut_xy:
3223 if (child_index % 2 == 1)
3224 point[0] += 1.0;
3225 if (child_index / 2 == 1)
3226 point[1] += 1.0;
3227 point[0] *= 0.5;
3228 point[1] *= 0.5;
3229 break;
3230 case RefinementCase<3>::cut_xz:
3231 // careful, this is slightly
3232 // different from xy and yz due to
3233 // different internal numbering of
3234 // children!
3235 if (child_index / 2 == 1)
3236 point[0] += 1.0;
3237 if (child_index % 2 == 1)
3238 point[2] += 1.0;
3239 point[0] *= 0.5;
3240 point[2] *= 0.5;
3241 break;
3242 case RefinementCase<3>::cut_yz:
3243 if (child_index % 2 == 1)
3244 point[1] += 1.0;
3245 if (child_index / 2 == 1)
3246 point[2] += 1.0;
3247 point[1] *= 0.5;
3248 point[2] *= 0.5;
3249 break;
3250 case RefinementCase<3>::cut_xyz:
3251 point += unit_cell_vertex(child_index);
3252 point *= 0.5;
3253 break;
3254 default:
3255 DEAL_II_ASSERT_UNREACHABLE();
3256 }
3257
3258 return point;
3259}
3260
3261
3262
3263template <>
3264inline Point<2>
3265GeometryInfo<2>::child_to_cell_coordinates(const Point<2> &p,
3266 const unsigned int child_index,
3267 const RefinementCase<2> refine_case)
3268{
3269 AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3270
3271 Point<2> point = p;
3272 switch (refine_case)
3273 {
3274 case RefinementCase<2>::cut_x:
3275 if (child_index == 1)
3276 point[0] += 1.0;
3277 point[0] *= 0.5;
3278 break;
3279 case RefinementCase<2>::cut_y:
3280 if (child_index == 1)
3281 point[1] += 1.0;
3282 point[1] *= 0.5;
3283 break;
3284 case RefinementCase<2>::cut_xy:
3285 point += unit_cell_vertex(child_index);
3286 point *= 0.5;
3287 break;
3288 default:
3289 DEAL_II_ASSERT_UNREACHABLE();
3290 }
3291
3292 return point;
3293}
3294
3295
3296
3297template <int dim>
3298inline Point<dim>
3299GeometryInfo<dim>::child_to_cell_coordinates(
3300 const Point<dim> & /*p*/,
3301 const unsigned int /*child_index*/,
3302 const RefinementCase<dim> /*refine_case*/)
3303{
3304 DEAL_II_NOT_IMPLEMENTED();
3305 return {};
3306}
3307
3308
3309
3310template <int dim>
3311inline bool
3312GeometryInfo<dim>::is_inside_unit_cell(const Point<dim> &)
3313{
3314 DEAL_II_NOT_IMPLEMENTED();
3315 return false;
3316}
3317
3318template <>
3319inline bool
3320GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p)
3321{
3322 return (p[0] >= 0.) && (p[0] <= 1.);
3323}
3324
3325
3326
3327template <>
3328inline bool
3329GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p)
3330{
3331 return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
3332}
3333
3334
3335
3336template <>
3337inline bool
3338GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p)
3339{
3340 return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
3341 (p[2] >= 0.) && (p[2] <= 1.);
3342}
3343
3344
3345
3346template <int dim>
3347inline bool
3348GeometryInfo<dim>::is_inside_unit_cell(const Point<dim> &, const double)
3349{
3350 DEAL_II_NOT_IMPLEMENTED();
3351 return false;
3352}
3353
3354template <>
3355inline bool
3356GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
3357{
3358 return (p[0] >= -eps) && (p[0] <= 1. + eps);
3359}
3360
3361
3362
3363template <>
3364inline bool
3365GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
3366{
3367 const double l = -eps, u = 1 + eps;
3368 return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
3369}
3370
3371
3372
3373template <>
3374inline bool
3375GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
3376{
3377 const double l = -eps, u = 1.0 + eps;
3378 return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
3379 (p[2] >= l) && (p[2] <= u);
3380}
3381
3382
3383
3384template <>
3385inline unsigned int
3386GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
3387 const unsigned int vertex)
3388{
3389 (void)line;
3390 AssertIndexRange(line, lines_per_cell);
3391 AssertIndexRange(vertex, 2);
3392
3393 return vertex;
3394}
3395
3396
3397template <>
3398inline unsigned int
3399GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
3400 const unsigned int vertex)
3401{
3402 constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
3403 return cell_vertices[line][vertex];
3404}
3405
3406
3407
3408template <>
3409inline unsigned int
3410GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
3411 const unsigned int vertex)
3412{
3413 AssertIndexRange(line, lines_per_cell);
3414 AssertIndexRange(vertex, 2);
3415
3416 constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3417 {1, 3},
3418 {0, 1},
3419 {2, 3},
3420 {4, 6}, // top face
3421 {5, 7},
3422 {4, 5},
3423 {6, 7},
3424 {0,
3425 4}, // connects of bottom
3426 {1, 5}, // top face
3427 {2, 6},
3428 {3, 7}};
3429 return vertices[line][vertex];
3430}
3431
3432
3433template <>
3434inline unsigned int
3435GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3436{
3437 DEAL_II_NOT_IMPLEMENTED();
3438 return numbers::invalid_unsigned_int;
3439}
3440
3441template <>
3442inline unsigned int
3443GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3444 const bool face_orientation,
3445 const bool face_flip,
3446 const bool face_rotation)
3447{
3448 AssertIndexRange(vertex, GeometryInfo<3>::vertices_per_face);
3449
3450 // set up a table to make sure that
3451 // we handle non-standard faces correctly
3452 //
3453 // so set up a table that for each vertex (of
3454 // a quad in standard position) describes
3455 // which vertex to take
3456 //
3457 // first index: four vertices 0...3
3458 //
3459 // second index: face_orientation; 0:
3460 // opposite normal, 1: standard
3461 //
3462 // third index: face_flip; 0: standard, 1:
3463 // face rotated by 180 degrees
3464 //
3465 // forth index: face_rotation: 0: standard,
3466 // 1: face rotated by 90 degrees
3467
3468 constexpr unsigned int vertex_translation[4][2][2][2] = {
3469 {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3470 // face_rotation=false and true
3471 {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3472 // face_rotation=false and true
3473 {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3474 // face_rotation=false and true
3475 {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3476 // face_rotation=false and true
3477
3478 {{{2, 3}, // vertex 1 ...
3479 {1, 0}},
3480 {{1, 0}, {2, 3}}},
3481
3482 {{{1, 0}, // vertex 2 ...
3483 {2, 3}},
3484 {{2, 3}, {1, 0}}},
3485
3486 {{{3, 1}, // vertex 3 ...
3487 {0, 2}},
3488 {{3, 1}, {0, 2}}}};
3489
3490 return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3491}
3492
3493
3494
3495template <int dim>
3496inline unsigned int
3497GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3498 const bool,
3499 const bool,
3500 const bool)
3501{
3502 Assert(dim > 1, ExcImpossibleInDim(dim));
3503 AssertIndexRange(vertex, GeometryInfo<dim>::vertices_per_face);
3504 return vertex;
3505}
3506
3507template <int dim>
3508inline unsigned int
3509GeometryInfo<dim>::n_children(const RefinementCase<dim> &ref_case)
3510{
3511 constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3512 0, 2, 2, 4, 2, 4, 4, 8};
3513
3514 return n_children[ref_case];
3515}
3516
3517
3518
3519template <int dim>
3520inline unsigned int
3521GeometryInfo<dim>::n_subfaces(const internal::SubfaceCase<dim> &)
3522{
3523 DEAL_II_NOT_IMPLEMENTED();
3524 return 0;
3525}
3526
3527template <>
3528inline unsigned int
3529GeometryInfo<1>::n_subfaces(const internal::SubfaceCase<1> &)
3530{
3531 Assert(false, ExcImpossibleInDim(1));
3532 return 0;
3533}
3534
3535template <>
3536inline unsigned int
3537GeometryInfo<2>::n_subfaces(const internal::SubfaceCase<2> &subface_case)
3538{
3539 return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3540}
3541
3542
3543
3544template <>
3545inline unsigned int
3546GeometryInfo<3>::n_subfaces(const internal::SubfaceCase<3> &subface_case)
3547{
3548 const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3549 0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3550 return nsubs[subface_case];
3551}
3552
3553
3554
3555template <int dim>
3556inline double
3557GeometryInfo<dim>::subface_ratio(const internal::SubfaceCase<dim> &,
3558 const unsigned int)
3559{
3560 DEAL_II_NOT_IMPLEMENTED();
3561 return 0.;
3562}
3563
3564template <>
3565inline double
3566GeometryInfo<1>::subface_ratio(const internal::SubfaceCase<1> &,
3567 const unsigned int)
3568{
3569 return 1;
3570}
3571
3572
3573template <>
3574inline double
3575GeometryInfo<2>::subface_ratio(const internal::SubfaceCase<2> &subface_case,
3576 const unsigned int)
3577{
3578 double ratio = 1;
3579 switch (subface_case)
3580 {
3581 case internal::SubfaceCase<2>::case_none:
3582 // Here, an
3583 // Assert(false,ExcInternalError())
3584 // would be the right
3585 // choice, but
3586 // unfortunately the
3587 // current function is
3588 // also called for faces
3589 // without children (see
3590 // tests/fe/mapping.cc).
3591 // Assert(false, ExcMessage("Face has no subfaces."));
3592 // Furthermore, assign
3593 // following value as
3594 // otherwise the
3595 // bits/volume_x tests
3596 // break
3597 ratio = 1. / GeometryInfo<2>::max_children_per_face;
3598 break;
3599 case internal::SubfaceCase<2>::case_x:
3600 ratio = 0.5;
3601 break;
3602 default:
3603 // there should be no
3604 // cases left
3605 DEAL_II_ASSERT_UNREACHABLE();
3606 break;
3607 }
3608
3609 return ratio;
3610}
3611
3612
3613template <>
3614inline double
3615GeometryInfo<3>::subface_ratio(const internal::SubfaceCase<3> &subface_case,
3616 const unsigned int subface_no)
3617{
3618 double ratio = 1;
3619 switch (subface_case)
3620 {
3621 case internal::SubfaceCase<3>::case_none:
3622 // Here, an
3623 // Assert(false,ExcInternalError())
3624 // would be the right
3625 // choice, but
3626 // unfortunately the
3627 // current function is
3628 // also called for faces
3629 // without children (see
3630 // tests/bits/mesh_3d_16.cc). Add
3631 // following switch to
3632 // avoid diffs in
3633 // tests/bits/mesh_3d_16
3634 ratio = 1. / GeometryInfo<3>::max_children_per_face;
3635 break;
3636 case internal::SubfaceCase<3>::case_x:
3637 case internal::SubfaceCase<3>::case_y:
3638 ratio = 0.5;
3639 break;
3640 case internal::SubfaceCase<3>::case_xy:
3641 case internal::SubfaceCase<3>::case_x1y2y:
3642 case internal::SubfaceCase<3>::case_y1x2x:
3643 ratio = 0.25;
3644 break;
3645 case internal::SubfaceCase<3>::case_x1y:
3646 case internal::SubfaceCase<3>::case_y1x:
3647 if (subface_no < 2)
3648 ratio = 0.25;
3649 else
3650 ratio = 0.5;
3651 break;
3652 case internal::SubfaceCase<3>::case_x2y:
3653 case internal::SubfaceCase<3>::case_y2x:
3654 if (subface_no == 0)
3655 ratio = 0.5;
3656 else
3657 ratio = 0.25;
3658 break;
3659 default:
3660 // there should be no
3661 // cases left
3662 DEAL_II_ASSERT_UNREACHABLE();
3663 break;
3664 }
3665
3666 return ratio;
3667}
3668
3669
3670
3671template <int dim>
3672RefinementCase<dim - 1> inline GeometryInfo<dim>::face_refinement_case(
3673 const RefinementCase<dim> &,
3674 const unsigned int,
3675 const bool,
3676 const bool,
3677 const bool)
3678{
3679 DEAL_II_NOT_IMPLEMENTED();
3680 return RefinementCase<dim - 1>::no_refinement;
3681}
3682
3683template <>
3684RefinementCase<0> inline GeometryInfo<1>::face_refinement_case(
3685 const RefinementCase<1> &,
3686 const unsigned int,
3687 const bool,
3688 const bool,
3689 const bool)
3690{
3691 Assert(false, ExcImpossibleInDim(1));
3692
3693 return RefinementCase<0>::no_refinement;
3694}
3695
3696
3697template <>
3698inline RefinementCase<1>
3699GeometryInfo<2>::face_refinement_case(
3700 const RefinementCase<2> &cell_refinement_case,
3701 const unsigned int face_no,
3702 const bool,
3703 const bool,
3704 const bool)
3705{
3706 const unsigned int dim = 2;
3707 AssertIndexRange(cell_refinement_case,
3708 RefinementCase<dim>::isotropic_refinement + 1);
3709 AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
3710
3711 // simple special case
3712 if (cell_refinement_case == RefinementCase<dim>::cut_xy)
3713 return RefinementCase<1>::cut_x;
3714
3715 const RefinementCase<dim - 1>
3716 ref_cases[RefinementCase<dim>::isotropic_refinement +
3717 1][GeometryInfo<dim>::faces_per_cell / 2] = {
3718 {RefinementCase<dim - 1>::no_refinement, // no_refinement
3719 RefinementCase<dim - 1>::no_refinement},
3720
3721 {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3722
3723 {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3724
3725 {RefinementCase<dim - 1>::cut_x, // cut_xy
3726 RefinementCase<dim - 1>::cut_x}};
3727
3728 return ref_cases[cell_refinement_case][face_no / 2];
3729}
3730
3731
3732template <>
3733inline RefinementCase<2>
3734GeometryInfo<3>::face_refinement_case(
3735 const RefinementCase<3> &cell_refinement_case,
3736 const unsigned int face_no,
3737 const bool face_orientation,
3738 const bool /*face_flip*/,
3739 const bool face_rotation)
3740{
3741 const unsigned int dim = 3;
3742 AssertIndexRange(cell_refinement_case,
3743 RefinementCase<dim>::isotropic_refinement + 1);
3744 AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
3745
3746 // simple special case
3747 if (cell_refinement_case == RefinementCase<dim>::cut_xyz)
3748 return RefinementCase<dim - 1>::cut_xy;
3749
3750 const RefinementCase<dim - 1>
3751 ref_cases[RefinementCase<dim>::isotropic_refinement + 1]
3752 [GeometryInfo<dim>::faces_per_cell / 2] = {
3753 {RefinementCase<dim - 1>::no_refinement, // no_refinement
3754 RefinementCase<dim - 1>::no_refinement,
3755 RefinementCase<dim - 1>::no_refinement},
3756
3757 {RefinementCase<dim - 1>::no_refinement, // cut_x
3758 RefinementCase<dim - 1>::cut_y,
3759 RefinementCase<dim - 1>::cut_x},
3760
3761 {RefinementCase<dim - 1>::cut_x, // cut_y
3762 RefinementCase<dim - 1>::no_refinement,
3763 RefinementCase<dim - 1>::cut_y},
3764
3765 {RefinementCase<dim - 1>::cut_x, // cut_xy
3766 RefinementCase<dim - 1>::cut_y,
3767 RefinementCase<dim - 1>::cut_xy},
3768
3769 {RefinementCase<dim - 1>::cut_y, // cut_z
3770 RefinementCase<dim - 1>::cut_x,
3771 RefinementCase<dim - 1>::no_refinement},
3772
3773 {RefinementCase<dim - 1>::cut_y, // cut_xz
3774 RefinementCase<dim - 1>::cut_xy,
3775 RefinementCase<dim - 1>::cut_x},
3776
3777 {RefinementCase<dim - 1>::cut_xy, // cut_yz
3778 RefinementCase<dim - 1>::cut_x,
3779 RefinementCase<dim - 1>::cut_y},
3780
3781 {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3782 RefinementCase<dim - 1>::cut_xy,
3783 RefinementCase<dim - 1>::cut_xy},
3784 };
3785
3786 const RefinementCase<dim - 1> ref_case =
3787 ref_cases[cell_refinement_case][face_no / 2];
3788
3789 const RefinementCase<dim - 1> flip[4] = {
3790 RefinementCase<dim - 1>::no_refinement,
3791 RefinementCase<dim - 1>::cut_y,
3792 RefinementCase<dim - 1>::cut_x,
3793 RefinementCase<dim - 1>::cut_xy};
3794
3795 // correct the ref_case for face_orientation
3796 // and face_rotation. for face_orientation,
3797 // 'true' is the default value whereas for
3798 // face_rotation, 'false' is standard. If
3799 // <tt>face_rotation==face_orientation</tt>,
3800 // then one of them is non-standard and we
3801 // have to swap cut_x and cut_y, otherwise no
3802 // change is necessary. face_flip has no
3803 // influence. however, in order to keep the
3804 // interface consistent with other functions,
3805 // we still include it as an argument to this
3806 // function
3807 return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3808}
3809
3810
3811
3812template <int dim>
3813inline RefinementCase<1>
3814GeometryInfo<dim>::line_refinement_case(const RefinementCase<dim> &,
3815 const unsigned int)
3816{
3817 DEAL_II_NOT_IMPLEMENTED();
3818 return RefinementCase<1>::no_refinement;
3819}
3820
3821template <>
3822inline RefinementCase<1>
3823GeometryInfo<1>::line_refinement_case(
3824 const RefinementCase<1> &cell_refinement_case,
3825 const unsigned int line_no)
3826{
3827 (void)line_no;
3828 const unsigned int dim = 1;
3829 (void)dim;
3830 AssertIndexRange(cell_refinement_case,
3831 RefinementCase<dim>::isotropic_refinement + 1);
3832 AssertIndexRange(line_no, GeometryInfo<dim>::lines_per_cell);
3833
3834 return cell_refinement_case;
3835}
3836
3837
3838template <>
3839inline RefinementCase<1>
3840GeometryInfo<2>::line_refinement_case(
3841 const RefinementCase<2> &cell_refinement_case,
3842 const unsigned int line_no)
3843{
3844 // Assertions are in face_refinement_case()
3845 return face_refinement_case(cell_refinement_case, line_no);
3846}
3847
3848
3849template <>
3850inline RefinementCase<1>
3851GeometryInfo<3>::line_refinement_case(
3852 const RefinementCase<3> &cell_refinement_case,
3853 const unsigned int line_no)
3854{
3855 const unsigned int dim = 3;
3856 AssertIndexRange(cell_refinement_case,
3857 RefinementCase<dim>::isotropic_refinement + 1);
3858 AssertIndexRange(line_no, GeometryInfo<dim>::lines_per_cell);
3859
3860 // simple special case
3861 if (cell_refinement_case == RefinementCase<dim>::cut_xyz)
3862 return RefinementCase<1>::cut_x;
3863
3864 // array indicating, which simple refine
3865 // case cuts a line in direction x, y or
3866 // z. For example, cut_y and everything
3867 // containing cut_y (cut_xy, cut_yz,
3868 // cut_xyz) cuts lines, which are in y
3869 // direction.
3870 const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3871 RefinementCase<dim>::cut_y,
3872 RefinementCase<dim>::cut_z};
3873
3874 // order the direction of lines
3875 // 0->x, 1->y, 2->z
3876 const unsigned int direction[lines_per_cell] = {
3877 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3878
3879 return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3880 RefinementCase<1>::cut_x :
3881 RefinementCase<1>::no_refinement);
3882}
3883
3884
3885
3886template <int dim>
3887inline RefinementCase<dim>
3888GeometryInfo<dim>::min_cell_refinement_case_for_face_refinement(
3889 const RefinementCase<dim - 1> &,
3890 const unsigned int,
3891 const bool,
3892 const bool,
3893 const bool)
3894{
3895 DEAL_II_NOT_IMPLEMENTED();
3896
3897 return RefinementCase<dim>::no_refinement;
3898}
3899
3900template <>
3901inline RefinementCase<1>
3902GeometryInfo<1>::min_cell_refinement_case_for_face_refinement(
3903 const RefinementCase<0> &,
3904 const unsigned int,
3905 const bool,
3906 const bool,
3907 const bool)
3908{
3909 const unsigned int dim = 1;
3910 Assert(false, ExcImpossibleInDim(dim));
3911
3912 return RefinementCase<dim>::no_refinement;
3913}
3914
3915
3916template <>
3917inline RefinementCase<2>
3918GeometryInfo<2>::min_cell_refinement_case_for_face_refinement(
3919 const RefinementCase<1> &face_refinement_case,
3920 const unsigned int face_no,
3921 const bool,
3922 const bool,
3923 const bool)
3924{
3925 const unsigned int dim = 2;
3926 AssertIndexRange(face_refinement_case,
3927 RefinementCase<dim - 1>::isotropic_refinement + 1);
3928 AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
3929
3930 if (face_refinement_case == RefinementCase<dim>::cut_x)
3931 return (face_no / 2) != 0u ? RefinementCase<dim>::cut_x :
3932 RefinementCase<dim>::cut_y;
3933 else
3934 return RefinementCase<dim>::no_refinement;
3935}
3936
3937
3938template <>
3939inline RefinementCase<3>
3940GeometryInfo<3>::min_cell_refinement_case_for_face_refinement(
3941 const RefinementCase<2> &face_refinement_case,
3942 const unsigned int face_no,
3943 const bool face_orientation,
3944 const bool /*face_flip*/,
3945 const bool face_rotation)
3946{
3947 const unsigned int dim = 3;
3948 AssertIndexRange(face_refinement_case,
3949 RefinementCase<dim - 1>::isotropic_refinement + 1);
3950 AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
3951
3952 const RefinementCase<2> flip[4] = {RefinementCase<2>::no_refinement,
3953 RefinementCase<2>::cut_y,
3954 RefinementCase<2>::cut_x,
3955 RefinementCase<2>::cut_xy};
3956
3957 // correct the face_refinement_case for
3958 // face_orientation and face_rotation. for
3959 // face_orientation, 'true' is the default
3960 // value whereas for face_rotation, 'false'
3961 // is standard. If
3962 // <tt>face_rotation==face_orientation</tt>,
3963 // then one of them is non-standard and we
3964 // have to swap cut_x and cut_y, otherwise no
3965 // change is necessary. face_flip has no
3966 // influence. however, in order to keep the
3967 // interface consistent with other functions,
3968 // we still include it as an argument to this
3969 // function
3970 const RefinementCase<dim - 1> std_face_ref =
3971 (face_orientation == face_rotation) ? flip[face_refinement_case] :
3972 face_refinement_case;
3973
3974 const RefinementCase<dim> face_to_cell[3][4] = {
3975 {RefinementCase<dim>::no_refinement, // faces 0 and 1
3976 RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3977 RefinementCase<dim>::cut_z,
3978 RefinementCase<dim>::cut_yz},
3979
3980 {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3981 // "exchanged on faces 2 and 3")
3982 RefinementCase<dim>::cut_z,
3983 RefinementCase<dim>::cut_x,
3984 RefinementCase<dim>::cut_xz},
3985
3986 {RefinementCase<dim>::no_refinement, // faces 4 and 5
3987 RefinementCase<dim>::cut_x,
3988 RefinementCase<dim>::cut_y,
3989 RefinementCase<dim>::cut_xy}};
3990
3991 return face_to_cell[face_no / 2][std_face_ref];
3992}
3993
3994
3995
3996template <int dim>
3997inline RefinementCase<dim>
3998GeometryInfo<dim>::min_cell_refinement_case_for_line_refinement(
3999 const unsigned int)
4000{
4001 DEAL_II_NOT_IMPLEMENTED();
4002
4003 return RefinementCase<dim>::no_refinement;
4004}
4005
4006template <>
4007inline RefinementCase<1>
4008GeometryInfo<1>::min_cell_refinement_case_for_line_refinement(
4009 const unsigned int line_no)
4010{
4011 (void)line_no;
4012 AssertIndexRange(line_no, 1);
4013
4014 return RefinementCase<1>::cut_x;
4015}
4016
4017
4018template <>
4019inline RefinementCase<2>
4020GeometryInfo<2>::min_cell_refinement_case_for_line_refinement(
4021 const unsigned int line_no)
4022{
4023 const unsigned int dim = 2;
4024 (void)dim;
4025 AssertIndexRange(line_no, GeometryInfo<dim>::lines_per_cell);
4026
4027 return (line_no / 2) != 0u ? RefinementCase<2>::cut_x :
4028 RefinementCase<2>::cut_y;
4029}
4030
4031
4032template <>
4033inline RefinementCase<3>
4034GeometryInfo<3>::min_cell_refinement_case_for_line_refinement(
4035 const unsigned int line_no)
4036{
4037 const unsigned int dim = 3;
4038 AssertIndexRange(line_no, GeometryInfo<dim>::lines_per_cell);
4039
4040 const RefinementCase<dim> ref_cases[6] = {
4041 RefinementCase<dim>::cut_y, // lines 0 and 1
4042 RefinementCase<dim>::cut_x, // lines 2 and 3
4043 RefinementCase<dim>::cut_y, // lines 4 and 5
4044 RefinementCase<dim>::cut_x, // lines 6 and 7
4045 RefinementCase<dim>::cut_z, // lines 8 and 9
4046 RefinementCase<dim>::cut_z}; // lines 10 and 11
4047
4048 return ref_cases[line_no / 2];
4049}
4050
4051
4052
4053template <>
4054inline unsigned int
4055GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
4056 const bool face_orientation,
4057 const bool face_flip,
4058 const bool face_rotation)
4059{
4060 AssertIndexRange(vertex, GeometryInfo<3>::vertices_per_face);
4061
4062 // set up a table to make sure that
4063 // we handle non-standard faces correctly
4064 //
4065 // so set up a table that for each vertex (of
4066 // a quad in standard position) describes
4067 // which vertex to take
4068 //
4069 // first index: four vertices 0...3
4070 //
4071 // second index: face_orientation; 0:
4072 // opposite normal, 1: standard
4073 //
4074 // third index: face_flip; 0: standard, 1:
4075 // face rotated by 180 degrees
4076 //
4077 // forth index: face_rotation: 0: standard,
4078 // 1: face rotated by 90 degrees
4079
4080 const unsigned int vertex_translation[4][2][2][2] = {
4081 {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
4082 // face_rotation=false and true
4083 {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
4084 // face_rotation=false and true
4085 {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
4086 // face_rotation=false and true
4087 {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
4088 // face_rotation=false and true
4089
4090 {{{2, 3}, // vertex 1 ...
4091 {1, 0}},
4092 {{1, 3}, {2, 0}}},
4093
4094 {{{1, 0}, // vertex 2 ...
4095 {2, 3}},
4096 {{2, 0}, {1, 3}}},
4097
4098 {{{3, 1}, // vertex 3 ...
4099 {0, 2}},
4100 {{3, 2}, {0, 1}}}};
4101
4102 return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
4103}
4104
4105
4106
4107template <int dim>
4108inline unsigned int
4109GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
4110 const bool,
4111 const bool,
4112 const bool)
4113{
4114 Assert(dim > 1, ExcImpossibleInDim(dim));
4115 AssertIndexRange(vertex, GeometryInfo<dim>::vertices_per_face);
4116 return vertex;
4117}
4118
4119
4120
4121template <>
4122inline unsigned int
4123GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
4124 const bool face_orientation,
4125 const bool face_flip,
4126 const bool face_rotation)
4127{
4128 AssertIndexRange(line, GeometryInfo<3>::lines_per_face);
4129
4130
4131 // make sure we handle
4132 // non-standard faces correctly
4133 //
4134 // so set up a table that for each line (of a
4135 // quad) describes which line to take
4136 //
4137 // first index: four lines 0...3
4138 //
4139 // second index: face_orientation; 0:
4140 // opposite normal, 1: standard
4141 //
4142 // third index: face_flip; 0: standard, 1:
4143 // face rotated by 180 degrees
4144 //
4145 // forth index: face_rotation: 0: standard,
4146 // 1: face rotated by 90 degrees
4147
4148 const unsigned int line_translation[4][2][2][2] = {
4149 {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4150 // face_rotation=false and true
4151 {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4152 // face_rotation=false and true
4153 {{0, 3}, // line 0, face_orientation=true, face_flip=false,
4154 // face_rotation=false and true
4155 {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
4156 // face_rotation=false and true
4157
4158 {{{3, 1}, // line 1 ...
4159 {2, 0}},
4160 {{1, 2}, {0, 3}}},
4161
4162 {{{0, 3}, // line 2 ...
4163 {1, 2}},
4164 {{2, 0}, {3, 1}}},
4165
4166 {{{1, 2}, // line 3 ...
4167 {0, 3}},
4168 {{3, 1}, {2, 0}}}};
4169
4170 return line_translation[line][face_orientation][face_flip][face_rotation];
4171}
4172
4173
4174
4175template <int dim>
4176inline unsigned int
4177GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
4178 const bool,
4179 const bool,
4180 const bool)
4181{
4182 DEAL_II_NOT_IMPLEMENTED();
4183 return line;
4184}
4185
4186
4187
4188template <>
4189inline unsigned int
4190GeometryInfo<2>::standard_to_real_line_vertex(const unsigned int vertex,
4191 const bool line_orientation)
4192{
4193 return line_orientation ? vertex : (1 - vertex);
4194}
4195
4196
4197
4198template <int dim>
4199inline unsigned int
4200GeometryInfo<dim>::standard_to_real_line_vertex(const unsigned int vertex,
4201 const bool)
4202{
4203 DEAL_II_NOT_IMPLEMENTED();
4204 return vertex;
4205}
4206
4207
4208
4209template <>
4210inline std::array<unsigned int, 2>
4211GeometryInfo<2>::standard_quad_vertex_to_line_vertex_index(
4212 const unsigned int vertex)
4213{
4214 return {{vertex % 2, vertex / 2}};
4215}
4216
4217
4218
4219template <int dim>
4220inline std::array<unsigned int, 2>
4221GeometryInfo<dim>::standard_quad_vertex_to_line_vertex_index(
4222 const unsigned int vertex)
4223{
4224 DEAL_II_NOT_IMPLEMENTED();
4225 (void)vertex;
4226 return {{0, 0}};
4227}
4228
4229
4230
4231template <>
4232inline std::array<unsigned int, 2>
4233GeometryInfo<3>::standard_hex_line_to_quad_line_index(const unsigned int i)
4234{
4235 // set up a table that for each
4236 // line describes a) from which
4237 // quad to take it, b) which line
4238 // therein it is if the face is
4239 // oriented correctly
4240 static const unsigned int lookup_table[GeometryInfo<3>::lines_per_cell][2] = {
4241 {4, 0}, // take first four lines from bottom face
4242 {4, 1},
4243 {4, 2},
4244 {4, 3},
4245
4246 {5, 0}, // second four lines from top face
4247 {5, 1},
4248 {5, 2},
4249 {5, 3},
4250
4251 {0, 0}, // the rest randomly
4252 {1, 0},
4253 {0, 1},
4254 {1, 1}};
4255
4256 return {{lookup_table[i][0], lookup_table[i][1]}};
4257}
4258
4259
4260
4261template <int dim>
4262inline std::array<unsigned int, 2>
4263GeometryInfo<dim>::standard_hex_line_to_quad_line_index(const unsigned int line)
4264{
4265 DEAL_II_NOT_IMPLEMENTED();
4266 (void)line;
4267 return {{0, 0}};
4268}
4269
4270
4271
4272template <>
4273inline std::array<unsigned int, 2>
4274GeometryInfo<3>::standard_hex_vertex_to_quad_vertex_index(
4275 const unsigned int vertex)
4276{
4277 // get the corner indices by asking either the bottom or the top face for its
4278 // vertices. handle non-standard faces by calling the vertex reordering
4279 // function from GeometryInfo
4280
4281 // bottom face (4) for first four vertices, top face (5) for the rest
4282 return {{4 + vertex / 4, vertex % 4}};
4283}
4284
4285
4286
4287template <int dim>
4288inline std::array<unsigned int, 2>
4289GeometryInfo<dim>::standard_hex_vertex_to_quad_vertex_index(
4290 const unsigned int vertex)
4291{
4292 DEAL_II_NOT_IMPLEMENTED();
4293 (void)vertex;
4294 return {{0, 0}};
4295}
4296
4297
4298
4299template <>
4300inline unsigned int
4301GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
4302 const bool face_orientation,
4303 const bool face_flip,
4304 const bool face_rotation)
4305{
4306 AssertIndexRange(line, GeometryInfo<3>::lines_per_face);
4307
4308
4309 // make sure we handle
4310 // non-standard faces correctly
4311 //
4312 // so set up a table that for each line (of a
4313 // quad) describes which line to take
4314 //
4315 // first index: four lines 0...3
4316 //
4317 // second index: face_orientation; 0:
4318 // opposite normal, 1: standard
4319 //
4320 // third index: face_flip; 0: standard, 1:
4321 // face rotated by 180 degrees
4322 //
4323 // forth index: face_rotation: 0: standard,
4324 // 1: face rotated by 90 degrees
4325
4326 const unsigned int line_translation[4][2][2][2] = {
4327 {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4328 // face_rotation=false and true
4329 {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4330 // face_rotation=false and true
4331 {{0, 2}, // line 0, face_orientation=true, face_flip=false,
4332 // face_rotation=false and true
4333 {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
4334 // face_rotation=false and true
4335
4336 {{{3, 1}, // line 1 ...
4337 {2, 0}},
4338 {{1, 3}, {0, 2}}},
4339
4340 {{{0, 3}, // line 2 ...
4341 {1, 2}},
4342 {{2, 1}, {3, 0}}},
4343
4344 {{{1, 2}, // line 3 ...
4345 {0, 3}},
4346 {{3, 0}, {2, 1}}}};
4347
4348 return line_translation[line][face_orientation][face_flip][face_rotation];
4349}
4350
4351
4352
4353template <int dim>
4354inline unsigned int
4355GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
4356 const bool,
4357 const bool,
4358 const bool)
4359{
4360 DEAL_II_NOT_IMPLEMENTED();
4361 return line;
4362}
4363
4364
4365
4366template <>
4367inline unsigned int
4368GeometryInfo<1>::child_cell_on_face(const RefinementCase<1> &,
4369 const unsigned int face,
4370 const unsigned int subface,
4371 const bool,
4372 const bool,
4373 const bool,
4374 const RefinementCase<0> &)
4375{
4376 (void)subface;
4377 AssertIndexRange(face, faces_per_cell);
4378 AssertIndexRange(subface, max_children_per_face);
4379
4380 return face;
4381}
4382
4383
4384
4385template <>
4386inline unsigned int
4387GeometryInfo<2>::child_cell_on_face(const RefinementCase<2> &ref_case,
4388 const unsigned int face,
4389 const unsigned int subface,
4390 const bool face_orientation,
4391 const bool /*face_flip*/,
4392 const bool /*face_rotation*/,
4393 const RefinementCase<1> &)
4394{
4395 AssertIndexRange(face, faces_per_cell);
4396 AssertIndexRange(subface, max_children_per_face);
4397
4398 // always return the child adjacent to the specified
4399 // subface. if the face of a cell is not refined, don't
4400 // throw an assertion but deliver the child adjacent to
4401 // the face nevertheless, i.e. deliver the child of
4402 // this cell adjacent to the subface of a possibly
4403 // refined neighbor. this simplifies setting neighbor
4404 // information in execute_refinement.
4405 constexpr unsigned int
4406 subcells[/* possible face orientation */ 2]
4407 [/* number of different ways to refine a cell */ 4]
4408 [/* faces_per_cell */ 4][/* max_children_per_face */ 2] = {
4409 {
4410 // Flipped orientation (face_orientation = false)
4411 {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
4412 {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
4413 {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
4414 },
4415 {
4416 // Normal orientation (face_orientation = true)
4417 {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
4418 {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
4419 {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
4420 }};
4421
4422 return subcells[face_orientation][ref_case - 1][face][subface];
4423}
4424
4425
4426
4427template <>
4428inline unsigned int
4429GeometryInfo<3>::child_cell_on_face(const RefinementCase<3> &ref_case,
4430 const unsigned int face,
4431 const unsigned int subface,
4432 const bool face_orientation,
4433 const bool face_flip,
4434 const bool face_rotation,
4435 const RefinementCase<2> &face_ref_case)
4436{
4437 const unsigned int dim = 3;
4438
4439 Assert(ref_case > RefinementCase<dim - 1>::no_refinement,
4440 ExcMessage("Cell has no children."));
4441 AssertIndexRange(face, faces_per_cell);
4442 if (!(subface == 0 &&
4443 face_ref_case == RefinementCase<dim - 1>::no_refinement))
4444 {
4445 AssertIndexRange(subface,
4446 GeometryInfo<dim - 1>::n_children(face_ref_case));
4447 }
4448
4449 // invalid number used for invalid cases,
4450 // e.g. when the children are more refined at
4451 // a given face than the face itself
4452 const unsigned int e = numbers::invalid_unsigned_int;
4453
4454 // the whole process of finding a child cell
4455 // at a given subface considering the
4456 // possibly anisotropic refinement cases of
4457 // the cell and the face as well as
4458 // orientation, flip and rotation of the face
4459 // is quite complicated. thus, we break it
4460 // down into several steps.
4461
4462 // first step: convert the given face refine
4463 // case to a face refine case concerning the
4464 // face in standard orientation (, flip and
4465 // rotation). This only affects cut_x and
4466 // cut_y
4467 const RefinementCase<dim - 1> flip[4] = {
4468 RefinementCase<dim - 1>::no_refinement,
4469 RefinementCase<dim - 1>::cut_y,
4470 RefinementCase<dim - 1>::cut_x,
4471 RefinementCase<dim - 1>::cut_xy};
4472 // for face_orientation, 'true' is the
4473 // default value whereas for face_rotation,
4474 // 'false' is standard. If
4475 // <tt>face_rotation==face_orientation</tt>,
4476 // then one of them is non-standard and we
4477 // have to swap cut_x and cut_y, otherwise no
4478 // change is necessary.
4479 const RefinementCase<dim - 1> std_face_ref =
4480 (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4481
4482 // second step: convert the given subface
4483 // index to the one for a standard face
4484 // respecting face_orientation, face_flip and
4485 // face_rotation
4486
4487 // first index: face_ref_case
4488 // second index: face_orientation
4489 // third index: face_flip
4490 // forth index: face_rotation
4491 // fifth index: subface index
4492 const unsigned int subface_exchange[4][2][2][2][4] = {
4493 // no_refinement (subface 0 stays 0,
4494 // all others are invalid)
4495 {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4496 {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4497 // cut_x (here, if the face is only
4498 // rotated OR only falsely oriented,
4499 // then subface 0 of the non-standard
4500 // face does NOT correspond to one of
4501 // the subfaces of a standard
4502 // face. Thus we indicate the subface
4503 // which is located at the lower left
4504 // corner (the origin of the face's
4505 // local coordinate system) with
4506 // '0'. The rest of this issue is
4507 // taken care of using the above
4508 // conversion to a 'standard face
4509 // refine case')
4510 {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4511 {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4512 // cut_y (the same applies as for
4513 // cut_x)
4514 {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4515 {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4516 // cut_xyz: this information is
4517 // identical to the information
4518 // returned by
4519 // GeometryInfo<3>::real_to_standard_face_vertex()
4520 {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4521 // face_rotation=false, subfaces 0,1,2,3
4522 {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4523 // face_rotation=true, subfaces 0,1,2,3
4524 {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4525 // face_rotation=false, subfaces 0,1,2,3
4526 {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4527 // face_rotation=true, subfaces 0,1,2,3
4528 {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4529 // face_rotation=false, subfaces 0,1,2,3
4530 {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4531 // face_rotation=true, subfaces 0,1,2,3
4532 {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4533 // face_rotation=false, subfaces 0,1,2,3
4534 {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4535 // face_rotation=true, subfaces 0,1,2,3
4536
4537 const unsigned int std_subface =
4538 subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4539 [subface];
4540 Assert(std_subface != e, ExcInternalError());
4541
4542 // third step: these are the children, which
4543 // can be found at the given subfaces of an
4544 // isotropically refined (standard) face
4545 //
4546 // first index: (refinement_case-1)
4547 // second index: face_index
4548 // third index: subface_index (isotropic refinement)
4549 const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4550 [max_children_per_face] = {
4551 // cut_x
4552 {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4553 {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4554 {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4555 {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4556 {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4557 {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4558 // cut_y
4559 {{0, 1, 0, 1},
4560 {0, 1, 0, 1},
4561 {0, 0, 0, 0},
4562 {1, 1, 1, 1},
4563 {0, 0, 1, 1},
4564 {0, 0, 1, 1}},
4565 // cut_xy
4566 {{0, 2, 0, 2},
4567 {1, 3, 1, 3},
4568 {0, 0, 1, 1},
4569 {2, 2, 3, 3},
4570 {0, 1, 2, 3},
4571 {0, 1, 2, 3}},
4572 // cut_z
4573 {{0, 0, 1, 1},
4574 {0, 0, 1, 1},
4575 {0, 1, 0, 1},
4576 {0, 1, 0, 1},
4577 {0, 0, 0, 0},
4578 {1, 1, 1, 1}},
4579 // cut_xz
4580 {{0, 0, 1, 1},
4581 {2, 2, 3, 3},
4582 {0, 1, 2, 3},
4583 {0, 1, 2, 3},
4584 {0, 2, 0, 2},
4585 {1, 3, 1, 3}},
4586 // cut_yz
4587 {{0, 1, 2, 3},
4588 {0, 1, 2, 3},
4589 {0, 2, 0, 2},
4590 {1, 3, 1, 3},
4591 {0, 0, 1, 1},
4592 {2, 2, 3, 3}},
4593 // cut_xyz
4594 {{0, 2, 4, 6},
4595 {1, 3, 5, 7},
4596 {0, 4, 1, 5},
4597 {2, 6, 3, 7},
4598 {0, 1, 2, 3},
4599 {4, 5, 6, 7}}};
4600
4601 // forth step: check, whether the given face
4602 // refine case is valid for the given cell
4603 // refine case. this is the case, if the
4604 // given face refine case is at least as
4605 // refined as the face is for the given cell
4606 // refine case
4607
4608 // note, that we are considering standard
4609 // face refinement cases here and thus must
4610 // not pass the given orientation, flip and
4611 // rotation flags
4612 if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4613 face_refinement_case(ref_case, face))
4614 {
4615 // all is fine. for anisotropic face
4616 // refine cases, select one of the
4617 // isotropic subfaces which neighbors the
4618 // same child
4619
4620 // first index: (standard) face refine case
4621 // second index: subface index
4622 const unsigned int equivalent_iso_subface[4][4] = {
4623 {0, e, e, e}, // no_refinement
4624 {0, 3, e, e}, // cut_x
4625 {0, 3, e, e}, // cut_y
4626 {0, 1, 2, 3}}; // cut_xy
4627
4628 const unsigned int equ_std_subface =
4629 equivalent_iso_subface[std_face_ref][std_subface];
4630 Assert(equ_std_subface != e, ExcInternalError());
4631
4632 return iso_children[ref_case - 1][face][equ_std_subface];
4633 }
4634 else
4635 {
4636 // the face_ref_case was too coarse,
4637 // throw an error
4638 Assert(false,
4639 ExcMessage("The face RefineCase is too coarse "
4640 "for the given cell RefineCase."));
4641 }
4642 // we only get here in case of an error
4643 return e;
4644}
4645
4646
4647
4648template <>
4649inline unsigned int
4650GeometryInfo<4>::child_cell_on_face(const RefinementCase<4> &,
4651 const unsigned int,
4652 const unsigned int,
4653 const bool,
4654 const bool,
4655 const bool,
4656 const RefinementCase<3> &)
4657{
4658 DEAL_II_NOT_IMPLEMENTED();
4659 return numbers::invalid_unsigned_int;
4660}
4661
4662
4663
4664template <>
4665inline unsigned int
4666GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4667 const unsigned int line,
4668 const bool,
4669 const bool,
4670 const bool)
4671{
4672 (void)line;
4673 AssertIndexRange(face, faces_per_cell);
4674 AssertIndexRange(line, lines_per_face);
4675
4676 // The face is a line itself.
4677 return face;
4678}
4679
4680
4681
4682template <>
4683inline unsigned int
4684GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4685 const unsigned int line,
4686 const bool face_orientation,
4687 const bool face_flip,
4688 const bool face_rotation)
4689{
4690 AssertIndexRange(face, faces_per_cell);
4691 AssertIndexRange(line, lines_per_face);
4692
4693 const unsigned lines[faces_per_cell][lines_per_face] = {
4694 {8, 10, 0, 4}, // left face
4695 {9, 11, 1, 5}, // right face
4696 {2, 6, 8, 9}, // front face
4697 {3, 7, 10, 11}, // back face
4698 {0, 1, 2, 3}, // bottom face
4699 {4, 5, 6, 7}}; // top face
4700 return lines[face][real_to_standard_face_line(
4701 line, face_orientation, face_flip, face_rotation)];
4702}
4703
4704
4705
4706inline unsigned int
4707GeometryInfo<0>::face_to_cell_lines(const unsigned int,
4708 const unsigned int,
4709 const bool,
4710 const bool,
4711 const bool)
4712{
4713 DEAL_II_NOT_IMPLEMENTED();
4714 return numbers::invalid_unsigned_int;
4715}
4716
4717
4718
4719template <int dim>
4720inline unsigned int
4721GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4722 const unsigned int,
4723 const bool,
4724 const bool,
4725 const bool)
4726{
4727 DEAL_II_NOT_IMPLEMENTED();
4728 return numbers::invalid_unsigned_int;
4729}
4730
4731
4732
4733template <int dim>
4734inline unsigned int
4735GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4736 const unsigned int vertex,
4737 const bool face_orientation,
4738 const bool face_flip,
4739 const bool face_rotation)
4740{
4741 return child_cell_on_face(RefinementCase<dim>::isotropic_refinement,
4742 face,
4743 vertex,
4744 face_orientation,
4745 face_flip,
4746 face_rotation);
4747}
4748
4749
4750
4751inline unsigned int
4752GeometryInfo<0>::face_to_cell_vertices(const unsigned int,
4753 const unsigned int,
4754 const bool,
4755 const bool,
4756 const bool)
4757{
4758 DEAL_II_NOT_IMPLEMENTED();
4759 return numbers::invalid_unsigned_int;
4760}
4761
4762
4763
4764template <int dim>
4765template <typename Number>
4766inline Point<dim, Number>
4767GeometryInfo<dim>::project_to_unit_cell(const Point<dim, Number> &q)
4768{
4769 Point<dim, Number> p;
4770 for (unsigned int i = 0; i < dim; ++i)
4771 p[i] = std::min(std::max(q[i], Number(0.)), Number(1.));
4772
4773 return p;
4774}
4775
4776
4777
4778template <int dim>
4779inline double
4780GeometryInfo<dim>::distance_to_unit_cell(const Point<dim> &p)
4781{
4782 double result = 0.0;
4783
4784 for (unsigned int i = 0; i < dim; ++i)
4785 {
4786 result = std::max(result, -p[i]);
4787 result = std::max(result, p[i] - 1.);
4788 }
4789
4790 return result;
4791}
4792
4793
4794
4795template <int dim>
4796inline double
4797GeometryInfo<dim>::d_linear_shape_function(const Point<dim> &xi,
4798 const unsigned int i)
4799{
4800 AssertIndexRange(i, GeometryInfo<dim>::vertices_per_cell);
4801
4802 switch (dim)
4803 {
4804 case 1:
4805 {
4806 const double x = xi[0];
4807 switch (i)
4808 {
4809 case 0:
4810 return 1 - x;
4811 case 1:
4812 return x;
4813 }
4814 break;
4815 }
4816
4817 case 2:
4818 {
4819 const double x = xi[0];
4820 const double y = xi[1];
4821 switch (i)
4822 {
4823 case 0:
4824 return (1 - x) * (1 - y);
4825 case 1:
4826 return x * (1 - y);
4827 case 2:
4828 return (1 - x) * y;
4829 case 3:
4830 return x * y;
4831 }
4832 break;
4833 }
4834
4835 case 3:
4836 {
4837 const double x = xi[0];
4838 const double y = xi[1];
4839 const double z = xi[2];
4840 switch (i)
4841 {
4842 case 0:
4843 return (1 - x) * (1 - y) * (1 - z);
4844 case 1:
4845 return x * (1 - y) * (1 - z);
4846 case 2:
4847 return (1 - x) * y * (1 - z);
4848 case 3:
4849 return x * y * (1 - z);
4850 case 4:
4851 return (1 - x) * (1 - y) * z;
4852 case 5:
4853 return x * (1 - y) * z;
4854 case 6:
4855 return (1 - x) * y * z;
4856 case 7:
4857 return x * y * z;
4858 }
4859 break;
4860 }
4861
4862 default:
4863 DEAL_II_NOT_IMPLEMENTED();
4864 }
4865 return -1e9;
4866}
4867
4868
4869
4870template <>
4871Tensor<1, 1> inline GeometryInfo<1>::d_linear_shape_function_gradient(
4872 const Point<1> &,
4873 const unsigned int i)
4874{
4875 AssertIndexRange(i, GeometryInfo<1>::vertices_per_cell);
4876
4877 switch (i)
4878 {
4879 case 0:
4880 return Point<1>(-1.);
4881 case 1:
4882 return Point<1>(1.);
4883 }
4884
4885 return Point<1>(-1e9);
4886}
4887
4888
4889
4890template <>
4891Tensor<1, 2> inline GeometryInfo<2>::d_linear_shape_function_gradient(
4892 const Point<2> &xi,
4893 const unsigned int i)
4894{
4895 AssertIndexRange(i, GeometryInfo<2>::vertices_per_cell);
4896
4897 const double x = xi[0];
4898 const double y = xi[1];
4899 switch (i)
4900 {
4901 case 0:
4902 return Point<2>(-(1 - y), -(1 - x));
4903 case 1:
4904 return Point<2>(1 - y, -x);
4905 case 2:
4906 return Point<2>(-y, 1 - x);
4907 case 3:
4908 return Point<2>(y, x);
4909 }
4910 return Point<2>(-1e9, -1e9);
4911}
4912
4913
4914
4915template <>
4916Tensor<1, 3> inline GeometryInfo<3>::d_linear_shape_function_gradient(
4917 const Point<3> &xi,
4918 const unsigned int i)
4919{
4920 AssertIndexRange(i, GeometryInfo<3>::vertices_per_cell);
4921
4922 const double x = xi[0];
4923 const double y = xi[1];
4924 const double z = xi[2];
4925 switch (i)
4926 {
4927 case 0:
4928 return Point<3>(-(1 - y) * (1 - z),
4929 -(1 - x) * (1 - z),
4930 -(1 - x) * (1 - y));
4931 case 1:
4932 return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4933 case 2:
4934 return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4935 case 3:
4936 return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4937 case 4:
4938 return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4939 case 5:
4940 return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4941 case 6:
4942 return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4943 case 7:
4944 return Point<3>(y * z, x * z, x * y);
4945 }
4946
4947 return Point<3>(-1e9, -1e9, -1e9);
4948}
4949
4950
4951
4952template <int dim>
4953inline Tensor<1, dim>
4954GeometryInfo<dim>::d_linear_shape_function_gradient(const Point<dim> &,
4955 const unsigned int)
4956{
4957 DEAL_II_NOT_IMPLEMENTED();
4958 return Tensor<1, dim>();
4959}
4960
4961
4962unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4963{
4964 return 0;
4965}
4966
4967
4968namespace internal
4969{
4970 namespace GeometryInfoHelper
4971 {
4972 // wedge product of a single
4973 // vector in 2d: we just have to
4974 // rotate it by 90 degrees to the
4975 // right
4976 inline Tensor<1, 2>
4977 wedge_product(const Tensor<1, 2> (&derivative)[1])
4978 {
4979 Tensor<1, 2> result;
4980 result[0] = derivative[0][1];
4981 result[1] = -derivative[0][0];
4982
4983 return result;
4984 }
4985
4986
4987 // wedge product of 2 vectors in
4988 // 3d is the cross product
4989 inline Tensor<1, 3>
4990 wedge_product(const Tensor<1, 3> (&derivative)[2])
4991 {
4992 return cross_product_3d(derivative[0], derivative[1]);
4993 }
4994
4995
4996 // wedge product of dim vectors
4997 // in dim-d: that's the
4998 // determinant of the matrix
4999 template <int dim>
5000 inline Tensor<0, dim>
5001 wedge_product(const Tensor<1, dim> (&derivative)[dim])
5002 {
5003 Tensor<2, dim> jacobian;
5004 for (unsigned int i = 0; i < dim; ++i)
5005 jacobian[i] = derivative[i];
5006
5007 return determinant(jacobian);
5008 }
5009 } // namespace GeometryInfoHelper
5010} // namespace internal
5011
5012
5013template <int dim>
5014template <int spacedim>
5015inline void
5016GeometryInfo<dim>::alternating_form_at_vertices
5017# ifndef DEAL_II_CXX14_CONSTEXPR_BUG
5018 (const Point<spacedim> (&vertices)[vertices_per_cell],
5019 Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell])
5020# else
5021 (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms)
5022# endif
5023{
5024 // for each of the vertices,
5025 // compute the alternating form
5026 // of the mapped unit
5027 // vectors. consider for
5028 // example the case of a quad
5029 // in spacedim==3: to do so, we
5030 // need to see how the
5031 // infinitesimal vectors
5032 // (d\xi_1,0) and (0,d\xi_2)
5033 // are transformed into
5034 // spacedim-dimensional space
5035 // and then form their cross
5036 // product (i.e. the wedge product
5037 // of two vectors). to this end, note
5038 // that
5039 // \vec x = sum_i \vec v_i phi_i(\vec xi)
5040 // so the transformed vectors are
5041 // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
5042 // and
5043 // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
5044 // which boils down to the columns
5045 // of the 3x2 matrix \grad_\xi x(\xi)
5046 //
5047 // a similar reasoning would
5048 // hold for all dim,spacedim
5049 // pairs -- we only have to
5050 // compute the wedge product of
5051 // the columns of the
5052 // derivatives
5053 for (unsigned int i = 0; i < vertices_per_cell; ++i)
5054 {
5055 Tensor<1, spacedim> derivatives[dim];
5056
5057 for (unsigned int j = 0; j < vertices_per_cell; ++j)
5058 {
5059 const Tensor<1, dim> grad_phi_j =
5060 d_linear_shape_function_gradient(unit_cell_vertex(i), j);
5061 for (unsigned int l = 0; l < dim; ++l)
5062 derivatives[l] += vertices[j] * grad_phi_j[l];
5063 }
5064
5065 forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
5066 }
5067}
5068
5069#endif // DOXYGEN
5070DEAL_II_NAMESPACE_CLOSE
5071
5072#endif
GeometryPrimitive::GeometryPrimitive
GeometryPrimitive(const unsigned int object_dimension)
GeometryPrimitive::GeometryPrimitive
GeometryPrimitive(const Object object)
Point
Definition point.h:111
RefinementCase::operator|
RefinementCase operator|(const RefinementCase &r) const
RefinementCase::memory_consumption
static std::size_t memory_consumption()
RefinementCase::operator~
RefinementCase operator~() const
RefinementCase::operator&
RefinementCase operator&(const RefinementCase &r) const
RefinementCase::RefinementCase
RefinementCase(const typename RefinementPossibilities< dim >::Possibilities refinement_case)
RefinementCase::n_refinement_cases
static constexpr unsigned int n_refinement_cases
Definition geometry_info.h:784
RefinementCase::cut_axis
static RefinementCase cut_axis(const unsigned int i)
RefinementCase::value
std::uint8_t value
Definition geometry_info.h:892
RefinementCase::all_refinement_cases
static std::array< RefinementCase< dim >, n_refinement_cases > all_refinement_cases()
RefinementCase::serialize
void serialize(Archive &ar, const unsigned int version)
RefinementCase::RefinementCase
RefinementCase(const std::uint8_t refinement_case)
internal::SubfaceCase::memory_consumption
static constexpr std::size_t memory_consumption()
internal::SubfaceCase::SubfaceCase
SubfaceCase(const typename SubfacePossibilities< dim >::Possibilities subface_possibility)
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
Assert
#define Assert(cond, exc)
Definition exceptions.h:1659
StandardExceptions::ExcImpossibleInDim
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
GeometryInfo::ExcInvalidCoordinate
static ::ExceptionBase & ExcInvalidCoordinate(double arg1)
GeometryInfo::ExcInvalidSubface
static ::ExceptionBase & ExcInvalidSubface(int arg1, int arg2, int arg3)
RefinementCase::ExcInvalidRefinementCase
static ::ExceptionBase & ExcInvalidRefinementCase(int arg1)
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:2115
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
DeclException3
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition exceptions.h:557
DeclException1
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:511
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
internal::SubfaceCase::ExcInvalidSubfaceCase
static ::ExceptionBase & ExcInvalidSubfaceCase(int arg1)
DEAL_II_ASSERT_UNREACHABLE
#define DEAL_II_ASSERT_UNREACHABLE()
Definition exceptions.h:1938
DEAL_II_NOT_IMPLEMENTED
#define DEAL_II_NOT_IMPLEMENTED()
Definition exceptions.h:1855
iota_view.h
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition types.h:220
std_cxx20::ranges::iota_view
boost::integer_range< IncrementableType > iota_view
Definition iota_view.h:45
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6960
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6943
ndarray
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:107
DEAL_II_HOST_DEVICE
#define DEAL_II_HOST_DEVICE
Definition numbers.h:30
GeometryInfo< 0 >::n_children
static unsigned int n_children(const RefinementCase< 0 > &refinement_case)
GeometryInfo< 0 >::face_to_cell_vertices
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)
GeometryInfo< 0 >::face_indices
static std::array< unsigned int, 0 > face_indices()
GeometryInfo< 0 >::vertex_indices
static std::array< unsigned int, vertices_per_cell > vertex_indices()
GeometryInfo< 0 >::face_to_cell_lines
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)
GeometryInfo::child_cell_from_point
static unsigned int child_cell_from_point(const Point< dim > &p)
GeometryInfo::standard_to_real_line_vertex
static unsigned int standard_to_real_line_vertex(const unsigned int vertex, const bool line_orientation=true)
GeometryInfo::ucd_to_deal
static constexpr std::array< unsigned int, vertices_per_cell > ucd_to_deal
Definition geometry_info.h:2097
GeometryInfo::unit_tangential_vectors
static constexpr ndarray< Tensor< 1, dim >, faces_per_cell, dim - 1 > unit_tangential_vectors
Definition geometry_info.h:2629
GeometryInfo::is_inside_unit_cell
static bool is_inside_unit_cell(const Point< dim > &p, const double eps)
GeometryInfo::distance_to_unit_cell
static double distance_to_unit_cell(const Point< dim > &p)
GeometryInfo::child_cell_on_face
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)
GeometryInfo::quads_per_face
static constexpr unsigned int quads_per_face
Definition geometry_info.h:2046
GeometryInfo::face_refinement_case
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)
GeometryInfo::standard_to_real_face_vertex
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)
GeometryInfo::min_cell_refinement_case_for_line_refinement
static RefinementCase< dim > min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
GeometryInfo::opposite_face
static constexpr std::array< unsigned int, faces_per_cell > opposite_face
Definition geometry_info.h:2637
GeometryInfo::max_children_per_face
static constexpr unsigned int max_children_per_face
Definition geometry_info.h:2005
GeometryInfo::face_to_cell_vertices
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)
GeometryInfo::vertices_per_cell
static constexpr unsigned int vertices_per_cell
Definition geometry_info.h:2011
GeometryInfo::dx_to_deal
static constexpr std::array< unsigned int, vertices_per_cell > dx_to_deal
Definition geometry_info.h:2113
GeometryInfo::lines_per_cell
static constexpr unsigned int lines_per_cell
Definition geometry_info.h:2058
GeometryInfo::unit_normal_direction
static constexpr std::array< unsigned int, faces_per_cell > unit_normal_direction
Definition geometry_info.h:2579
GeometryInfo::standard_hex_line_to_quad_line_index
static std::array< unsigned int, 2 > standard_hex_line_to_quad_line_index(const unsigned int line)
GeometryInfo::faces_per_cell
static constexpr unsigned int faces_per_cell
Definition geometry_info.h:1977
GeometryInfo::standard_to_real_face_line
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)
GeometryInfo::face_indices
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
GeometryInfo::subface_ratio
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
GeometryInfo::hexes_per_cell
static constexpr unsigned int hexes_per_cell
Definition geometry_info.h:2076
GeometryInfo::project_to_unit_cell
static Point< dim, Number > project_to_unit_cell(const Point< dim, Number > &p)
GeometryInfo::line_refinement_case
static RefinementCase< 1 > line_refinement_case(const RefinementCase< dim > &cell_refinement_case, const unsigned int line_no)
GeometryInfo::standard_hex_vertex_to_quad_vertex_index
static std::array< unsigned int, 2 > standard_hex_vertex_to_quad_vertex_index(const unsigned int vertex)
GeometryInfo::standard_quad_vertex_to_line_vertex_index
static std::array< unsigned int, 2 > standard_quad_vertex_to_line_vertex_index(const unsigned int vertex)
GeometryInfo::d_linear_shape_function
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
GeometryInfo::unit_normal_vector
static constexpr std::array< Tensor< 1, dim >, faces_per_cell > unit_normal_vector
Definition geometry_info.h:2612
GeometryInfo::vertices_per_face
static constexpr unsigned int vertices_per_face
Definition geometry_info.h:2034
GeometryInfo::real_to_standard_face_line
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)
GeometryInfo::line_to_cell_vertices
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
GeometryInfo::vertex_indices
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
GeometryInfo::is_inside_unit_cell
static bool is_inside_unit_cell(const Point< dim > &p)
GeometryInfo::min_cell_refinement_case_for_face_refinement
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)
GeometryInfo::n_children
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
GeometryInfo::real_to_standard_face_vertex
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)
GeometryInfo::face_to_cell_lines
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)
GeometryInfo::cell_to_child_coordinates
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)
GeometryInfo::vertex_to_face
static constexpr ndarray< unsigned int, vertices_per_cell, dim > vertex_to_face
Definition geometry_info.h:2127
GeometryInfo::d_linear_shape_function_gradient
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
GeometryInfo::alternating_form_at_vertices
static void alternating_form_at_vertices(const Point< spacedim >(&vertices)[vertices_per_cell], Tensor< spacedim - dim, spacedim >(&forms)[vertices_per_cell])
GeometryInfo::unit_cell_vertex
static Point< dim > unit_cell_vertex(const unsigned int vertex)
GeometryInfo::unit_normal_orientation
static constexpr std::array< int, faces_per_cell > unit_normal_orientation
Definition geometry_info.h:2598
GeometryInfo::quads_per_cell
static constexpr unsigned int quads_per_cell
Definition geometry_info.h:2069
GeometryInfo::child_to_cell_coordinates
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)
GeometryInfo::max_children_per_cell
static constexpr unsigned int max_children_per_cell
Definition geometry_info.h:1972
GeometryInfo::n_subfaces
static unsigned int n_subfaces(const internal::SubfaceCase< dim > &subface_case)
GeometryInfo::lines_per_face
static constexpr unsigned int lines_per_face
Definition geometry_info.h:2040
determinant
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
Definition symmetric_tensor.h:2639
Generated by doxygen 1.9.8