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