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