Reference documentation for deal.II version Git f81eda9982 2020-03-28 21:30:57 -0400
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
geometry_info.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_geometry_info_h
17 #define dealii_geometry_info_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/point.h>
24 
25 #include <boost/range/irange.hpp>
26 
27 #include <array>
28 #include <cstdint>
29 
30 
31 
32 DEAL_II_NAMESPACE_OPEN
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 std::array<std::array<Tensor<1, 1>, 0>, 2>
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 std::array<std::array<unsigned int, 1>, 2>
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 std::array<std::array<Tensor<1, 2>, 1>, 4>
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 std::array<std::array<unsigned int, 2>, 4>
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 std::array<std::array<Tensor<1, 3>, 2>, 6>
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 std::array<std::array<unsigned int, 3>, 8>
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 std::array<std::array<Tensor<1, 4>, 3>, 8>
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 std::array<std::array<unsigned int, 4>, 16>
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 
420 {
421 public:
428  enum Object
429  {
433  vertex = 0,
437  line = 1,
441  quad = 2,
445  hex = 3
446  };
447 
452  GeometryPrimitive(const Object object);
453 
459  GeometryPrimitive(const unsigned int object_dimension);
460 
465  operator unsigned int() const;
466 
467 private:
472 };
473 
474 
475 
489 template <int dim>
491 {
527  {
531  no_refinement = 0,
532 
544  isotropic_refinement = 0xFF
545  };
546 };
547 
548 
549 
560 template <>
562 {
598  {
602  no_refinement = 0,
606  cut_x = 1,
610  isotropic_refinement = cut_x
611  };
612 };
613 
614 
615 
627 template <>
629 {
665  {
669  no_refinement = 0,
673  cut_x = 1,
677  cut_y = 2,
681  cut_xy = cut_x | cut_y,
682 
686  isotropic_refinement = cut_xy
687  };
688 };
689 
690 
691 
703 template <>
705 {
741  {
745  no_refinement = 0,
749  cut_x = 1,
753  cut_y = 2,
757  cut_xy = cut_x | cut_y,
761  cut_z = 4,
765  cut_xz = cut_x | cut_z,
769  cut_yz = cut_y | cut_z,
773  cut_xyz = cut_x | cut_y | cut_z,
774 
778  isotropic_refinement = cut_xyz
779  };
780 };
781 
782 
783 
795 template <int dim>
797 {
798 public:
802  RefinementCase();
803 
809  const typename RefinementPossibilities<dim>::Possibilities refinement_case);
810 
816  explicit RefinementCase(const std::uint8_t refinement_case);
817 
829  operator std::uint8_t() const;
830 
836  operator|(const RefinementCase &r) const;
837 
842  RefinementCase operator&(const RefinementCase &r) const;
843 
852  operator~() const;
853 
854 
860  static RefinementCase
861  cut_axis(const unsigned int i);
862 
866  static std::size_t
867  memory_consumption();
868 
873  template <class Archive>
874  void
875  serialize(Archive &ar, const unsigned int version);
876 
881  ExcInvalidRefinementCase,
882  int,
883  << "The refinement flags given (" << arg1
884  << ") contain set bits that do not "
885  << "make sense for the space dimension of the object to which they are applied.");
886 
887 private:
892  std::uint8_t value : (dim > 0 ? dim : 1);
893 };
894 
895 
896 namespace internal
897 {
918  template <int dim>
920  {
925  {
929  case_none = 0,
930 
934  case_isotropic = static_cast<std::uint8_t>(-1)
935  };
936  };
937 
938 
948  template <>
950  {
957  {
961  case_none = 0,
962 
966  case_isotropic = case_none
967  };
968  };
969 
970 
971 
982  template <>
984  {
991  {
995  case_none = 0,
996 
1000  case_isotropic = case_none
1001  };
1002  };
1003 
1004 
1005 
1017  template <>
1019  {
1027  {
1031  case_none = 0,
1035  case_x = 1,
1039  case_isotropic = case_x
1040  };
1041  };
1042 
1043 
1044 
1137  template <>
1139  {
1147  {
1148  case_none = 0,
1149  case_x = 1,
1150  case_x1y = 2,
1151  case_x2y = 3,
1152  case_x1y2y = 4,
1153  case_y = 5,
1154  case_y1x = 6,
1155  case_y2x = 7,
1156  case_y1x2x = 8,
1157  case_xy = 9,
1158 
1159  case_isotropic = case_xy
1160  };
1161  };
1162 
1163 
1164 
1172  template <int dim>
1173  class SubfaceCase : public SubfacePossibilities<dim>
1174  {
1175  public:
1182  subface_possibility);
1183 
1195  operator std::uint8_t() const;
1196 
1200  static constexpr std::size_t
1201  memory_consumption();
1202 
1207  ExcInvalidSubfaceCase,
1208  int,
1209  << "The subface case given (" << arg1 << ") does not make sense "
1210  << "for the space dimension of the object to which they are applied.");
1211 
1212  private:
1217  std::uint8_t value : (dim == 3 ? 4 : 1);
1218  };
1219 
1220 } // namespace internal
1221 
1222 
1223 
1224 template <int dim>
1226 
1227 
1228 
1246 template <>
1247 struct GeometryInfo<0>
1248 {
1256  static constexpr unsigned int max_children_per_cell = 1;
1257 
1261  static constexpr unsigned int faces_per_cell = 0;
1262 
1279  static std::array<unsigned int, 0>
1280  face_indices();
1281 
1289  static constexpr unsigned int max_children_per_face = 0;
1290 
1296  static unsigned int
1297  n_children(const RefinementCase<0> &refinement_case);
1298 
1302  static constexpr unsigned int vertices_per_cell = 1;
1303 
1322  static std::array<unsigned int, vertices_per_cell>
1323  vertex_indices();
1324 
1331  static constexpr unsigned int vertices_per_face = 0;
1332 
1336  static constexpr unsigned int lines_per_face = 0;
1337 
1341  static constexpr unsigned int quads_per_face = 0;
1342 
1346  static constexpr unsigned int lines_per_cell = 0;
1347 
1351  static constexpr unsigned int quads_per_cell = 0;
1352 
1356  static constexpr unsigned int hexes_per_cell = 0;
1357 
1375  static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1376 
1390  static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1391 };
1392 
1393 
1394 
1919 template <int dim>
1920 struct GeometryInfo
1921 {
1929  static constexpr unsigned int max_children_per_cell = 1 << dim;
1930 
1934  static constexpr unsigned int faces_per_cell = 2 * dim;
1935 
1950  static boost::integer_range<unsigned int>
1951  face_indices();
1952 
1960  static constexpr unsigned int max_children_per_face =
1961  GeometryInfo<dim - 1>::max_children_per_cell;
1962 
1966  static constexpr unsigned int vertices_per_cell = 1 << dim;
1967 
1981  static boost::integer_range<unsigned int>
1982  vertex_indices();
1983 
1987  static constexpr unsigned int vertices_per_face =
1988  GeometryInfo<dim - 1>::vertices_per_cell;
1989 
1993  static constexpr unsigned int lines_per_face =
1994  GeometryInfo<dim - 1>::lines_per_cell;
1995 
1999  static constexpr unsigned int quads_per_face =
2000  GeometryInfo<dim - 1>::quads_per_cell;
2001 
2011  static constexpr unsigned int lines_per_cell =
2012  (2 * GeometryInfo<dim - 1>::lines_per_cell +
2013  GeometryInfo<dim - 1>::vertices_per_cell);
2014 
2022  static constexpr unsigned int quads_per_cell =
2023  (2 * GeometryInfo<dim - 1>::quads_per_cell +
2024  GeometryInfo<dim - 1>::lines_per_cell);
2025 
2029  static constexpr unsigned int hexes_per_cell =
2030  (2 * GeometryInfo<dim - 1>::hexes_per_cell +
2031  GeometryInfo<dim - 1>::quads_per_cell);
2032 
2050  static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
2051  internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal();
2052 
2066  static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
2067  internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal();
2068 
2079  static constexpr std::array<std::array<unsigned int, dim>, vertices_per_cell>
2080  vertex_to_face =
2081  internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face();
2082 
2087  static unsigned int
2088  n_children(const RefinementCase<dim> &refinement_case);
2089 
2094  static unsigned int
2095  n_subfaces(const internal::SubfaceCase<dim> &subface_case);
2096 
2106  static double
2107  subface_ratio(const internal::SubfaceCase<dim> &subface_case,
2108  const unsigned int subface_no);
2109 
2115  static RefinementCase<dim - 1>
2116  face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2117  const unsigned int face_no,
2118  const bool face_orientation = true,
2119  const bool face_flip = false,
2120  const bool face_rotation = false);
2121 
2127  static RefinementCase<dim>
2128  min_cell_refinement_case_for_face_refinement(
2129  const RefinementCase<dim - 1> &face_refinement_case,
2130  const unsigned int face_no,
2131  const bool face_orientation = true,
2132  const bool face_flip = false,
2133  const bool face_rotation = false);
2134 
2139  static RefinementCase<1>
2140  line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2141  const unsigned int line_no);
2142 
2147  static RefinementCase<dim>
2148  min_cell_refinement_case_for_line_refinement(const unsigned int line_no);
2149 
2196  static unsigned int
2197  child_cell_on_face(const RefinementCase<dim> & ref_case,
2198  const unsigned int face,
2199  const unsigned int subface,
2200  const bool face_orientation = true,
2201  const bool face_flip = false,
2202  const bool face_rotation = false,
2203  const RefinementCase<dim - 1> &face_refinement_case =
2205 
2219  static unsigned int
2220  line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
2221 
2242  static unsigned int
2243  face_to_cell_vertices(const unsigned int face,
2244  const unsigned int vertex,
2245  const bool face_orientation = true,
2246  const bool face_flip = false,
2247  const bool face_rotation = false);
2248 
2260  static unsigned int
2261  face_to_cell_lines(const unsigned int face,
2262  const unsigned int line,
2263  const bool face_orientation = true,
2264  const bool face_flip = false,
2265  const bool face_rotation = false);
2266 
2276  static unsigned int
2277  standard_to_real_face_vertex(const unsigned int vertex,
2278  const bool face_orientation = true,
2279  const bool face_flip = false,
2280  const bool face_rotation = false);
2281 
2291  static unsigned int
2292  real_to_standard_face_vertex(const unsigned int vertex,
2293  const bool face_orientation = true,
2294  const bool face_flip = false,
2295  const bool face_rotation = false);
2296 
2306  static unsigned int
2307  standard_to_real_face_line(const unsigned int line,
2308  const bool face_orientation = true,
2309  const bool face_flip = false,
2310  const bool face_rotation = false);
2311 
2321  static unsigned int
2322  real_to_standard_face_line(const unsigned int line,
2323  const bool face_orientation = true,
2324  const bool face_flip = false,
2325  const bool face_rotation = false);
2326 
2332  static Point<dim>
2333  unit_cell_vertex(const unsigned int vertex);
2334 
2344  static unsigned int
2345  child_cell_from_point(const Point<dim> &p);
2346 
2354  static Point<dim>
2355  cell_to_child_coordinates(const Point<dim> & p,
2356  const unsigned int child_index,
2357  const RefinementCase<dim> refine_case =
2359 
2365  static Point<dim>
2366  child_to_cell_coordinates(const Point<dim> & p,
2367  const unsigned int child_index,
2368  const RefinementCase<dim> refine_case =
2370 
2375  static bool
2376  is_inside_unit_cell(const Point<dim> &p);
2377 
2389  static bool
2390  is_inside_unit_cell(const Point<dim> &p, const double eps);
2391 
2396  static Point<dim>
2397  project_to_unit_cell(const Point<dim> &p);
2398 
2404  static double
2405  distance_to_unit_cell(const Point<dim> &p);
2406 
2411  static double
2412  d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2413 
2418  static Tensor<1, dim>
2419  d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2420 
2472  template <int spacedim>
2473  static void
2474  alternating_form_at_vertices
2475 #ifndef DEAL_II_CONSTEXPR_BUG
2476  (const Point<spacedim> (&vertices)[vertices_per_cell],
2477  Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2478 #else
2479  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2480 #endif
2481 
2491  static constexpr std::array<unsigned int, faces_per_cell>
2492  unit_normal_direction =
2493  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction();
2494 
2511  static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2512  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation();
2513 
2524  static constexpr std::array<Tensor<1, dim>, faces_per_cell>
2525  unit_normal_vector =
2526  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_vector();
2527 
2541  static constexpr std::array<std::array<Tensor<1, dim>, dim - 1>,
2542  faces_per_cell>
2543  unit_tangential_vectors = internal::GeometryInfoHelper::Initializers<
2544  dim>::unit_tangential_vectors();
2545 
2551  static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2552  internal::GeometryInfoHelper::Initializers<dim>::opposite_face();
2553 
2554 
2558  DeclException1(ExcInvalidCoordinate,
2559  double,
2560  << "The coordinates must satisfy 0 <= x_i <= 1, "
2561  << "but here we have x_i=" << arg1);
2562 
2566  DeclException3(ExcInvalidSubface,
2567  int,
2568  int,
2569  int,
2570  << "RefinementCase<dim> " << arg1 << ": face " << arg2
2571  << " has no subface " << arg3);
2572 };
2573 
2574 
2575 
2576 #ifndef DOXYGEN
2577 
2578 
2579 /* -------------- declaration of explicit specializations ------------- */
2580 
2581 
2582 template <>
2585  const unsigned int i);
2586 template <>
2589  const unsigned int i);
2590 template <>
2593  const unsigned int i);
2594 
2595 
2596 
2597 /* -------------- inline functions ------------- */
2598 
2599 
2600 inline GeometryPrimitive::GeometryPrimitive(const Object object)
2601  : object(object)
2602 {}
2603 
2604 
2605 
2606 inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2607  : object(static_cast<Object>(object_dimension))
2608 {}
2609 
2610 
2611 inline GeometryPrimitive::operator unsigned int() const
2612 {
2613  return static_cast<unsigned int>(object);
2614 }
2615 
2616 
2617 
2618 namespace internal
2619 {
2620  template <int dim>
2621  inline SubfaceCase<dim>::SubfaceCase(
2622  const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2623  : value(subface_possibility)
2624  {}
2625 
2626 
2627  template <int dim>
2628  inline SubfaceCase<dim>::operator std::uint8_t() const
2629  {
2630  return value;
2631  }
2632 
2633 
2634 } // namespace internal
2635 
2636 
2637 template <int dim>
2638 inline RefinementCase<dim>
2639 RefinementCase<dim>::cut_axis(const unsigned int)
2640 {
2641  Assert(false, ExcInternalError());
2642  return static_cast<std::uint8_t>(-1);
2643 }
2644 
2645 
2646 template <>
2647 inline RefinementCase<1>
2648 RefinementCase<1>::cut_axis(const unsigned int i)
2649 {
2650  AssertIndexRange(i, 1);
2651 
2653  return options[i];
2654 }
2655 
2656 
2657 
2658 template <>
2659 inline RefinementCase<2>
2660 RefinementCase<2>::cut_axis(const unsigned int i)
2661 {
2662  AssertIndexRange(i, 2);
2663 
2666  return options[i];
2667 }
2668 
2669 
2670 
2671 template <>
2672 inline RefinementCase<3>
2673 RefinementCase<3>::cut_axis(const unsigned int i)
2674 {
2675  AssertIndexRange(i, 3);
2676 
2680  return options[i];
2681 }
2682 
2683 
2684 
2685 template <int dim>
2687  : value(RefinementPossibilities<dim>::no_refinement)
2688 {}
2689 
2690 
2691 
2692 template <int dim>
2694  const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2695  : value(refinement_case)
2696 {
2697  // check that only those bits of
2698  // the given argument are set that
2699  // make sense for a given space
2700  // dimension
2701  Assert((refinement_case &
2703  refinement_case,
2704  ExcInvalidRefinementCase(refinement_case));
2705 }
2706 
2707 
2708 
2709 template <int dim>
2710 inline RefinementCase<dim>::RefinementCase(const std::uint8_t refinement_case)
2711  : value(refinement_case)
2712 {
2713  // check that only those bits of
2714  // the given argument are set that
2715  // make sense for a given space
2716  // dimension
2717  Assert((refinement_case &
2719  refinement_case,
2720  ExcInvalidRefinementCase(refinement_case));
2721 }
2722 
2723 
2724 
2725 template <int dim>
2726 inline RefinementCase<dim>::operator std::uint8_t() const
2727 {
2728  return value;
2729 }
2730 
2731 
2732 
2733 template <int dim>
2734 inline RefinementCase<dim>
2736 {
2737  return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2738 }
2739 
2740 
2741 
2742 template <int dim>
2744  operator&(const RefinementCase<dim> &r) const
2745 {
2746  return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2747 }
2748 
2749 
2750 
2751 template <int dim>
2752 inline RefinementCase<dim>
2754 {
2755  return RefinementCase<dim>(static_cast<std::uint8_t>(
2757 }
2758 
2759 
2760 
2761 template <int dim>
2762 inline std::size_t
2764 {
2765  return sizeof(RefinementCase<dim>);
2766 }
2767 
2768 
2769 
2770 template <int dim>
2771 template <class Archive>
2772 inline void
2773 RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2774 {
2775  // serialization can't deal with bitfields, so copy from/to a full sized
2776  // std::uint8_t
2777  std::uint8_t uchar_value = value;
2778  ar & uchar_value;
2779  value = uchar_value;
2780 }
2781 
2782 
2783 
2784 template <>
2785 inline Point<1>
2786 GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2787 {
2789 
2790  return Point<1>(static_cast<double>(vertex));
2791 }
2792 
2793 
2794 
2795 template <>
2796 inline Point<2>
2797 GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2798 {
2800 
2801  return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2802 }
2803 
2804 
2805 
2806 template <>
2807 inline Point<3>
2808 GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2809 {
2811 
2812  return {static_cast<double>(vertex % 2),
2813  static_cast<double>(vertex / 2 % 2),
2814  static_cast<double>(vertex / 4)};
2815 }
2816 
2817 
2818 
2819 inline std::array<unsigned int, 0>
2821 {
2822  return {{}};
2823 }
2824 
2825 
2826 
2827 inline std::array<unsigned int, 1>
2829 {
2830  return {{0}};
2831 }
2832 
2833 
2834 
2835 template <int dim>
2836 inline boost::integer_range<unsigned int>
2838 {
2839  return boost::irange(0U, faces_per_cell);
2840 }
2841 
2842 
2843 
2844 template <int dim>
2845 inline boost::integer_range<unsigned int>
2847 {
2848  return boost::irange(0U, vertices_per_cell);
2849 }
2850 
2851 
2852 
2853 template <int dim>
2854 inline Point<dim>
2855 GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2856 {
2857  Assert(false, ExcNotImplemented());
2858 
2859  return Point<dim>();
2860 }
2861 
2862 
2863 
2864 template <>
2865 inline unsigned int
2867 {
2868  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2869 
2870  return (p[0] <= 0.5 ? 0 : 1);
2871 }
2872 
2873 
2874 
2875 template <>
2876 inline unsigned int
2878 {
2879  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2880  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2881 
2882  return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
2883 }
2884 
2885 
2886 
2887 template <>
2888 inline unsigned int
2890 {
2891  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2892  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2893  Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
2894 
2895  return (p[0] <= 0.5 ?
2896  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
2897  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
2898 }
2899 
2900 
2901 template <int dim>
2902 inline unsigned int
2904 {
2905  Assert(false, ExcNotImplemented());
2906 
2907  return 0;
2908 }
2909 
2910 
2911 
2912 template <>
2913 inline Point<1>
2915  const unsigned int child_index,
2916  const RefinementCase<1> refine_case)
2917 
2918 {
2919  AssertIndexRange(child_index, 2);
2921  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
2922 
2923  return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
2924 }
2925 
2926 
2927 
2928 template <>
2929 inline Point<2>
2931  const unsigned int child_index,
2932  const RefinementCase<2> refine_case)
2933 
2934 {
2935  AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
2936 
2937  Point<2> point = p;
2938  switch (refine_case)
2939  {
2941  point[0] *= 2.0;
2942  if (child_index == 1)
2943  point[0] -= 1.0;
2944  break;
2946  point[1] *= 2.0;
2947  if (child_index == 1)
2948  point[1] -= 1.0;
2949  break;
2951  point *= 2.0;
2952  point -= unit_cell_vertex(child_index);
2953  break;
2954  default:
2955  Assert(false, ExcInternalError());
2956  }
2957 
2958  return point;
2959 }
2960 
2961 
2962 
2963 template <>
2964 inline Point<3>
2966  const unsigned int child_index,
2967  const RefinementCase<3> refine_case)
2968 
2969 {
2970  AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
2971 
2972  Point<3> point = p;
2973  // there might be a cleverer way to do
2974  // this, but since this function is called
2975  // in very few places for initialization
2976  // purposes only, I don't care at the
2977  // moment
2978  switch (refine_case)
2979  {
2981  point[0] *= 2.0;
2982  if (child_index == 1)
2983  point[0] -= 1.0;
2984  break;
2986  point[1] *= 2.0;
2987  if (child_index == 1)
2988  point[1] -= 1.0;
2989  break;
2991  point[2] *= 2.0;
2992  if (child_index == 1)
2993  point[2] -= 1.0;
2994  break;
2996  point[0] *= 2.0;
2997  point[1] *= 2.0;
2998  if (child_index % 2 == 1)
2999  point[0] -= 1.0;
3000  if (child_index / 2 == 1)
3001  point[1] -= 1.0;
3002  break;
3004  // careful, this is slightly
3005  // different from xy and yz due to
3006  // different internal numbering of
3007  // children!
3008  point[0] *= 2.0;
3009  point[2] *= 2.0;
3010  if (child_index / 2 == 1)
3011  point[0] -= 1.0;
3012  if (child_index % 2 == 1)
3013  point[2] -= 1.0;
3014  break;
3016  point[1] *= 2.0;
3017  point[2] *= 2.0;
3018  if (child_index % 2 == 1)
3019  point[1] -= 1.0;
3020  if (child_index / 2 == 1)
3021  point[2] -= 1.0;
3022  break;
3024  point *= 2.0;
3025  point -= unit_cell_vertex(child_index);
3026  break;
3027  default:
3028  Assert(false, ExcInternalError());
3029  }
3030 
3031  return point;
3032 }
3033 
3034 
3035 
3036 template <int dim>
3037 inline Point<dim>
3039  const Point<dim> & /*p*/,
3040  const unsigned int /*child_index*/,
3041  const RefinementCase<dim> /*refine_case*/)
3042 
3043 {
3044  Assert(false, ExcNotImplemented());
3045  return Point<dim>();
3046 }
3047 
3048 
3049 
3050 template <>
3051 inline Point<1>
3053  const unsigned int child_index,
3054  const RefinementCase<1> refine_case)
3055 
3056 {
3057  AssertIndexRange(child_index, 2);
3059  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3060 
3061  return (p + unit_cell_vertex(child_index)) * 0.5;
3062 }
3063 
3064 
3065 
3066 template <>
3067 inline Point<3>
3069  const unsigned int child_index,
3070  const RefinementCase<3> refine_case)
3071 
3072 {
3073  AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3074 
3075  Point<3> point = p;
3076  // there might be a cleverer way to do
3077  // this, but since this function is called
3078  // in very few places for initialization
3079  // purposes only, I don't care at the
3080  // moment
3081  switch (refine_case)
3082  {
3084  if (child_index == 1)
3085  point[0] += 1.0;
3086  point[0] *= 0.5;
3087  break;
3089  if (child_index == 1)
3090  point[1] += 1.0;
3091  point[1] *= 0.5;
3092  break;
3094  if (child_index == 1)
3095  point[2] += 1.0;
3096  point[2] *= 0.5;
3097  break;
3099  if (child_index % 2 == 1)
3100  point[0] += 1.0;
3101  if (child_index / 2 == 1)
3102  point[1] += 1.0;
3103  point[0] *= 0.5;
3104  point[1] *= 0.5;
3105  break;
3107  // careful, this is slightly
3108  // different from xy and yz due to
3109  // different internal numbering of
3110  // children!
3111  if (child_index / 2 == 1)
3112  point[0] += 1.0;
3113  if (child_index % 2 == 1)
3114  point[2] += 1.0;
3115  point[0] *= 0.5;
3116  point[2] *= 0.5;
3117  break;
3119  if (child_index % 2 == 1)
3120  point[1] += 1.0;
3121  if (child_index / 2 == 1)
3122  point[2] += 1.0;
3123  point[1] *= 0.5;
3124  point[2] *= 0.5;
3125  break;
3127  point += unit_cell_vertex(child_index);
3128  point *= 0.5;
3129  break;
3130  default:
3131  Assert(false, ExcInternalError());
3132  }
3133 
3134  return point;
3135 }
3136 
3137 
3138 
3139 template <>
3140 inline Point<2>
3142  const unsigned int child_index,
3143  const RefinementCase<2> refine_case)
3144 {
3145  AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3146 
3147  Point<2> point = p;
3148  switch (refine_case)
3149  {
3151  if (child_index == 1)
3152  point[0] += 1.0;
3153  point[0] *= 0.5;
3154  break;
3156  if (child_index == 1)
3157  point[1] += 1.0;
3158  point[1] *= 0.5;
3159  break;
3161  point += unit_cell_vertex(child_index);
3162  point *= 0.5;
3163  break;
3164  default:
3165  Assert(false, ExcInternalError());
3166  }
3167 
3168  return point;
3169 }
3170 
3171 
3172 
3173 template <int dim>
3174 inline Point<dim>
3176  const Point<dim> & /*p*/,
3177  const unsigned int /*child_index*/,
3178  const RefinementCase<dim> /*refine_case*/)
3179 {
3180  Assert(false, ExcNotImplemented());
3181  return Point<dim>();
3182 }
3183 
3184 
3185 
3186 template <int dim>
3187 inline bool
3189 {
3190  Assert(false, ExcNotImplemented());
3191  return false;
3192 }
3193 
3194 template <>
3195 inline bool
3197 {
3198  return (p[0] >= 0.) && (p[0] <= 1.);
3199 }
3200 
3201 
3202 
3203 template <>
3204 inline bool
3206 {
3207  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
3208 }
3209 
3210 
3211 
3212 template <>
3213 inline bool
3215 {
3216  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
3217  (p[2] >= 0.) && (p[2] <= 1.);
3218 }
3219 
3220 
3221 
3222 template <int dim>
3223 inline bool
3225 {
3226  Assert(false, ExcNotImplemented());
3227  return false;
3228 }
3229 
3230 template <>
3231 inline bool
3232 GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
3233 {
3234  return (p[0] >= -eps) && (p[0] <= 1. + eps);
3235 }
3236 
3237 
3238 
3239 template <>
3240 inline bool
3241 GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
3242 {
3243  const double l = -eps, u = 1 + eps;
3244  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
3245 }
3246 
3247 
3248 
3249 template <>
3250 inline bool
3251 GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
3252 {
3253  const double l = -eps, u = 1.0 + eps;
3254  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
3255  (p[2] >= l) && (p[2] <= u);
3256 }
3257 
3258 
3259 
3260 template <>
3261 inline unsigned int
3262 GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
3263  const unsigned int vertex)
3264 {
3265  (void)line;
3267  AssertIndexRange(vertex, 2);
3268 
3269  return vertex;
3270 }
3271 
3272 
3273 template <>
3274 inline unsigned int
3275 GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
3276  const unsigned int vertex)
3277 {
3278  constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
3279  return cell_vertices[line][vertex];
3280 }
3281 
3282 
3283 
3284 template <>
3285 inline unsigned int
3286 GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
3287  const unsigned int vertex)
3288 {
3290  AssertIndexRange(vertex, 2);
3291 
3292  constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3293  {1, 3},
3294  {0, 1},
3295  {2, 3},
3296  {4, 6}, // top face
3297  {5, 7},
3298  {4, 5},
3299  {6, 7},
3300  {0,
3301  4}, // connects of bottom
3302  {1, 5}, // top face
3303  {2, 6},
3304  {3, 7}};
3305  return vertices[line][vertex];
3306 }
3307 
3308 
3309 template <>
3310 inline unsigned int
3311 GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3312 {
3313  Assert(false, ExcNotImplemented());
3315 }
3316 
3317 template <>
3318 inline unsigned int
3319 GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3320  const bool face_orientation,
3321  const bool face_flip,
3322  const bool face_rotation)
3323 {
3325 
3326  // set up a table to make sure that
3327  // we handle non-standard faces correctly
3328  //
3329  // so set up a table that for each vertex (of
3330  // a quad in standard position) describes
3331  // which vertex to take
3332  //
3333  // first index: four vertices 0...3
3334  //
3335  // second index: face_orientation; 0:
3336  // opposite normal, 1: standard
3337  //
3338  // third index: face_flip; 0: standard, 1:
3339  // face rotated by 180 degrees
3340  //
3341  // forth index: face_rotation: 0: standard,
3342  // 1: face rotated by 90 degrees
3343 
3344  constexpr unsigned int vertex_translation[4][2][2][2] = {
3345  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3346  // face_rotation=false and true
3347  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3348  // face_rotation=false and true
3349  {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3350  // face_rotation=false and true
3351  {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3352  // face_rotation=false and true
3353 
3354  {{{2, 3}, // vertex 1 ...
3355  {1, 0}},
3356  {{1, 0}, {2, 3}}},
3357 
3358  {{{1, 0}, // vertex 2 ...
3359  {2, 3}},
3360  {{2, 3}, {1, 0}}},
3361 
3362  {{{3, 1}, // vertex 3 ...
3363  {0, 2}},
3364  {{3, 1}, {0, 2}}}};
3365 
3366  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3367 }
3368 
3369 
3370 
3371 template <int dim>
3372 inline unsigned int
3373 GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3374  const bool,
3375  const bool,
3376  const bool)
3377 {
3378  Assert(dim > 1, ExcImpossibleInDim(dim));
3380  return vertex;
3381 }
3382 
3383 template <int dim>
3384 inline unsigned int
3386 {
3387  constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3388  0, 2, 2, 4, 2, 4, 4, 8};
3389 
3390  return n_children[ref_case];
3391 }
3392 
3393 
3394 
3395 template <int dim>
3396 inline unsigned int
3398 {
3399  Assert(false, ExcNotImplemented());
3400  return 0;
3401 }
3402 
3403 template <>
3404 inline unsigned int
3406 {
3407  Assert(false, ExcImpossibleInDim(1));
3408  return 0;
3409 }
3410 
3411 template <>
3412 inline unsigned int
3414 {
3415  return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3416 }
3417 
3418 
3419 
3420 template <>
3421 inline unsigned int
3423 {
3424  const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3425  0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3426  return nsubs[subface_case];
3427 }
3428 
3429 
3430 
3431 template <int dim>
3432 inline double
3434  const unsigned int)
3435 {
3436  Assert(false, ExcNotImplemented());
3437  return 0.;
3438 }
3439 
3440 template <>
3441 inline double
3443  const unsigned int)
3444 {
3445  return 1;
3446 }
3447 
3448 
3449 template <>
3450 inline double
3452  const unsigned int)
3453 {
3454  double ratio = 1;
3455  switch (subface_case)
3456  {
3458  // Here, an
3459  // Assert(false,ExcInternalError())
3460  // would be the right
3461  // choice, but
3462  // unfortunately the
3463  // current function is
3464  // also called for faces
3465  // without children (see
3466  // tests/fe/mapping.cc).
3467  // Assert(false, ExcMessage("Face has no subfaces."));
3468  // Furthermore, assign
3469  // following value as
3470  // otherwise the
3471  // bits/volume_x tests
3472  // break
3474  break;
3476  ratio = 0.5;
3477  break;
3478  default:
3479  // there should be no
3480  // cases left
3481  Assert(false, ExcInternalError());
3482  break;
3483  }
3484 
3485  return ratio;
3486 }
3487 
3488 
3489 template <>
3490 inline double
3492  const unsigned int subface_no)
3493 {
3494  double ratio = 1;
3495  switch (subface_case)
3496  {
3498  // Here, an
3499  // Assert(false,ExcInternalError())
3500  // would be the right
3501  // choice, but
3502  // unfortunately the
3503  // current function is
3504  // also called for faces
3505  // without children (see
3506  // tests/bits/mesh_3d_16.cc). Add
3507  // following switch to
3508  // avoid diffs in
3509  // tests/bits/mesh_3d_16
3511  break;
3514  ratio = 0.5;
3515  break;
3519  ratio = 0.25;
3520  break;
3523  if (subface_no < 2)
3524  ratio = 0.25;
3525  else
3526  ratio = 0.5;
3527  break;
3530  if (subface_no == 0)
3531  ratio = 0.5;
3532  else
3533  ratio = 0.25;
3534  break;
3535  default:
3536  // there should be no
3537  // cases left
3538  Assert(false, ExcInternalError());
3539  break;
3540  }
3541 
3542  return ratio;
3543 }
3544 
3545 
3546 
3547 template <int dim>
3549  const RefinementCase<dim> &,
3550  const unsigned int,
3551  const bool,
3552  const bool,
3553  const bool)
3554 {
3555  Assert(false, ExcNotImplemented());
3556  return RefinementCase<dim - 1>::no_refinement;
3557 }
3558 
3559 template <>
3561  const RefinementCase<1> &,
3562  const unsigned int,
3563  const bool,
3564  const bool,
3565  const bool)
3566 {
3567  Assert(false, ExcImpossibleInDim(1));
3568 
3570 }
3571 
3572 
3573 template <>
3574 inline RefinementCase<1>
3576  const RefinementCase<2> &cell_refinement_case,
3577  const unsigned int face_no,
3578  const bool,
3579  const bool,
3580  const bool)
3581 {
3582  const unsigned int dim = 2;
3583  AssertIndexRange(cell_refinement_case,
3586 
3587  const RefinementCase<dim - 1>
3590  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3591  RefinementCase<dim - 1>::no_refinement},
3592 
3593  {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3594 
3595  {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3596 
3597  {RefinementCase<dim - 1>::cut_x, // cut_xy
3598  RefinementCase<dim - 1>::cut_x}};
3599 
3600  return ref_cases[cell_refinement_case][face_no / 2];
3601 }
3602 
3603 
3604 template <>
3605 inline RefinementCase<2>
3607  const RefinementCase<3> &cell_refinement_case,
3608  const unsigned int face_no,
3609  const bool face_orientation,
3610  const bool /*face_flip*/,
3611  const bool face_rotation)
3612 {
3613  const unsigned int dim = 3;
3614  AssertIndexRange(cell_refinement_case,
3617 
3618  const RefinementCase<dim - 1>
3621  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3622  RefinementCase<dim - 1>::no_refinement,
3623  RefinementCase<dim - 1>::no_refinement},
3624 
3625  {RefinementCase<dim - 1>::no_refinement, // cut_x
3626  RefinementCase<dim - 1>::cut_y,
3627  RefinementCase<dim - 1>::cut_x},
3628 
3629  {RefinementCase<dim - 1>::cut_x, // cut_y
3630  RefinementCase<dim - 1>::no_refinement,
3631  RefinementCase<dim - 1>::cut_y},
3632 
3633  {RefinementCase<dim - 1>::cut_x, // cut_xy
3634  RefinementCase<dim - 1>::cut_y,
3635  RefinementCase<dim - 1>::cut_xy},
3636 
3637  {RefinementCase<dim - 1>::cut_y, // cut_z
3638  RefinementCase<dim - 1>::cut_x,
3639  RefinementCase<dim - 1>::no_refinement},
3640 
3641  {RefinementCase<dim - 1>::cut_y, // cut_xz
3642  RefinementCase<dim - 1>::cut_xy,
3643  RefinementCase<dim - 1>::cut_x},
3644 
3645  {RefinementCase<dim - 1>::cut_xy, // cut_yz
3646  RefinementCase<dim - 1>::cut_x,
3647  RefinementCase<dim - 1>::cut_y},
3648 
3649  {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3650  RefinementCase<dim - 1>::cut_xy,
3651  RefinementCase<dim - 1>::cut_xy},
3652  };
3653 
3654  const RefinementCase<dim - 1> ref_case =
3655  ref_cases[cell_refinement_case][face_no / 2];
3656 
3657  const RefinementCase<dim - 1> flip[4] = {
3658  RefinementCase<dim - 1>::no_refinement,
3659  RefinementCase<dim - 1>::cut_y,
3660  RefinementCase<dim - 1>::cut_x,
3661  RefinementCase<dim - 1>::cut_xy};
3662 
3663  // correct the ref_case for face_orientation
3664  // and face_rotation. for face_orientation,
3665  // 'true' is the default value whereas for
3666  // face_rotation, 'false' is standard. If
3667  // <tt>face_rotation==face_orientation</tt>,
3668  // then one of them is non-standard and we
3669  // have to swap cut_x and cut_y, otherwise no
3670  // change is necessary. face_flip has no
3671  // influence. however, in order to keep the
3672  // interface consistent with other functions,
3673  // we still include it as an argument to this
3674  // function
3675  return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3676 }
3677 
3678 
3679 
3680 template <int dim>
3681 inline RefinementCase<1>
3683  const unsigned int)
3684 {
3685  Assert(false, ExcNotImplemented());
3687 }
3688 
3689 template <>
3690 inline RefinementCase<1>
3692  const RefinementCase<1> &cell_refinement_case,
3693  const unsigned int line_no)
3694 {
3695  (void)line_no;
3696  const unsigned int dim = 1;
3697  (void)dim;
3698  AssertIndexRange(cell_refinement_case,
3701 
3702  return cell_refinement_case;
3703 }
3704 
3705 
3706 template <>
3707 inline RefinementCase<1>
3709  const RefinementCase<2> &cell_refinement_case,
3710  const unsigned int line_no)
3711 {
3712  // Assertions are in face_refinement_case()
3713  return face_refinement_case(cell_refinement_case, line_no);
3714 }
3715 
3716 
3717 template <>
3718 inline RefinementCase<1>
3720  const RefinementCase<3> &cell_refinement_case,
3721  const unsigned int line_no)
3722 {
3723  const unsigned int dim = 3;
3724  AssertIndexRange(cell_refinement_case,
3727 
3728  // array indicating, which simple refine
3729  // case cuts a line in direction x, y or
3730  // z. For example, cut_y and everything
3731  // containing cut_y (cut_xy, cut_yz,
3732  // cut_xyz) cuts lines, which are in y
3733  // direction.
3734  const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3737 
3738  // order the direction of lines
3739  // 0->x, 1->y, 2->z
3740  const unsigned int direction[lines_per_cell] = {
3741  1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3742 
3743  return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3746 }
3747 
3748 
3749 
3750 template <int dim>
3751 inline RefinementCase<dim>
3753  const RefinementCase<dim - 1> &,
3754  const unsigned int,
3755  const bool,
3756  const bool,
3757  const bool)
3758 {
3759  Assert(false, ExcNotImplemented());
3760 
3762 }
3763 
3764 template <>
3765 inline RefinementCase<1>
3767  const RefinementCase<0> &,
3768  const unsigned int,
3769  const bool,
3770  const bool,
3771  const bool)
3772 {
3773  const unsigned int dim = 1;
3774  Assert(false, ExcImpossibleInDim(dim));
3775 
3777 }
3778 
3779 
3780 template <>
3781 inline RefinementCase<2>
3784  const unsigned int face_no,
3785  const bool,
3786  const bool,
3787  const bool)
3788 {
3789  const unsigned int dim = 2;
3790  AssertIndexRange(face_refinement_case,
3793 
3794  if (face_refinement_case == RefinementCase<dim>::cut_x)
3795  return (face_no / 2) ? RefinementCase<dim>::cut_x :
3797  else
3799 }
3800 
3801 
3802 template <>
3803 inline RefinementCase<3>
3805  const RefinementCase<2> &face_refinement_case,
3806  const unsigned int face_no,
3807  const bool face_orientation,
3808  const bool /*face_flip*/,
3809  const bool face_rotation)
3810 {
3811  const unsigned int dim = 3;
3812  AssertIndexRange(face_refinement_case,
3815 
3820 
3821  // correct the face_refinement_case for
3822  // face_orientation and face_rotation. for
3823  // face_orientation, 'true' is the default
3824  // value whereas for face_rotation, 'false'
3825  // is standard. If
3826  // <tt>face_rotation==face_orientation</tt>,
3827  // then one of them is non-standard and we
3828  // have to swap cut_x and cut_y, otherwise no
3829  // change is necessary. face_flip has no
3830  // influence. however, in order to keep the
3831  // interface consistent with other functions,
3832  // we still include it as an argument to this
3833  // function
3834  const RefinementCase<dim - 1> std_face_ref =
3835  (face_orientation == face_rotation) ? flip[face_refinement_case] :
3836  face_refinement_case;
3837 
3838  const RefinementCase<dim> face_to_cell[3][4] = {
3839  {RefinementCase<dim>::no_refinement, // faces 0 and 1
3840  RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3841  RefinementCase<dim>::cut_z,
3843 
3844  {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3845  // "exchanged on faces 2 and 3")
3846  RefinementCase<dim>::cut_z,
3849 
3850  {RefinementCase<dim>::no_refinement, // faces 4 and 5
3854 
3855  return face_to_cell[face_no / 2][std_face_ref];
3856 }
3857 
3858 
3859 
3860 template <int dim>
3861 inline RefinementCase<dim>
3863  const unsigned int)
3864 {
3865  Assert(false, ExcNotImplemented());
3866 
3868 }
3869 
3870 template <>
3871 inline RefinementCase<1>
3873  const unsigned int line_no)
3874 {
3875  (void)line_no;
3876  AssertIndexRange(line_no, 1);
3877 
3878  return RefinementCase<1>::cut_x;
3879 }
3880 
3881 
3882 template <>
3883 inline RefinementCase<2>
3885  const unsigned int line_no)
3886 {
3887  const unsigned int dim = 2;
3888  (void)dim;
3890 
3891  return (line_no / 2) ? RefinementCase<2>::cut_x : RefinementCase<2>::cut_y;
3892 }
3893 
3894 
3895 template <>
3896 inline RefinementCase<3>
3898  const unsigned int line_no)
3899 {
3900  const unsigned int dim = 3;
3902 
3903  const RefinementCase<dim> ref_cases[6] = {
3904  RefinementCase<dim>::cut_y, // lines 0 and 1
3905  RefinementCase<dim>::cut_x, // lines 2 and 3
3906  RefinementCase<dim>::cut_y, // lines 4 and 5
3907  RefinementCase<dim>::cut_x, // lines 6 and 7
3908  RefinementCase<dim>::cut_z, // lines 8 and 9
3909  RefinementCase<dim>::cut_z}; // lines 10 and 11
3910 
3911  return ref_cases[line_no / 2];
3912 }
3913 
3914 
3915 
3916 template <>
3917 inline unsigned int
3918 GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
3919  const bool face_orientation,
3920  const bool face_flip,
3921  const bool face_rotation)
3922 {
3924 
3925  // set up a table to make sure that
3926  // we handle non-standard faces correctly
3927  //
3928  // so set up a table that for each vertex (of
3929  // a quad in standard position) describes
3930  // which vertex to take
3931  //
3932  // first index: four vertices 0...3
3933  //
3934  // second index: face_orientation; 0:
3935  // opposite normal, 1: standard
3936  //
3937  // third index: face_flip; 0: standard, 1:
3938  // face rotated by 180 degrees
3939  //
3940  // forth index: face_rotation: 0: standard,
3941  // 1: face rotated by 90 degrees
3942 
3943  const unsigned int vertex_translation[4][2][2][2] = {
3944  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3945  // face_rotation=false and true
3946  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3947  // face_rotation=false and true
3948  {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
3949  // face_rotation=false and true
3950  {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
3951  // face_rotation=false and true
3952 
3953  {{{2, 3}, // vertex 1 ...
3954  {1, 0}},
3955  {{1, 3}, {2, 0}}},
3956 
3957  {{{1, 0}, // vertex 2 ...
3958  {2, 3}},
3959  {{2, 0}, {1, 3}}},
3960 
3961  {{{3, 1}, // vertex 3 ...
3962  {0, 2}},
3963  {{3, 2}, {0, 1}}}};
3964 
3965  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3966 }
3967 
3968 
3969 
3970 template <int dim>
3971 inline unsigned int
3972 GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
3973  const bool,
3974  const bool,
3975  const bool)
3976 {
3977  Assert(dim > 1, ExcImpossibleInDim(dim));
3979  return vertex;
3980 }
3981 
3982 
3983 
3984 template <>
3985 inline unsigned int
3986 GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
3987  const bool face_orientation,
3988  const bool face_flip,
3989  const bool face_rotation)
3990 {
3992 
3993 
3994  // make sure we handle
3995  // non-standard faces correctly
3996  //
3997  // so set up a table that for each line (of a
3998  // quad) describes which line to take
3999  //
4000  // first index: four lines 0...3
4001  //
4002  // second index: face_orientation; 0:
4003  // opposite normal, 1: standard
4004  //
4005  // third index: face_flip; 0: standard, 1:
4006  // face rotated by 180 degrees
4007  //
4008  // forth index: face_rotation: 0: standard,
4009  // 1: face rotated by 90 degrees
4010 
4011  const unsigned int line_translation[4][2][2][2] = {
4012  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4013  // face_rotation=false and true
4014  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4015  // face_rotation=false and true
4016  {{0, 3}, // line 0, face_orientation=true, face_flip=false,
4017  // face_rotation=false and true
4018  {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
4019  // face_rotation=false and true
4020 
4021  {{{3, 1}, // line 1 ...
4022  {2, 0}},
4023  {{1, 2}, {0, 3}}},
4024 
4025  {{{0, 3}, // line 2 ...
4026  {1, 2}},
4027  {{2, 0}, {3, 1}}},
4028 
4029  {{{1, 2}, // line 3 ...
4030  {0, 3}},
4031  {{3, 1}, {2, 0}}}};
4032 
4033  return line_translation[line][face_orientation][face_flip][face_rotation];
4034 }
4035 
4036 
4037 
4038 template <int dim>
4039 inline unsigned int
4040 GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
4041  const bool,
4042  const bool,
4043  const bool)
4044 {
4045  Assert(false, ExcNotImplemented());
4046  return line;
4047 }
4048 
4049 
4050 
4051 template <>
4052 inline unsigned int
4053 GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
4054  const bool face_orientation,
4055  const bool face_flip,
4056  const bool face_rotation)
4057 {
4059 
4060 
4061  // make sure we handle
4062  // non-standard faces correctly
4063  //
4064  // so set up a table that for each line (of a
4065  // quad) describes which line to take
4066  //
4067  // first index: four lines 0...3
4068  //
4069  // second index: face_orientation; 0:
4070  // opposite normal, 1: standard
4071  //
4072  // third index: face_flip; 0: standard, 1:
4073  // face rotated by 180 degrees
4074  //
4075  // forth index: face_rotation: 0: standard,
4076  // 1: face rotated by 90 degrees
4077 
4078  const unsigned int line_translation[4][2][2][2] = {
4079  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4080  // face_rotation=false and true
4081  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4082  // face_rotation=false and true
4083  {{0, 2}, // line 0, face_orientation=true, face_flip=false,
4084  // face_rotation=false and true
4085  {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
4086  // face_rotation=false and true
4087 
4088  {{{3, 1}, // line 1 ...
4089  {2, 0}},
4090  {{1, 3}, {0, 2}}},
4091 
4092  {{{0, 3}, // line 2 ...
4093  {1, 2}},
4094  {{2, 1}, {3, 0}}},
4095 
4096  {{{1, 2}, // line 3 ...
4097  {0, 3}},
4098  {{3, 0}, {2, 1}}}};
4099 
4100  return line_translation[line][face_orientation][face_flip][face_rotation];
4101 }
4102 
4103 
4104 
4105 template <int dim>
4106 inline unsigned int
4107 GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
4108  const bool,
4109  const bool,
4110  const bool)
4111 {
4112  Assert(false, ExcNotImplemented());
4113  return line;
4114 }
4115 
4116 
4117 
4118 template <>
4119 inline unsigned int
4121  const unsigned int face,
4122  const unsigned int subface,
4123  const bool,
4124  const bool,
4125  const bool,
4126  const RefinementCase<0> &)
4127 {
4128  (void)subface;
4131 
4132  return face;
4133 }
4134 
4135 
4136 
4137 template <>
4138 inline unsigned int
4140  const unsigned int face,
4141  const unsigned int subface,
4142  const bool /*face_orientation*/,
4143  const bool face_flip,
4144  const bool /*face_rotation*/,
4145  const RefinementCase<1> &)
4146 {
4149 
4150  // always return the child adjacent to the specified
4151  // subface. if the face of a cell is not refined, don't
4152  // throw an assertion but deliver the child adjacent to
4153  // the face nevertheless, i.e. deliver the child of
4154  // this cell adjacent to the subface of a possibly
4155  // refined neighbor. this simplifies setting neighbor
4156  // information in execute_refinement.
4157  const unsigned int
4159  [max_children_per_face] = {
4160  {
4161  // Normal orientation (face_flip = false)
4162  {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
4163  {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
4164  {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
4165  },
4166  {
4167  // Flipped orientation (face_flip = true)
4168  {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
4169  {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
4170  {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
4171  }};
4172 
4173  return subcells[face_flip][ref_case - 1][face][subface];
4174 }
4175 
4176 
4177 
4178 template <>
4179 inline unsigned int
4181  const unsigned int face,
4182  const unsigned int subface,
4183  const bool face_orientation,
4184  const bool face_flip,
4185  const bool face_rotation,
4186  const RefinementCase<2> &face_ref_case)
4187 {
4188  const unsigned int dim = 3;
4189 
4191  ExcMessage("Cell has no children."));
4193  if (!(subface == 0 &&
4194  face_ref_case == RefinementCase<dim - 1>::no_refinement))
4195  {
4196  AssertIndexRange(subface,
4197  GeometryInfo<dim - 1>::n_children(face_ref_case));
4198  }
4199 
4200  // invalid number used for invalid cases,
4201  // e.g. when the children are more refined at
4202  // a given face than the face itself
4203  const unsigned int e = numbers::invalid_unsigned_int;
4204 
4205  // the whole process of finding a child cell
4206  // at a given subface considering the
4207  // possibly anisotropic refinement cases of
4208  // the cell and the face as well as
4209  // orientation, flip and rotation of the face
4210  // is quite complicated. thus, we break it
4211  // down into several steps.
4212 
4213  // first step: convert the given face refine
4214  // case to a face refine case concerning the
4215  // face in standard orientation (, flip and
4216  // rotation). This only affects cut_x and
4217  // cut_y
4218  const RefinementCase<dim - 1> flip[4] = {
4219  RefinementCase<dim - 1>::no_refinement,
4220  RefinementCase<dim - 1>::cut_y,
4221  RefinementCase<dim - 1>::cut_x,
4222  RefinementCase<dim - 1>::cut_xy};
4223  // for face_orientation, 'true' is the
4224  // default value whereas for face_rotation,
4225  // 'false' is standard. If
4226  // <tt>face_rotation==face_orientation</tt>,
4227  // then one of them is non-standard and we
4228  // have to swap cut_x and cut_y, otherwise no
4229  // change is necessary.
4230  const RefinementCase<dim - 1> std_face_ref =
4231  (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4232 
4233  // second step: convert the given subface
4234  // index to the one for a standard face
4235  // respecting face_orientation, face_flip and
4236  // face_rotation
4237 
4238  // first index: face_ref_case
4239  // second index: face_orientation
4240  // third index: face_flip
4241  // forth index: face_rotation
4242  // fifth index: subface index
4243  const unsigned int subface_exchange[4][2][2][2][4] = {
4244  // no_refinement (subface 0 stays 0,
4245  // all others are invalid)
4246  {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4247  {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4248  // cut_x (here, if the face is only
4249  // rotated OR only falsely oriented,
4250  // then subface 0 of the non-standard
4251  // face does NOT correspond to one of
4252  // the subfaces of a standard
4253  // face. Thus we indicate the subface
4254  // which is located at the lower left
4255  // corner (the origin of the face's
4256  // local coordinate system) with
4257  // '0'. The rest of this issue is
4258  // taken care of using the above
4259  // conversion to a 'standard face
4260  // refine case')
4261  {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4262  {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4263  // cut_y (the same applies as for
4264  // cut_x)
4265  {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4266  {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4267  // cut_xyz: this information is
4268  // identical to the information
4269  // returned by
4270  // GeometryInfo<3>::real_to_standard_face_vertex()
4271  {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4272  // face_rotation=false, subfaces 0,1,2,3
4273  {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4274  // face_rotation=true, subfaces 0,1,2,3
4275  {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4276  // face_rotation=false, subfaces 0,1,2,3
4277  {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4278  // face_rotation=true, subfaces 0,1,2,3
4279  {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4280  // face_rotation=false, subfaces 0,1,2,3
4281  {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4282  // face_rotation=true, subfaces 0,1,2,3
4283  {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4284  // face_rotation=false, subfaces 0,1,2,3
4285  {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4286  // face_rotation=true, subfaces 0,1,2,3
4287 
4288  const unsigned int std_subface =
4289  subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4290  [subface];
4291  Assert(std_subface != e, ExcInternalError());
4292 
4293  // third step: these are the children, which
4294  // can be found at the given subfaces of an
4295  // isotropically refined (standard) face
4296  //
4297  // first index: (refinement_case-1)
4298  // second index: face_index
4299  // third index: subface_index (isotropic refinement)
4300  const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4301  [max_children_per_face] = {
4302  // cut_x
4303  {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4304  {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4305  {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4306  {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4307  {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4308  {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4309  // cut_y
4310  {{0, 1, 0, 1},
4311  {0, 1, 0, 1},
4312  {0, 0, 0, 0},
4313  {1, 1, 1, 1},
4314  {0, 0, 1, 1},
4315  {0, 0, 1, 1}},
4316  // cut_xy
4317  {{0, 2, 0, 2},
4318  {1, 3, 1, 3},
4319  {0, 0, 1, 1},
4320  {2, 2, 3, 3},
4321  {0, 1, 2, 3},
4322  {0, 1, 2, 3}},
4323  // cut_z
4324  {{0, 0, 1, 1},
4325  {0, 0, 1, 1},
4326  {0, 1, 0, 1},
4327  {0, 1, 0, 1},
4328  {0, 0, 0, 0},
4329  {1, 1, 1, 1}},
4330  // cut_xz
4331  {{0, 0, 1, 1},
4332  {2, 2, 3, 3},
4333  {0, 1, 2, 3},
4334  {0, 1, 2, 3},
4335  {0, 2, 0, 2},
4336  {1, 3, 1, 3}},
4337  // cut_yz
4338  {{0, 1, 2, 3},
4339  {0, 1, 2, 3},
4340  {0, 2, 0, 2},
4341  {1, 3, 1, 3},
4342  {0, 0, 1, 1},
4343  {2, 2, 3, 3}},
4344  // cut_xyz
4345  {{0, 2, 4, 6},
4346  {1, 3, 5, 7},
4347  {0, 4, 1, 5},
4348  {2, 6, 3, 7},
4349  {0, 1, 2, 3},
4350  {4, 5, 6, 7}}};
4351 
4352  // forth step: check, whether the given face
4353  // refine case is valid for the given cell
4354  // refine case. this is the case, if the
4355  // given face refine case is at least as
4356  // refined as the face is for the given cell
4357  // refine case
4358 
4359  // note, that we are considering standard
4360  // face refinement cases here and thus must
4361  // not pass the given orientation, flip and
4362  // rotation flags
4363  if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4364  face_refinement_case(ref_case, face))
4365  {
4366  // all is fine. for anisotropic face
4367  // refine cases, select one of the
4368  // isotropic subfaces which neighbors the
4369  // same child
4370 
4371  // first index: (standard) face refine case
4372  // second index: subface index
4373  const unsigned int equivalent_iso_subface[4][4] = {
4374  {0, e, e, e}, // no_refinement
4375  {0, 3, e, e}, // cut_x
4376  {0, 3, e, e}, // cut_y
4377  {0, 1, 2, 3}}; // cut_xy
4378 
4379  const unsigned int equ_std_subface =
4380  equivalent_iso_subface[std_face_ref][std_subface];
4381  Assert(equ_std_subface != e, ExcInternalError());
4382 
4383  return iso_children[ref_case - 1][face][equ_std_subface];
4384  }
4385  else
4386  {
4387  // the face_ref_case was too coarse,
4388  // throw an error
4389  Assert(false,
4390  ExcMessage("The face RefineCase is too coarse "
4391  "for the given cell RefineCase."));
4392  }
4393  // we only get here in case of an error
4394  return e;
4395 }
4396 
4397 
4398 
4399 template <>
4400 inline unsigned int
4402  const unsigned int,
4403  const unsigned int,
4404  const bool,
4405  const bool,
4406  const bool,
4407  const RefinementCase<3> &)
4408 {
4409  Assert(false, ExcNotImplemented());
4411 }
4412 
4413 
4414 
4415 template <>
4416 inline unsigned int
4417 GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4418  const unsigned int line,
4419  const bool,
4420  const bool,
4421  const bool)
4422 {
4423  (void)line;
4426 
4427  // The face is a line itself.
4428  return face;
4429 }
4430 
4431 
4432 
4433 template <>
4434 inline unsigned int
4435 GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4436  const unsigned int line,
4437  const bool face_orientation,
4438  const bool face_flip,
4439  const bool face_rotation)
4440 {
4443 
4444  const unsigned lines[faces_per_cell][lines_per_face] = {
4445  {8, 10, 0, 4}, // left face
4446  {9, 11, 1, 5}, // right face
4447  {2, 6, 8, 9}, // front face
4448  {3, 7, 10, 11}, // back face
4449  {0, 1, 2, 3}, // bottom face
4450  {4, 5, 6, 7}}; // top face
4451  return lines[face][real_to_standard_face_line(
4452  line, face_orientation, face_flip, face_rotation)];
4453 }
4454 
4455 
4456 
4457 template <int dim>
4458 inline unsigned int
4459 GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4460  const unsigned int,
4461  const bool,
4462  const bool,
4463  const bool)
4464 {
4465  Assert(false, ExcNotImplemented());
4467 }
4468 
4469 
4470 
4471 template <int dim>
4472 inline unsigned int
4473 GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4474  const unsigned int vertex,
4475  const bool face_orientation,
4476  const bool face_flip,
4477  const bool face_rotation)
4478 {
4480  face,
4481  vertex,
4482  face_orientation,
4483  face_flip,
4484  face_rotation);
4485 }
4486 
4487 
4488 
4489 template <int dim>
4490 inline Point<dim>
4492 {
4493  Point<dim> p = q;
4494  for (unsigned int i = 0; i < dim; i++)
4495  if (p[i] < 0.)
4496  p[i] = 0.;
4497  else if (p[i] > 1.)
4498  p[i] = 1.;
4499 
4500  return p;
4501 }
4502 
4503 
4504 
4505 template <int dim>
4506 inline double
4508 {
4509  double result = 0.0;
4510 
4511  for (unsigned int i = 0; i < dim; i++)
4512  if ((-p[i]) > result)
4513  result = -p[i];
4514  else if ((p[i] - 1.) > result)
4515  result = (p[i] - 1.);
4516 
4517  return result;
4518 }
4519 
4520 
4521 
4522 template <int dim>
4523 inline double
4525  const unsigned int i)
4526 {
4528 
4529  switch (dim)
4530  {
4531  case 1:
4532  {
4533  const double x = xi[0];
4534  switch (i)
4535  {
4536  case 0:
4537  return 1 - x;
4538  case 1:
4539  return x;
4540  }
4541  break;
4542  }
4543 
4544  case 2:
4545  {
4546  const double x = xi[0];
4547  const double y = xi[1];
4548  switch (i)
4549  {
4550  case 0:
4551  return (1 - x) * (1 - y);
4552  case 1:
4553  return x * (1 - y);
4554  case 2:
4555  return (1 - x) * y;
4556  case 3:
4557  return x * y;
4558  }
4559  break;
4560  }
4561 
4562  case 3:
4563  {
4564  const double x = xi[0];
4565  const double y = xi[1];
4566  const double z = xi[2];
4567  switch (i)
4568  {
4569  case 0:
4570  return (1 - x) * (1 - y) * (1 - z);
4571  case 1:
4572  return x * (1 - y) * (1 - z);
4573  case 2:
4574  return (1 - x) * y * (1 - z);
4575  case 3:
4576  return x * y * (1 - z);
4577  case 4:
4578  return (1 - x) * (1 - y) * z;
4579  case 5:
4580  return x * (1 - y) * z;
4581  case 6:
4582  return (1 - x) * y * z;
4583  case 7:
4584  return x * y * z;
4585  }
4586  break;
4587  }
4588 
4589  default:
4590  Assert(false, ExcNotImplemented());
4591  }
4592  return -1e9;
4593 }
4594 
4595 
4596 
4597 template <>
4599  const Point<1> &,
4600  const unsigned int i)
4601 {
4603 
4604  switch (i)
4605  {
4606  case 0:
4607  return Point<1>(-1.);
4608  case 1:
4609  return Point<1>(1.);
4610  }
4611 
4612  return Point<1>(-1e9);
4613 }
4614 
4615 
4616 
4617 template <>
4619  const Point<2> & xi,
4620  const unsigned int i)
4621 {
4623 
4624  const double x = xi[0];
4625  const double y = xi[1];
4626  switch (i)
4627  {
4628  case 0:
4629  return Point<2>(-(1 - y), -(1 - x));
4630  case 1:
4631  return Point<2>(1 - y, -x);
4632  case 2:
4633  return Point<2>(-y, 1 - x);
4634  case 3:
4635  return Point<2>(y, x);
4636  }
4637  return Point<2>(-1e9, -1e9);
4638 }
4639 
4640 
4641 
4642 template <>
4644  const Point<3> & xi,
4645  const unsigned int i)
4646 {
4648 
4649  const double x = xi[0];
4650  const double y = xi[1];
4651  const double z = xi[2];
4652  switch (i)
4653  {
4654  case 0:
4655  return Point<3>(-(1 - y) * (1 - z),
4656  -(1 - x) * (1 - z),
4657  -(1 - x) * (1 - y));
4658  case 1:
4659  return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4660  case 2:
4661  return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4662  case 3:
4663  return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4664  case 4:
4665  return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4666  case 5:
4667  return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4668  case 6:
4669  return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4670  case 7:
4671  return Point<3>(y * z, x * z, x * y);
4672  }
4673 
4674  return Point<3>(-1e9, -1e9, -1e9);
4675 }
4676 
4677 
4678 
4679 template <int dim>
4680 inline Tensor<1, dim>
4682  const unsigned int)
4683 {
4684  Assert(false, ExcNotImplemented());
4685  return Tensor<1, dim>();
4686 }
4687 
4688 
4689 unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4690 {
4691  return 0;
4692 }
4693 
4694 
4695 namespace internal
4696 {
4697  namespace GeometryInfoHelper
4698  {
4699  // wedge product of a single
4700  // vector in 2d: we just have to
4701  // rotate it by 90 degrees to the
4702  // right
4703  inline Tensor<1, 2>
4704  wedge_product(const Tensor<1, 2> (&derivative)[1])
4705  {
4706  Tensor<1, 2> result;
4707  result[0] = derivative[0][1];
4708  result[1] = -derivative[0][0];
4709 
4710  return result;
4711  }
4712 
4713 
4714  // wedge product of 2 vectors in
4715  // 3d is the cross product
4716  inline Tensor<1, 3>
4717  wedge_product(const Tensor<1, 3> (&derivative)[2])
4718  {
4719  return cross_product_3d(derivative[0], derivative[1]);
4720  }
4721 
4722 
4723  // wedge product of dim vectors
4724  // in dim-d: that's the
4725  // determinant of the matrix
4726  template <int dim>
4727  inline Tensor<0, dim>
4728  wedge_product(const Tensor<1, dim> (&derivative)[dim])
4729  {
4730  Tensor<2, dim> jacobian;
4731  for (unsigned int i = 0; i < dim; ++i)
4732  jacobian[i] = derivative[i];
4733 
4734  return determinant(jacobian);
4735  }
4736  } // namespace GeometryInfoHelper
4737 } // namespace internal
4738 
4739 
4740 template <int dim>
4741 template <int spacedim>
4742 inline void
4744 # ifndef DEAL_II_CONSTEXPR_BUG
4745  (const Point<spacedim> (&vertices)[vertices_per_cell],
4747 # else
4748  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms)
4749 # endif
4750 {
4751  // for each of the vertices,
4752  // compute the alternating form
4753  // of the mapped unit
4754  // vectors. consider for
4755  // example the case of a quad
4756  // in spacedim==3: to do so, we
4757  // need to see how the
4758  // infinitesimal vectors
4759  // (d\xi_1,0) and (0,d\xi_2)
4760  // are transformed into
4761  // spacedim-dimensional space
4762  // and then form their cross
4763  // product (i.e. the wedge product
4764  // of two vectors). to this end, note
4765  // that
4766  // \vec x = sum_i \vec v_i phi_i(\vec xi)
4767  // so the transformed vectors are
4768  // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
4769  // and
4770  // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
4771  // which boils down to the columns
4772  // of the 3x2 matrix \grad_\xi x(\xi)
4773  //
4774  // a similar reasoning would
4775  // hold for all dim,spacedim
4776  // pairs -- we only have to
4777  // compute the wedge product of
4778  // the columns of the
4779  // derivatives
4780  for (unsigned int i = 0; i < vertices_per_cell; ++i)
4781  {
4782  Tensor<1, spacedim> derivatives[dim];
4783 
4784  for (unsigned int j = 0; j < vertices_per_cell; ++j)
4785  {
4786  const Tensor<1, dim> grad_phi_j =
4788  for (unsigned int l = 0; l < dim; ++l)
4789  derivatives[l] += vertices[j] * grad_phi_j[l];
4790  }
4791 
4792  forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
4793  }
4794 }
4795 
4796 #endif // DOXYGEN
4797 DEAL_II_NAMESPACE_CLOSE
4798 
4799 #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:190
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 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)
void serialize(Archive &ar, const unsigned int version)
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
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 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
static double distance_to_unit_cell(const Point< dim > &p)
UpdateFlags operator|(const UpdateFlags f1, const UpdateFlags f2)
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:518
static constexpr unsigned int lines_per_cell
#define Assert(cond, exc)
Definition: exceptions.h:1419
static boost::integer_range< unsigned int > vertex_indices()
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
static ::ExceptionBase & ExcInvalidCoordinate(double arg1)
static std::size_t memory_consumption()
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)
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 Point< dim > project_to_unit_cell(const Point< dim > &p)
static const std::array< unsigned int, vertices_per_cell > dx_to_deal
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
Definition: tensor.h:417
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)
static unsigned int n_subfaces(const internal::SubfaceCase< dim > &subface_case)
static boost::integer_range< unsigned int > face_indices()
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 ::ExceptionBase & ExcNotImplemented()
std::uint8_t value
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:564
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
static RefinementCase cut_axis(const unsigned int i)
static ::ExceptionBase & ExcInternalError()