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