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