Reference documentation for deal.II version Git d3aed38b93 2021-10-28 13:33:27 +0200
\(\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]]) ?
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) ? 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) ? RefinementCase<2>::cut_x : RefinementCase<2>::cut_y;
3979 }
3980 
3981 
3982 template <>
3983 inline RefinementCase<3>
3985  const unsigned int line_no)
3986 {
3987  const unsigned int dim = 3;
3989 
3990  const RefinementCase<dim> ref_cases[6] = {
3991  RefinementCase<dim>::cut_y, // lines 0 and 1
3992  RefinementCase<dim>::cut_x, // lines 2 and 3
3993  RefinementCase<dim>::cut_y, // lines 4 and 5
3994  RefinementCase<dim>::cut_x, // lines 6 and 7
3995  RefinementCase<dim>::cut_z, // lines 8 and 9
3996  RefinementCase<dim>::cut_z}; // lines 10 and 11
3997 
3998  return ref_cases[line_no / 2];
3999 }
4000 
4001 
4002 
4003 template <>
4004 inline unsigned int
4005 GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
4006  const bool face_orientation,
4007  const bool face_flip,
4008  const bool face_rotation)
4009 {
4011 
4012  // set up a table to make sure that
4013  // we handle non-standard faces correctly
4014  //
4015  // so set up a table that for each vertex (of
4016  // a quad in standard position) describes
4017  // which vertex to take
4018  //
4019  // first index: four vertices 0...3
4020  //
4021  // second index: face_orientation; 0:
4022  // opposite normal, 1: standard
4023  //
4024  // third index: face_flip; 0: standard, 1:
4025  // face rotated by 180 degrees
4026  //
4027  // forth index: face_rotation: 0: standard,
4028  // 1: face rotated by 90 degrees
4029 
4030  const unsigned int vertex_translation[4][2][2][2] = {
4031  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
4032  // face_rotation=false and true
4033  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
4034  // face_rotation=false and true
4035  {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
4036  // face_rotation=false and true
4037  {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
4038  // face_rotation=false and true
4039 
4040  {{{2, 3}, // vertex 1 ...
4041  {1, 0}},
4042  {{1, 3}, {2, 0}}},
4043 
4044  {{{1, 0}, // vertex 2 ...
4045  {2, 3}},
4046  {{2, 0}, {1, 3}}},
4047 
4048  {{{3, 1}, // vertex 3 ...
4049  {0, 2}},
4050  {{3, 2}, {0, 1}}}};
4051 
4052  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
4053 }
4054 
4055 
4056 
4057 template <int dim>
4058 inline unsigned int
4059 GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
4060  const bool,
4061  const bool,
4062  const bool)
4063 {
4064  Assert(dim > 1, ExcImpossibleInDim(dim));
4066  return vertex;
4067 }
4068 
4069 
4070 
4071 template <>
4072 inline unsigned int
4073 GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
4074  const bool face_orientation,
4075  const bool face_flip,
4076  const bool face_rotation)
4077 {
4079 
4080 
4081  // make sure we handle
4082  // non-standard faces correctly
4083  //
4084  // so set up a table that for each line (of a
4085  // quad) describes which line to take
4086  //
4087  // first index: four lines 0...3
4088  //
4089  // second index: face_orientation; 0:
4090  // opposite normal, 1: standard
4091  //
4092  // third index: face_flip; 0: standard, 1:
4093  // face rotated by 180 degrees
4094  //
4095  // forth index: face_rotation: 0: standard,
4096  // 1: face rotated by 90 degrees
4097 
4098  const unsigned int line_translation[4][2][2][2] = {
4099  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4100  // face_rotation=false and true
4101  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4102  // face_rotation=false and true
4103  {{0, 3}, // line 0, face_orientation=true, face_flip=false,
4104  // face_rotation=false and true
4105  {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
4106  // face_rotation=false and true
4107 
4108  {{{3, 1}, // line 1 ...
4109  {2, 0}},
4110  {{1, 2}, {0, 3}}},
4111 
4112  {{{0, 3}, // line 2 ...
4113  {1, 2}},
4114  {{2, 0}, {3, 1}}},
4115 
4116  {{{1, 2}, // line 3 ...
4117  {0, 3}},
4118  {{3, 1}, {2, 0}}}};
4119 
4120  return line_translation[line][face_orientation][face_flip][face_rotation];
4121 }
4122 
4123 
4124 
4125 template <int dim>
4126 inline unsigned int
4127 GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
4128  const bool,
4129  const bool,
4130  const bool)
4131 {
4132  Assert(false, ExcNotImplemented());
4133  return line;
4134 }
4135 
4136 
4137 
4138 template <>
4139 inline unsigned int
4140 GeometryInfo<2>::standard_to_real_line_vertex(const unsigned int vertex,
4141  const bool line_orientation)
4142 {
4143  return line_orientation ? vertex : (1 - vertex);
4144 }
4145 
4146 
4147 
4148 template <int dim>
4149 inline unsigned int
4150 GeometryInfo<dim>::standard_to_real_line_vertex(const unsigned int vertex,
4151  const bool)
4152 {
4153  Assert(false, ExcNotImplemented());
4154  return vertex;
4155 }
4156 
4157 
4158 
4159 template <>
4160 inline std::array<unsigned int, 2>
4162  const unsigned int vertex)
4163 {
4164  return {{vertex % 2, vertex / 2}};
4165 }
4166 
4167 
4168 
4169 template <int dim>
4170 inline std::array<unsigned int, 2>
4172  const unsigned int vertex)
4173 {
4174  Assert(false, ExcNotImplemented());
4175  (void)vertex;
4176  return {{0, 0}};
4177 }
4178 
4179 
4180 
4181 template <>
4182 inline std::array<unsigned int, 2>
4184 {
4185  // set up a table that for each
4186  // line describes a) from which
4187  // quad to take it, b) which line
4188  // therein it is if the face is
4189  // oriented correctly
4190  static const unsigned int lookup_table[GeometryInfo<3>::lines_per_cell][2] = {
4191  {4, 0}, // take first four lines from bottom face
4192  {4, 1},
4193  {4, 2},
4194  {4, 3},
4195 
4196  {5, 0}, // second four lines from top face
4197  {5, 1},
4198  {5, 2},
4199  {5, 3},
4200 
4201  {0, 0}, // the rest randomly
4202  {1, 0},
4203  {0, 1},
4204  {1, 1}};
4205 
4206  return {{lookup_table[i][0], lookup_table[i][1]}};
4207 }
4208 
4209 
4210 
4211 template <int dim>
4212 inline std::array<unsigned int, 2>
4214 {
4215  Assert(false, ExcNotImplemented());
4216  (void)line;
4217  return {{0, 0}};
4218 }
4219 
4220 
4221 
4222 template <>
4223 inline std::array<unsigned int, 2>
4225  const unsigned int vertex)
4226 {
4227  // get the corner indices by asking either the bottom or the top face for its
4228  // vertices. handle non-standard faces by calling the vertex reordering
4229  // function from GeometryInfo
4230 
4231  // bottom face (4) for first four vertices, top face (5) for the rest
4232  return {{4 + vertex / 4, vertex % 4}};
4233 }
4234 
4235 
4236 
4237 template <int dim>
4238 inline std::array<unsigned int, 2>
4240  const unsigned int vertex)
4241 {
4242  Assert(false, ExcNotImplemented());
4243  (void)vertex;
4244  return {{0, 0}};
4245 }
4246 
4247 
4248 
4249 template <>
4250 inline unsigned int
4251 GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
4252  const bool face_orientation,
4253  const bool face_flip,
4254  const bool face_rotation)
4255 {
4257 
4258 
4259  // make sure we handle
4260  // non-standard faces correctly
4261  //
4262  // so set up a table that for each line (of a
4263  // quad) describes which line to take
4264  //
4265  // first index: four lines 0...3
4266  //
4267  // second index: face_orientation; 0:
4268  // opposite normal, 1: standard
4269  //
4270  // third index: face_flip; 0: standard, 1:
4271  // face rotated by 180 degrees
4272  //
4273  // forth index: face_rotation: 0: standard,
4274  // 1: face rotated by 90 degrees
4275 
4276  const unsigned int line_translation[4][2][2][2] = {
4277  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4278  // face_rotation=false and true
4279  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4280  // face_rotation=false and true
4281  {{0, 2}, // line 0, face_orientation=true, face_flip=false,
4282  // face_rotation=false and true
4283  {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
4284  // face_rotation=false and true
4285 
4286  {{{3, 1}, // line 1 ...
4287  {2, 0}},
4288  {{1, 3}, {0, 2}}},
4289 
4290  {{{0, 3}, // line 2 ...
4291  {1, 2}},
4292  {{2, 1}, {3, 0}}},
4293 
4294  {{{1, 2}, // line 3 ...
4295  {0, 3}},
4296  {{3, 0}, {2, 1}}}};
4297 
4298  return line_translation[line][face_orientation][face_flip][face_rotation];
4299 }
4300 
4301 
4302 
4303 template <int dim>
4304 inline unsigned int
4305 GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
4306  const bool,
4307  const bool,
4308  const bool)
4309 {
4310  Assert(false, ExcNotImplemented());
4311  return line;
4312 }
4313 
4314 
4315 
4316 template <>
4317 inline unsigned int
4319  const unsigned int face,
4320  const unsigned int subface,
4321  const bool,
4322  const bool,
4323  const bool,
4324  const RefinementCase<0> &)
4325 {
4326  (void)subface;
4329 
4330  return face;
4331 }
4332 
4333 
4334 
4335 template <>
4336 inline unsigned int
4338  const unsigned int face,
4339  const unsigned int subface,
4340  const bool /*face_orientation*/,
4341  const bool face_flip,
4342  const bool /*face_rotation*/,
4343  const RefinementCase<1> &)
4344 {
4347 
4348  // always return the child adjacent to the specified
4349  // subface. if the face of a cell is not refined, don't
4350  // throw an assertion but deliver the child adjacent to
4351  // the face nevertheless, i.e. deliver the child of
4352  // this cell adjacent to the subface of a possibly
4353  // refined neighbor. this simplifies setting neighbor
4354  // information in execute_refinement.
4355  const unsigned int
4357  [max_children_per_face] = {
4358  {
4359  // Normal orientation (face_flip = false)
4360  {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
4361  {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
4362  {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
4363  },
4364  {
4365  // Flipped orientation (face_flip = true)
4366  {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
4367  {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
4368  {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
4369  }};
4370 
4371  return subcells[face_flip][ref_case - 1][face][subface];
4372 }
4373 
4374 
4375 
4376 template <>
4377 inline unsigned int
4379  const unsigned int face,
4380  const unsigned int subface,
4381  const bool face_orientation,
4382  const bool face_flip,
4383  const bool face_rotation,
4384  const RefinementCase<2> &face_ref_case)
4385 {
4386  const unsigned int dim = 3;
4387 
4389  ExcMessage("Cell has no children."));
4391  if (!(subface == 0 &&
4392  face_ref_case == RefinementCase<dim - 1>::no_refinement))
4393  {
4394  AssertIndexRange(subface,
4395  GeometryInfo<dim - 1>::n_children(face_ref_case));
4396  }
4397 
4398  // invalid number used for invalid cases,
4399  // e.g. when the children are more refined at
4400  // a given face than the face itself
4401  const unsigned int e = numbers::invalid_unsigned_int;
4402 
4403  // the whole process of finding a child cell
4404  // at a given subface considering the
4405  // possibly anisotropic refinement cases of
4406  // the cell and the face as well as
4407  // orientation, flip and rotation of the face
4408  // is quite complicated. thus, we break it
4409  // down into several steps.
4410 
4411  // first step: convert the given face refine
4412  // case to a face refine case concerning the
4413  // face in standard orientation (, flip and
4414  // rotation). This only affects cut_x and
4415  // cut_y
4416  const RefinementCase<dim - 1> flip[4] = {
4417  RefinementCase<dim - 1>::no_refinement,
4418  RefinementCase<dim - 1>::cut_y,
4419  RefinementCase<dim - 1>::cut_x,
4420  RefinementCase<dim - 1>::cut_xy};
4421  // for face_orientation, 'true' is the
4422  // default value whereas for face_rotation,
4423  // 'false' is standard. If
4424  // <tt>face_rotation==face_orientation</tt>,
4425  // then one of them is non-standard and we
4426  // have to swap cut_x and cut_y, otherwise no
4427  // change is necessary.
4428  const RefinementCase<dim - 1> std_face_ref =
4429  (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4430 
4431  // second step: convert the given subface
4432  // index to the one for a standard face
4433  // respecting face_orientation, face_flip and
4434  // face_rotation
4435 
4436  // first index: face_ref_case
4437  // second index: face_orientation
4438  // third index: face_flip
4439  // forth index: face_rotation
4440  // fifth index: subface index
4441  const unsigned int subface_exchange[4][2][2][2][4] = {
4442  // no_refinement (subface 0 stays 0,
4443  // all others are invalid)
4444  {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4445  {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4446  // cut_x (here, if the face is only
4447  // rotated OR only falsely oriented,
4448  // then subface 0 of the non-standard
4449  // face does NOT correspond to one of
4450  // the subfaces of a standard
4451  // face. Thus we indicate the subface
4452  // which is located at the lower left
4453  // corner (the origin of the face's
4454  // local coordinate system) with
4455  // '0'. The rest of this issue is
4456  // taken care of using the above
4457  // conversion to a 'standard face
4458  // refine case')
4459  {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4460  {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4461  // cut_y (the same applies as for
4462  // cut_x)
4463  {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4464  {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4465  // cut_xyz: this information is
4466  // identical to the information
4467  // returned by
4468  // GeometryInfo<3>::real_to_standard_face_vertex()
4469  {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4470  // face_rotation=false, subfaces 0,1,2,3
4471  {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4472  // face_rotation=true, subfaces 0,1,2,3
4473  {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4474  // face_rotation=false, subfaces 0,1,2,3
4475  {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4476  // face_rotation=true, subfaces 0,1,2,3
4477  {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4478  // face_rotation=false, subfaces 0,1,2,3
4479  {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4480  // face_rotation=true, subfaces 0,1,2,3
4481  {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4482  // face_rotation=false, subfaces 0,1,2,3
4483  {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4484  // face_rotation=true, subfaces 0,1,2,3
4485 
4486  const unsigned int std_subface =
4487  subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4488  [subface];
4489  Assert(std_subface != e, ExcInternalError());
4490 
4491  // third step: these are the children, which
4492  // can be found at the given subfaces of an
4493  // isotropically refined (standard) face
4494  //
4495  // first index: (refinement_case-1)
4496  // second index: face_index
4497  // third index: subface_index (isotropic refinement)
4498  const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4499  [max_children_per_face] = {
4500  // cut_x
4501  {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4502  {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4503  {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4504  {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4505  {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4506  {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4507  // cut_y
4508  {{0, 1, 0, 1},
4509  {0, 1, 0, 1},
4510  {0, 0, 0, 0},
4511  {1, 1, 1, 1},
4512  {0, 0, 1, 1},
4513  {0, 0, 1, 1}},
4514  // cut_xy
4515  {{0, 2, 0, 2},
4516  {1, 3, 1, 3},
4517  {0, 0, 1, 1},
4518  {2, 2, 3, 3},
4519  {0, 1, 2, 3},
4520  {0, 1, 2, 3}},
4521  // cut_z
4522  {{0, 0, 1, 1},
4523  {0, 0, 1, 1},
4524  {0, 1, 0, 1},
4525  {0, 1, 0, 1},
4526  {0, 0, 0, 0},
4527  {1, 1, 1, 1}},
4528  // cut_xz
4529  {{0, 0, 1, 1},
4530  {2, 2, 3, 3},
4531  {0, 1, 2, 3},
4532  {0, 1, 2, 3},
4533  {0, 2, 0, 2},
4534  {1, 3, 1, 3}},
4535  // cut_yz
4536  {{0, 1, 2, 3},
4537  {0, 1, 2, 3},
4538  {0, 2, 0, 2},
4539  {1, 3, 1, 3},
4540  {0, 0, 1, 1},
4541  {2, 2, 3, 3}},
4542  // cut_xyz
4543  {{0, 2, 4, 6},
4544  {1, 3, 5, 7},
4545  {0, 4, 1, 5},
4546  {2, 6, 3, 7},
4547  {0, 1, 2, 3},
4548  {4, 5, 6, 7}}};
4549 
4550  // forth step: check, whether the given face
4551  // refine case is valid for the given cell
4552  // refine case. this is the case, if the
4553  // given face refine case is at least as
4554  // refined as the face is for the given cell
4555  // refine case
4556 
4557  // note, that we are considering standard
4558  // face refinement cases here and thus must
4559  // not pass the given orientation, flip and
4560  // rotation flags
4561  if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4562  face_refinement_case(ref_case, face))
4563  {
4564  // all is fine. for anisotropic face
4565  // refine cases, select one of the
4566  // isotropic subfaces which neighbors the
4567  // same child
4568 
4569  // first index: (standard) face refine case
4570  // second index: subface index
4571  const unsigned int equivalent_iso_subface[4][4] = {
4572  {0, e, e, e}, // no_refinement
4573  {0, 3, e, e}, // cut_x
4574  {0, 3, e, e}, // cut_y
4575  {0, 1, 2, 3}}; // cut_xy
4576 
4577  const unsigned int equ_std_subface =
4578  equivalent_iso_subface[std_face_ref][std_subface];
4579  Assert(equ_std_subface != e, ExcInternalError());
4580 
4581  return iso_children[ref_case - 1][face][equ_std_subface];
4582  }
4583  else
4584  {
4585  // the face_ref_case was too coarse,
4586  // throw an error
4587  Assert(false,
4588  ExcMessage("The face RefineCase is too coarse "
4589  "for the given cell RefineCase."));
4590  }
4591  // we only get here in case of an error
4592  return e;
4593 }
4594 
4595 
4596 
4597 template <>
4598 inline unsigned int
4600  const unsigned int,
4601  const unsigned int,
4602  const bool,
4603  const bool,
4604  const bool,
4605  const RefinementCase<3> &)
4606 {
4607  Assert(false, ExcNotImplemented());
4609 }
4610 
4611 
4612 
4613 template <>
4614 inline unsigned int
4615 GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4616  const unsigned int line,
4617  const bool,
4618  const bool,
4619  const bool)
4620 {
4621  (void)line;
4624 
4625  // The face is a line itself.
4626  return face;
4627 }
4628 
4629 
4630 
4631 template <>
4632 inline unsigned int
4633 GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4634  const unsigned int line,
4635  const bool face_orientation,
4636  const bool face_flip,
4637  const bool face_rotation)
4638 {
4641 
4642  const unsigned lines[faces_per_cell][lines_per_face] = {
4643  {8, 10, 0, 4}, // left face
4644  {9, 11, 1, 5}, // right face
4645  {2, 6, 8, 9}, // front face
4646  {3, 7, 10, 11}, // back face
4647  {0, 1, 2, 3}, // bottom face
4648  {4, 5, 6, 7}}; // top face
4649  return lines[face][real_to_standard_face_line(
4650  line, face_orientation, face_flip, face_rotation)];
4651 }
4652 
4653 
4654 
4655 inline unsigned int
4656 GeometryInfo<0>::face_to_cell_lines(const unsigned int,
4657  const unsigned int,
4658  const bool,
4659  const bool,
4660  const bool)
4661 {
4662  Assert(false, ExcNotImplemented());
4664 }
4665 
4666 
4667 
4668 template <int dim>
4669 inline unsigned int
4670 GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4671  const unsigned int,
4672  const bool,
4673  const bool,
4674  const bool)
4675 {
4676  Assert(false, ExcNotImplemented());
4678 }
4679 
4680 
4681 
4682 template <int dim>
4683 inline unsigned int
4684 GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4685  const unsigned int vertex,
4686  const bool face_orientation,
4687  const bool face_flip,
4688  const bool face_rotation)
4689 {
4691  face,
4692  vertex,
4693  face_orientation,
4694  face_flip,
4695  face_rotation);
4696 }
4697 
4698 
4699 
4700 inline unsigned int
4701 GeometryInfo<0>::face_to_cell_vertices(const unsigned int,
4702  const unsigned int,
4703  const bool,
4704  const bool,
4705  const bool)
4706 {
4707  Assert(false, ExcNotImplemented());
4709 }
4710 
4711 
4712 
4713 template <int dim>
4714 template <typename Number>
4715 inline Point<dim, Number>
4717 {
4719  for (unsigned int i = 0; i < dim; ++i)
4720  p[i] = std::min(std::max(q[i], Number(0.)), Number(1.));
4721 
4722  return p;
4723 }
4724 
4725 
4726 
4727 template <int dim>
4728 inline double
4730 {
4731  double result = 0.0;
4732 
4733  for (unsigned int i = 0; i < dim; ++i)
4734  {
4735  result = std::max(result, -p[i]);
4736  result = std::max(result, p[i] - 1.);
4737  }
4738 
4739  return result;
4740 }
4741 
4742 
4743 
4744 template <int dim>
4745 inline double
4747  const unsigned int i)
4748 {
4750 
4751  switch (dim)
4752  {
4753  case 1:
4754  {
4755  const double x = xi[0];
4756  switch (i)
4757  {
4758  case 0:
4759  return 1 - x;
4760  case 1:
4761  return x;
4762  }
4763  break;
4764  }
4765 
4766  case 2:
4767  {
4768  const double x = xi[0];
4769  const double y = xi[1];
4770  switch (i)
4771  {
4772  case 0:
4773  return (1 - x) * (1 - y);
4774  case 1:
4775  return x * (1 - y);
4776  case 2:
4777  return (1 - x) * y;
4778  case 3:
4779  return x * y;
4780  }
4781  break;
4782  }
4783 
4784  case 3:
4785  {
4786  const double x = xi[0];
4787  const double y = xi[1];
4788  const double z = xi[2];
4789  switch (i)
4790  {
4791  case 0:
4792  return (1 - x) * (1 - y) * (1 - z);
4793  case 1:
4794  return x * (1 - y) * (1 - z);
4795  case 2:
4796  return (1 - x) * y * (1 - z);
4797  case 3:
4798  return x * y * (1 - z);
4799  case 4:
4800  return (1 - x) * (1 - y) * z;
4801  case 5:
4802  return x * (1 - y) * z;
4803  case 6:
4804  return (1 - x) * y * z;
4805  case 7:
4806  return x * y * z;
4807  }
4808  break;
4809  }
4810 
4811  default:
4812  Assert(false, ExcNotImplemented());
4813  }
4814  return -1e9;
4815 }
4816 
4817 
4818 
4819 template <>
4821  const Point<1> &,
4822  const unsigned int i)
4823 {
4825 
4826  switch (i)
4827  {
4828  case 0:
4829  return Point<1>(-1.);
4830  case 1:
4831  return Point<1>(1.);
4832  }
4833 
4834  return Point<1>(-1e9);
4835 }
4836 
4837 
4838 
4839 template <>
4841  const Point<2> & xi,
4842  const unsigned int i)
4843 {
4845 
4846  const double x = xi[0];
4847  const double y = xi[1];
4848  switch (i)
4849  {
4850  case 0:
4851  return Point<2>(-(1 - y), -(1 - x));
4852  case 1:
4853  return Point<2>(1 - y, -x);
4854  case 2:
4855  return Point<2>(-y, 1 - x);
4856  case 3:
4857  return Point<2>(y, x);
4858  }
4859  return Point<2>(-1e9, -1e9);
4860 }
4861 
4862 
4863 
4864 template <>
4866  const Point<3> & xi,
4867  const unsigned int i)
4868 {
4870 
4871  const double x = xi[0];
4872  const double y = xi[1];
4873  const double z = xi[2];
4874  switch (i)
4875  {
4876  case 0:
4877  return Point<3>(-(1 - y) * (1 - z),
4878  -(1 - x) * (1 - z),
4879  -(1 - x) * (1 - y));
4880  case 1:
4881  return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4882  case 2:
4883  return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4884  case 3:
4885  return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4886  case 4:
4887  return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4888  case 5:
4889  return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4890  case 6:
4891  return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4892  case 7:
4893  return Point<3>(y * z, x * z, x * y);
4894  }
4895 
4896  return Point<3>(-1e9, -1e9, -1e9);
4897 }
4898 
4899 
4900 
4901 template <int dim>
4902 inline Tensor<1, dim>
4904  const unsigned int)
4905 {
4906  Assert(false, ExcNotImplemented());
4907  return Tensor<1, dim>();
4908 }
4909 
4910 
4911 unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4912 {
4913  return 0;
4914 }
4915 
4916 
4917 namespace internal
4918 {
4919  namespace GeometryInfoHelper
4920  {
4921  // wedge product of a single
4922  // vector in 2d: we just have to
4923  // rotate it by 90 degrees to the
4924  // right
4925  inline Tensor<1, 2>
4926  wedge_product(const Tensor<1, 2> (&derivative)[1])
4927  {
4928  Tensor<1, 2> result;
4929  result[0] = derivative[0][1];
4930  result[1] = -derivative[0][0];
4931 
4932  return result;
4933  }
4934 
4935 
4936  // wedge product of 2 vectors in
4937  // 3d is the cross product
4938  inline Tensor<1, 3>
4939  wedge_product(const Tensor<1, 3> (&derivative)[2])
4940  {
4941  return cross_product_3d(derivative[0], derivative[1]);
4942  }
4943 
4944 
4945  // wedge product of dim vectors
4946  // in dim-d: that's the
4947  // determinant of the matrix
4948  template <int dim>
4949  inline Tensor<0, dim>
4950  wedge_product(const Tensor<1, dim> (&derivative)[dim])
4951  {
4952  Tensor<2, dim> jacobian;
4953  for (unsigned int i = 0; i < dim; ++i)
4954  jacobian[i] = derivative[i];
4955 
4956  return determinant(jacobian);
4957  }
4958  } // namespace GeometryInfoHelper
4959 } // namespace internal
4960 
4961 
4962 template <int dim>
4963 template <int spacedim>
4964 inline void
4966 # ifndef DEAL_II_CXX14_CONSTEXPR_BUG
4969 # else
4971 # endif
4972 {
4973  // for each of the vertices,
4974  // compute the alternating form
4975  // of the mapped unit
4976  // vectors. consider for
4977  // example the case of a quad
4978  // in spacedim==3: to do so, we
4979  // need to see how the
4980  // infinitesimal vectors
4981  // (d\xi_1,0) and (0,d\xi_2)
4982  // are transformed into
4983  // spacedim-dimensional space
4984  // and then form their cross
4985  // product (i.e. the wedge product
4986  // of two vectors). to this end, note
4987  // that
4988  // \vec x = sum_i \vec v_i phi_i(\vec xi)
4989  // so the transformed vectors are
4990  // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
4991  // and
4992  // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
4993  // which boils down to the columns
4994  // of the 3x2 matrix \grad_\xi x(\xi)
4995  //
4996  // a similar reasoning would
4997  // hold for all dim,spacedim
4998  // pairs -- we only have to
4999  // compute the wedge product of
5000  // the columns of the
5001  // derivatives
5002  for (unsigned int i = 0; i < vertices_per_cell; ++i)
5003  {
5004  Tensor<1, spacedim> derivatives[dim];
5005 
5006  for (unsigned int j = 0; j < vertices_per_cell; ++j)
5007  {
5008  const Tensor<1, dim> grad_phi_j =
5010  for (unsigned int l = 0; l < dim; ++l)
5011  derivatives[l] += vertices[j] * grad_phi_j[l];
5012  }
5013 
5014  forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
5015  }
5016 }
5017 
5018 #endif // DOXYGEN
5020 
5021 #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:2650
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:401
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:400
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()