Reference documentation for deal.II version GIT 1e9e64a91d 2022-09-28 19:20:02+00:00
\(\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 - 2022 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  {
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  {
326  }
327 
328  static constexpr ::ndarray<unsigned int, 16, 4>
329  vertex_to_face()
330  {
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 {
523  enum Possibilities : std::uint8_t
524  {
529 
541  isotropic_refinement = 0xFF
542  };
543 };
544 
545 
546 
556 template <>
558 {
593  enum Possibilities : std::uint8_t
594  {
602  cut_x = 1,
606  isotropic_refinement = cut_x
607  };
608 };
609 
610 
611 
622 template <>
624 {
659  enum Possibilities : std::uint8_t
660  {
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 {
734  enum Possibilities : std::uint8_t
735  {
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:
796 
802  const typename RefinementPossibilities<dim>::Possibilities refinement_case);
803 
809  explicit RefinementCase(const std::uint8_t refinement_case);
810 
822  operator std::uint8_t() const;
823 
829  operator|(const RefinementCase &r) const;
830 
836  operator&(const RefinementCase &r) const;
837 
846  operator~() const;
847 
848 
854  static RefinementCase
855  cut_axis(const unsigned int i);
856 
860  static std::size_t
862 
868  template <class Archive>
869  void
870  serialize(Archive &ar, const unsigned int version);
871 
877  int,
878  << "The refinement flags given (" << arg1
879  << ") contain set bits that do not "
880  << "make sense for the space dimension of the object to which they are applied.");
881 
882 private:
887  std::uint8_t value : (dim > 0 ? dim : 1);
888 };
889 
890 
891 namespace internal
892 {
912  template <int dim>
914  {
919  {
924 
928  case_isotropic = static_cast<std::uint8_t>(-1)
929  };
930  };
931 
932 
941  template <>
943  {
950  {
955 
960  };
961  };
962 
963 
964 
974  template <>
976  {
983  {
988 
993  };
994  };
995 
996 
997 
1008  template <>
1010  {
1018  {
1026  case_x = 1,
1030  case_isotropic = case_x
1031  };
1032  };
1033 
1034 
1035 
1127  template <>
1129  {
1137  {
1139  case_x = 1,
1140  case_x1y = 2,
1141  case_x2y = 3,
1142  case_x1y2y = 4,
1143  case_y = 5,
1144  case_y1x = 6,
1145  case_y2x = 7,
1146  case_y1x2x = 8,
1147  case_xy = 9,
1148 
1149  case_isotropic = case_xy
1150  };
1151  };
1152 
1153 
1154 
1161  template <int dim>
1162  class SubfaceCase : public SubfacePossibilities<dim>
1163  {
1164  public:
1171  subface_possibility);
1172 
1184  operator std::uint8_t() const;
1185 
1189  static constexpr std::size_t
1191 
1197  int,
1198  << "The subface case given (" << arg1 << ") does not make sense "
1199  << "for the space dimension of the object to which they are applied.");
1200 
1201  private:
1206  std::uint8_t value : (dim == 3 ? 4 : 1);
1207  };
1208 
1209 } // namespace internal
1210 
1211 
1212 
1213 template <int dim>
1214 struct GeometryInfo;
1215 
1216 
1217 
1237 template <>
1238 struct GeometryInfo<0>
1239 {
1247  static constexpr unsigned int max_children_per_cell = 1;
1248 
1252  static constexpr unsigned int faces_per_cell = 0;
1253 
1270  static std::array<unsigned int, 0>
1272 
1280  static constexpr unsigned int max_children_per_face = 0;
1281 
1287  static unsigned int
1288  n_children(const RefinementCase<0> &refinement_case);
1289 
1293  static constexpr unsigned int vertices_per_cell = 1;
1294 
1313  static std::array<unsigned int, vertices_per_cell>
1315 
1339  static unsigned int
1340  face_to_cell_vertices(const unsigned int face,
1341  const unsigned int vertex,
1342  const bool face_orientation = true,
1343  const bool face_flip = false,
1344  const bool face_rotation = false);
1345 
1360  static unsigned int
1361  face_to_cell_lines(const unsigned int face,
1362  const unsigned int line,
1363  const bool face_orientation = true,
1364  const bool face_flip = false,
1365  const bool face_rotation = false);
1366 
1373  static constexpr unsigned int vertices_per_face = 0;
1374 
1378  static constexpr unsigned int lines_per_face = 0;
1379 
1383  static constexpr unsigned int quads_per_face = 0;
1384 
1388  static constexpr unsigned int lines_per_cell = 0;
1389 
1393  static constexpr unsigned int quads_per_cell = 0;
1394 
1398  static constexpr unsigned int hexes_per_cell = 0;
1399 
1417  static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1418 
1432  static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1433 };
1434 
1435 
1436 
1967 template <int dim>
1969 {
1977  static constexpr unsigned int max_children_per_cell = 1 << dim;
1978 
1982  static constexpr unsigned int faces_per_cell = 2 * dim;
1983 
2002 
2010  static constexpr unsigned int max_children_per_face =
2012 
2016  static constexpr unsigned int vertices_per_cell = 1 << dim;
2017 
2035 
2039  static constexpr unsigned int vertices_per_face =
2041 
2045  static constexpr unsigned int lines_per_face =
2047 
2051  static constexpr unsigned int quads_per_face =
2053 
2063  static constexpr unsigned int lines_per_cell =
2064  (2 * GeometryInfo<dim - 1>::lines_per_cell +
2066 
2074  static constexpr unsigned int quads_per_cell =
2075  (2 * GeometryInfo<dim - 1>::quads_per_cell +
2076  GeometryInfo<dim - 1>::lines_per_cell);
2077 
2081  static constexpr unsigned int hexes_per_cell =
2082  (2 * GeometryInfo<dim - 1>::hexes_per_cell +
2083  GeometryInfo<dim - 1>::quads_per_cell);
2084 
2102  static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
2103  internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal();
2104 
2118  static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
2119  internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal();
2120 
2133  internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face();
2134 
2139  static unsigned int
2140  n_children(const RefinementCase<dim> &refinement_case);
2141 
2146  static unsigned int
2148 
2158  static double
2160  const unsigned int subface_no);
2161 
2167  static RefinementCase<dim - 1>
2168  face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2169  const unsigned int face_no,
2170  const bool face_orientation = true,
2171  const bool face_flip = false,
2172  const bool face_rotation = false);
2173 
2179  static RefinementCase<dim>
2182  const unsigned int face_no,
2183  const bool face_orientation = true,
2184  const bool face_flip = false,
2185  const bool face_rotation = false);
2186 
2191  static RefinementCase<1>
2192  line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2193  const unsigned int line_no);
2194 
2199  static RefinementCase<dim>
2201 
2248  static unsigned int
2250  const unsigned int face,
2251  const unsigned int subface,
2252  const bool face_orientation = true,
2253  const bool face_flip = false,
2254  const bool face_rotation = false,
2257 
2271  static unsigned int
2272  line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
2273 
2294  static unsigned int
2295  face_to_cell_vertices(const unsigned int face,
2296  const unsigned int vertex,
2297  const bool face_orientation = true,
2298  const bool face_flip = false,
2299  const bool face_rotation = false);
2300 
2312  static unsigned int
2313  face_to_cell_lines(const unsigned int face,
2314  const unsigned int line,
2315  const bool face_orientation = true,
2316  const bool face_flip = false,
2317  const bool face_rotation = false);
2318 
2328  static unsigned int
2329  standard_to_real_face_vertex(const unsigned int vertex,
2330  const bool face_orientation = true,
2331  const bool face_flip = false,
2332  const bool face_rotation = false);
2333 
2343  static unsigned int
2344  real_to_standard_face_vertex(const unsigned int vertex,
2345  const bool face_orientation = true,
2346  const bool face_flip = false,
2347  const bool face_rotation = false);
2348 
2358  static unsigned int
2359  standard_to_real_face_line(const unsigned int line,
2360  const bool face_orientation = true,
2361  const bool face_flip = false,
2362  const bool face_rotation = false);
2363 
2369  static unsigned int
2370  standard_to_real_line_vertex(const unsigned int vertex,
2371  const bool line_orientation = true);
2372 
2380  static std::array<unsigned int, 2>
2381  standard_quad_vertex_to_line_vertex_index(const unsigned int vertex);
2382 
2390  static std::array<unsigned int, 2>
2391  standard_hex_vertex_to_quad_vertex_index(const unsigned int vertex);
2392 
2400  static std::array<unsigned int, 2>
2401  standard_hex_line_to_quad_line_index(const unsigned int line);
2402 
2412  static unsigned int
2413  real_to_standard_face_line(const unsigned int line,
2414  const bool face_orientation = true,
2415  const bool face_flip = false,
2416  const bool face_rotation = false);
2417 
2423  static Point<dim>
2424  unit_cell_vertex(const unsigned int vertex);
2425 
2435  static unsigned int
2437 
2445  static Point<dim>
2447  const unsigned int child_index,
2448  const RefinementCase<dim> refine_case =
2450 
2456  static Point<dim>
2458  const unsigned int child_index,
2459  const RefinementCase<dim> refine_case =
2461 
2466  static bool
2468 
2480  static bool
2481  is_inside_unit_cell(const Point<dim> &p, const double eps);
2482 
2487  template <typename Number = double>
2488  static Point<dim, Number>
2490 
2496  static double
2498 
2503  static double
2504  d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2505 
2510  static Tensor<1, dim>
2511  d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2512 
2564  template <int spacedim>
2565  static void
2567 #ifndef DEAL_II_CXX14_CONSTEXPR_BUG
2569  Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2570 #else
2571  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2572 #endif
2573 
2583  static constexpr std::array<unsigned int, faces_per_cell>
2585  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction();
2586 
2603  static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2604  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation();
2605 
2616  static constexpr std::array<Tensor<1, dim>, faces_per_cell>
2618  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_vector();
2619 
2633  static constexpr ndarray<Tensor<1, dim>, faces_per_cell, dim - 1>
2634  unit_tangential_vectors = internal::GeometryInfoHelper::Initializers<
2636 
2642  static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2643  internal::GeometryInfoHelper::Initializers<dim>::opposite_face();
2644 
2645 
2650  double,
2651  << "The coordinates must satisfy 0 <= x_i <= 1, "
2652  << "but here we have x_i=" << arg1);
2653 
2658  int,
2659  int,
2660  int,
2661  << "RefinementCase<dim> " << arg1 << ": face " << arg2
2662  << " has no subface " << arg3);
2663 };
2664 
2665 
2666 
2667 #ifndef DOXYGEN
2668 
2669 
2670 /* -------------- declaration of explicit specializations ------------- */
2671 
2672 
2673 template <>
2676  const unsigned int i);
2677 template <>
2680  const unsigned int i);
2681 template <>
2684  const unsigned int i);
2685 
2686 
2687 
2688 /* -------------- inline functions ------------- */
2689 
2690 
2691 inline GeometryPrimitive::GeometryPrimitive(const Object object)
2692  : object(object)
2693 {}
2694 
2695 
2696 
2697 inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2698  : object(static_cast<Object>(object_dimension))
2699 {}
2700 
2701 
2702 inline GeometryPrimitive::operator unsigned int() const
2703 {
2704  return static_cast<unsigned int>(object);
2705 }
2706 
2707 
2708 
2709 namespace internal
2710 {
2711  template <int dim>
2713  const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2714  : value(subface_possibility)
2715  {}
2716 
2717 
2718  template <int dim>
2719  inline SubfaceCase<dim>::operator std::uint8_t() const
2720  {
2721  return value;
2722  }
2723 
2724 
2725 } // namespace internal
2726 
2727 
2728 template <int dim>
2729 inline RefinementCase<dim>
2730 RefinementCase<dim>::cut_axis(const unsigned int)
2731 {
2732  Assert(false, ExcInternalError());
2733  return static_cast<std::uint8_t>(-1);
2734 }
2735 
2736 
2737 template <>
2738 inline RefinementCase<1>
2739 RefinementCase<1>::cut_axis(const unsigned int i)
2740 {
2741  AssertIndexRange(i, 1);
2742 
2744  return options[i];
2745 }
2746 
2747 
2748 
2749 template <>
2750 inline RefinementCase<2>
2751 RefinementCase<2>::cut_axis(const unsigned int i)
2752 {
2753  AssertIndexRange(i, 2);
2754 
2757  return options[i];
2758 }
2759 
2760 
2761 
2762 template <>
2763 inline RefinementCase<3>
2764 RefinementCase<3>::cut_axis(const unsigned int i)
2765 {
2766  AssertIndexRange(i, 3);
2767 
2771  return options[i];
2772 }
2773 
2774 
2775 
2776 template <int dim>
2778  : value(RefinementPossibilities<dim>::no_refinement)
2779 {}
2780 
2781 
2782 
2783 template <int dim>
2785  const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2786  : value(refinement_case)
2787 {
2788  // check that only those bits of
2789  // the given argument are set that
2790  // make sense for a given space
2791  // dimension
2792  Assert((refinement_case &
2794  refinement_case,
2795  ExcInvalidRefinementCase(refinement_case));
2796 }
2797 
2798 
2799 
2800 template <int dim>
2801 inline RefinementCase<dim>::RefinementCase(const std::uint8_t refinement_case)
2802  : value(refinement_case)
2803 {
2804  // check that only those bits of
2805  // the given argument are set that
2806  // make sense for a given space
2807  // dimension
2808  Assert((refinement_case &
2810  refinement_case,
2811  ExcInvalidRefinementCase(refinement_case));
2812 }
2813 
2814 
2815 
2816 template <int dim>
2817 inline RefinementCase<dim>::operator std::uint8_t() const
2818 {
2819  return value;
2820 }
2821 
2822 
2823 
2824 template <int dim>
2825 inline RefinementCase<dim>
2827 {
2828  return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2829 }
2830 
2831 
2832 
2833 template <int dim>
2834 inline RefinementCase<dim>
2836 {
2837  return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2838 }
2839 
2840 
2841 
2842 template <int dim>
2843 inline RefinementCase<dim>
2845 {
2846  return RefinementCase<dim>(static_cast<std::uint8_t>(
2848 }
2849 
2850 
2851 
2852 template <int dim>
2853 inline std::size_t
2855 {
2856  return sizeof(RefinementCase<dim>);
2857 }
2858 
2859 
2860 
2861 template <int dim>
2862 template <class Archive>
2863 inline void
2864 RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2865 {
2866  // serialization can't deal with bitfields, so copy from/to a full sized
2867  // std::uint8_t
2868  std::uint8_t uchar_value = value;
2869  ar & uchar_value;
2870  value = uchar_value;
2871 }
2872 
2873 
2874 
2875 template <>
2876 inline Point<1>
2877 GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2878 {
2879  AssertIndexRange(vertex, vertices_per_cell);
2880 
2881  return Point<1>(static_cast<double>(vertex));
2882 }
2883 
2884 
2885 
2886 template <>
2887 inline Point<2>
2888 GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2889 {
2890  AssertIndexRange(vertex, vertices_per_cell);
2891 
2892  return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2893 }
2894 
2895 
2896 
2897 template <>
2898 inline Point<3>
2899 GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2900 {
2901  AssertIndexRange(vertex, vertices_per_cell);
2902 
2903  return {static_cast<double>(vertex % 2),
2904  static_cast<double>(vertex / 2 % 2),
2905  static_cast<double>(vertex / 4)};
2906 }
2907 
2908 
2909 
2910 inline std::array<unsigned int, 0>
2912 {
2913  return {{}};
2914 }
2915 
2916 
2917 
2918 inline std::array<unsigned int, 1>
2920 {
2921  return {{0}};
2922 }
2923 
2924 
2925 
2926 template <int dim>
2929 {
2930  return {0U, faces_per_cell};
2931 }
2932 
2933 
2934 
2935 template <int dim>
2938 {
2939  return {0U, vertices_per_cell};
2940 }
2941 
2942 
2943 
2944 template <int dim>
2945 inline Point<dim>
2946 GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2947 {
2948  Assert(false, ExcNotImplemented());
2949 
2950  return {};
2951 }
2952 
2953 
2954 
2955 template <>
2956 inline unsigned int
2958 {
2959  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2960 
2961  return (p[0] <= 0.5 ? 0 : 1);
2962 }
2963 
2964 
2965 
2966 template <>
2967 inline unsigned int
2969 {
2970  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2971  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2972 
2973  return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
2974 }
2975 
2976 
2977 
2978 template <>
2979 inline unsigned int
2981 {
2982  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2983  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2984  Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
2985 
2986  return (p[0] <= 0.5 ?
2987  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
2988  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
2989 }
2990 
2991 
2992 template <int dim>
2993 inline unsigned int
2995 {
2996  Assert(false, ExcNotImplemented());
2997 
2998  return 0;
2999 }
3000 
3001 
3002 
3003 template <>
3004 inline Point<1>
3006  const unsigned int child_index,
3007  const RefinementCase<1> refine_case)
3008 
3009 {
3010  AssertIndexRange(child_index, 2);
3012  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3013 
3014  return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
3015 }
3016 
3017 
3018 
3019 template <>
3020 inline Point<2>
3022  const unsigned int child_index,
3023  const RefinementCase<2> refine_case)
3024 
3025 {
3026  AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3027 
3028  Point<2> point = p;
3029  switch (refine_case)
3030  {
3032  point[0] *= 2.0;
3033  if (child_index == 1)
3034  point[0] -= 1.0;
3035  break;
3037  point[1] *= 2.0;
3038  if (child_index == 1)
3039  point[1] -= 1.0;
3040  break;
3042  point *= 2.0;
3043  point -= unit_cell_vertex(child_index);
3044  break;
3045  default:
3046  Assert(false, ExcInternalError());
3047  }
3048 
3049  return point;
3050 }
3051 
3052 
3053 
3054 template <>
3055 inline Point<3>
3057  const unsigned int child_index,
3058  const RefinementCase<3> refine_case)
3059 
3060 {
3061  AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3062 
3063  Point<3> point = p;
3064  // there might be a cleverer way to do
3065  // this, but since this function is called
3066  // in very few places for initialization
3067  // purposes only, I don't care at the
3068  // moment
3069  switch (refine_case)
3070  {
3072  point[0] *= 2.0;
3073  if (child_index == 1)
3074  point[0] -= 1.0;
3075  break;
3077  point[1] *= 2.0;
3078  if (child_index == 1)
3079  point[1] -= 1.0;
3080  break;
3082  point[2] *= 2.0;
3083  if (child_index == 1)
3084  point[2] -= 1.0;
3085  break;
3087  point[0] *= 2.0;
3088  point[1] *= 2.0;
3089  if (child_index % 2 == 1)
3090  point[0] -= 1.0;
3091  if (child_index / 2 == 1)
3092  point[1] -= 1.0;
3093  break;
3095  // careful, this is slightly
3096  // different from xy and yz due to
3097  // different internal numbering of
3098  // children!
3099  point[0] *= 2.0;
3100  point[2] *= 2.0;
3101  if (child_index / 2 == 1)
3102  point[0] -= 1.0;
3103  if (child_index % 2 == 1)
3104  point[2] -= 1.0;
3105  break;
3107  point[1] *= 2.0;
3108  point[2] *= 2.0;
3109  if (child_index % 2 == 1)
3110  point[1] -= 1.0;
3111  if (child_index / 2 == 1)
3112  point[2] -= 1.0;
3113  break;
3115  point *= 2.0;
3116  point -= unit_cell_vertex(child_index);
3117  break;
3118  default:
3119  Assert(false, ExcInternalError());
3120  }
3121 
3122  return point;
3123 }
3124 
3125 
3126 
3127 template <int dim>
3128 inline Point<dim>
3130  const Point<dim> & /*p*/,
3131  const unsigned int /*child_index*/,
3132  const RefinementCase<dim> /*refine_case*/)
3133 
3134 {
3135  Assert(false, ExcNotImplemented());
3136  return {};
3137 }
3138 
3139 
3140 
3141 template <>
3142 inline Point<1>
3144  const unsigned int child_index,
3145  const RefinementCase<1> refine_case)
3146 
3147 {
3148  AssertIndexRange(child_index, 2);
3150  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3151 
3152  return (p + unit_cell_vertex(child_index)) * 0.5;
3153 }
3154 
3155 
3156 
3157 template <>
3158 inline Point<3>
3160  const unsigned int child_index,
3161  const RefinementCase<3> refine_case)
3162 
3163 {
3164  AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3165 
3166  Point<3> point = p;
3167  // there might be a cleverer way to do
3168  // this, but since this function is called
3169  // in very few places for initialization
3170  // purposes only, I don't care at the
3171  // moment
3172  switch (refine_case)
3173  {
3175  if (child_index == 1)
3176  point[0] += 1.0;
3177  point[0] *= 0.5;
3178  break;
3180  if (child_index == 1)
3181  point[1] += 1.0;
3182  point[1] *= 0.5;
3183  break;
3185  if (child_index == 1)
3186  point[2] += 1.0;
3187  point[2] *= 0.5;
3188  break;
3190  if (child_index % 2 == 1)
3191  point[0] += 1.0;
3192  if (child_index / 2 == 1)
3193  point[1] += 1.0;
3194  point[0] *= 0.5;
3195  point[1] *= 0.5;
3196  break;
3198  // careful, this is slightly
3199  // different from xy and yz due to
3200  // different internal numbering of
3201  // children!
3202  if (child_index / 2 == 1)
3203  point[0] += 1.0;
3204  if (child_index % 2 == 1)
3205  point[2] += 1.0;
3206  point[0] *= 0.5;
3207  point[2] *= 0.5;
3208  break;
3210  if (child_index % 2 == 1)
3211  point[1] += 1.0;
3212  if (child_index / 2 == 1)
3213  point[2] += 1.0;
3214  point[1] *= 0.5;
3215  point[2] *= 0.5;
3216  break;
3218  point += unit_cell_vertex(child_index);
3219  point *= 0.5;
3220  break;
3221  default:
3222  Assert(false, ExcInternalError());
3223  }
3224 
3225  return point;
3226 }
3227 
3228 
3229 
3230 template <>
3231 inline Point<2>
3233  const unsigned int child_index,
3234  const RefinementCase<2> refine_case)
3235 {
3236  AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3237 
3238  Point<2> point = p;
3239  switch (refine_case)
3240  {
3242  if (child_index == 1)
3243  point[0] += 1.0;
3244  point[0] *= 0.5;
3245  break;
3247  if (child_index == 1)
3248  point[1] += 1.0;
3249  point[1] *= 0.5;
3250  break;
3252  point += unit_cell_vertex(child_index);
3253  point *= 0.5;
3254  break;
3255  default:
3256  Assert(false, ExcInternalError());
3257  }
3258 
3259  return point;
3260 }
3261 
3262 
3263 
3264 template <int dim>
3265 inline Point<dim>
3267  const Point<dim> & /*p*/,
3268  const unsigned int /*child_index*/,
3269  const RefinementCase<dim> /*refine_case*/)
3270 {
3271  Assert(false, ExcNotImplemented());
3272  return {};
3273 }
3274 
3275 
3276 
3277 template <int dim>
3278 inline bool
3280 {
3281  Assert(false, ExcNotImplemented());
3282  return false;
3283 }
3284 
3285 template <>
3286 inline bool
3288 {
3289  return (p[0] >= 0.) && (p[0] <= 1.);
3290 }
3291 
3292 
3293 
3294 template <>
3295 inline bool
3297 {
3298  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
3299 }
3300 
3301 
3302 
3303 template <>
3304 inline bool
3306 {
3307  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
3308  (p[2] >= 0.) && (p[2] <= 1.);
3309 }
3310 
3311 
3312 
3313 template <int dim>
3314 inline bool
3316 {
3317  Assert(false, ExcNotImplemented());
3318  return false;
3319 }
3320 
3321 template <>
3322 inline bool
3323 GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
3324 {
3325  return (p[0] >= -eps) && (p[0] <= 1. + eps);
3326 }
3327 
3328 
3329 
3330 template <>
3331 inline bool
3332 GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
3333 {
3334  const double l = -eps, u = 1 + eps;
3335  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
3336 }
3337 
3338 
3339 
3340 template <>
3341 inline bool
3342 GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
3343 {
3344  const double l = -eps, u = 1.0 + eps;
3345  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
3346  (p[2] >= l) && (p[2] <= u);
3347 }
3348 
3349 
3350 
3351 template <>
3352 inline unsigned int
3353 GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
3354  const unsigned int vertex)
3355 {
3356  (void)line;
3357  AssertIndexRange(line, lines_per_cell);
3358  AssertIndexRange(vertex, 2);
3359 
3360  return vertex;
3361 }
3362 
3363 
3364 template <>
3365 inline unsigned int
3366 GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
3367  const unsigned int vertex)
3368 {
3369  constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
3370  return cell_vertices[line][vertex];
3371 }
3372 
3373 
3374 
3375 template <>
3376 inline unsigned int
3377 GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
3378  const unsigned int vertex)
3379 {
3380  AssertIndexRange(line, lines_per_cell);
3381  AssertIndexRange(vertex, 2);
3382 
3383  constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3384  {1, 3},
3385  {0, 1},
3386  {2, 3},
3387  {4, 6}, // top face
3388  {5, 7},
3389  {4, 5},
3390  {6, 7},
3391  {0,
3392  4}, // connects of bottom
3393  {1, 5}, // top face
3394  {2, 6},
3395  {3, 7}};
3396  return vertices[line][vertex];
3397 }
3398 
3399 
3400 template <>
3401 inline unsigned int
3402 GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3403 {
3404  Assert(false, ExcNotImplemented());
3406 }
3407 
3408 template <>
3409 inline unsigned int
3410 GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3411  const bool face_orientation,
3412  const bool face_flip,
3413  const bool face_rotation)
3414 {
3416 
3417  // set up a table to make sure that
3418  // we handle non-standard faces correctly
3419  //
3420  // so set up a table that for each vertex (of
3421  // a quad in standard position) describes
3422  // which vertex to take
3423  //
3424  // first index: four vertices 0...3
3425  //
3426  // second index: face_orientation; 0:
3427  // opposite normal, 1: standard
3428  //
3429  // third index: face_flip; 0: standard, 1:
3430  // face rotated by 180 degrees
3431  //
3432  // forth index: face_rotation: 0: standard,
3433  // 1: face rotated by 90 degrees
3434 
3435  constexpr unsigned int vertex_translation[4][2][2][2] = {
3436  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3437  // face_rotation=false and true
3438  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3439  // face_rotation=false and true
3440  {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3441  // face_rotation=false and true
3442  {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3443  // face_rotation=false and true
3444 
3445  {{{2, 3}, // vertex 1 ...
3446  {1, 0}},
3447  {{1, 0}, {2, 3}}},
3448 
3449  {{{1, 0}, // vertex 2 ...
3450  {2, 3}},
3451  {{2, 3}, {1, 0}}},
3452 
3453  {{{3, 1}, // vertex 3 ...
3454  {0, 2}},
3455  {{3, 1}, {0, 2}}}};
3456 
3457  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3458 }
3459 
3460 
3461 
3462 template <int dim>
3463 inline unsigned int
3464 GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3465  const bool,
3466  const bool,
3467  const bool)
3468 {
3469  Assert(dim > 1, ExcImpossibleInDim(dim));
3471  return vertex;
3472 }
3473 
3474 template <int dim>
3475 inline unsigned int
3477 {
3478  constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3479  0, 2, 2, 4, 2, 4, 4, 8};
3480 
3481  return n_children[ref_case];
3482 }
3483 
3484 
3485 
3486 template <int dim>
3487 inline unsigned int
3489 {
3490  Assert(false, ExcNotImplemented());
3491  return 0;
3492 }
3493 
3494 template <>
3495 inline unsigned int
3497 {
3498  Assert(false, ExcImpossibleInDim(1));
3499  return 0;
3500 }
3501 
3502 template <>
3503 inline unsigned int
3505 {
3506  return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3507 }
3508 
3509 
3510 
3511 template <>
3512 inline unsigned int
3514 {
3515  const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3516  0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3517  return nsubs[subface_case];
3518 }
3519 
3520 
3521 
3522 template <int dim>
3523 inline double
3525  const unsigned int)
3526 {
3527  Assert(false, ExcNotImplemented());
3528  return 0.;
3529 }
3530 
3531 template <>
3532 inline double
3534  const unsigned int)
3535 {
3536  return 1;
3537 }
3538 
3539 
3540 template <>
3541 inline double
3543  const unsigned int)
3544 {
3545  double ratio = 1;
3546  switch (subface_case)
3547  {
3549  // Here, an
3550  // Assert(false,ExcInternalError())
3551  // would be the right
3552  // choice, but
3553  // unfortunately the
3554  // current function is
3555  // also called for faces
3556  // without children (see
3557  // tests/fe/mapping.cc).
3558  // Assert(false, ExcMessage("Face has no subfaces."));
3559  // Furthermore, assign
3560  // following value as
3561  // otherwise the
3562  // bits/volume_x tests
3563  // break
3565  break;
3567  ratio = 0.5;
3568  break;
3569  default:
3570  // there should be no
3571  // cases left
3572  Assert(false, ExcInternalError());
3573  break;
3574  }
3575 
3576  return ratio;
3577 }
3578 
3579 
3580 template <>
3581 inline double
3583  const unsigned int subface_no)
3584 {
3585  double ratio = 1;
3586  switch (subface_case)
3587  {
3589  // Here, an
3590  // Assert(false,ExcInternalError())
3591  // would be the right
3592  // choice, but
3593  // unfortunately the
3594  // current function is
3595  // also called for faces
3596  // without children (see
3597  // tests/bits/mesh_3d_16.cc). Add
3598  // following switch to
3599  // avoid diffs in
3600  // tests/bits/mesh_3d_16
3602  break;
3605  ratio = 0.5;
3606  break;
3610  ratio = 0.25;
3611  break;
3614  if (subface_no < 2)
3615  ratio = 0.25;
3616  else
3617  ratio = 0.5;
3618  break;
3621  if (subface_no == 0)
3622  ratio = 0.5;
3623  else
3624  ratio = 0.25;
3625  break;
3626  default:
3627  // there should be no
3628  // cases left
3629  Assert(false, ExcInternalError());
3630  break;
3631  }
3632 
3633  return ratio;
3634 }
3635 
3636 
3637 
3638 template <int dim>
3640  const RefinementCase<dim> &,
3641  const unsigned int,
3642  const bool,
3643  const bool,
3644  const bool)
3645 {
3646  Assert(false, ExcNotImplemented());
3647  return RefinementCase<dim - 1>::no_refinement;
3648 }
3649 
3650 template <>
3652  const RefinementCase<1> &,
3653  const unsigned int,
3654  const bool,
3655  const bool,
3656  const bool)
3657 {
3658  Assert(false, ExcImpossibleInDim(1));
3659 
3661 }
3662 
3663 
3664 template <>
3665 inline RefinementCase<1>
3667  const RefinementCase<2> &cell_refinement_case,
3668  const unsigned int face_no,
3669  const bool,
3670  const bool,
3671  const bool)
3672 {
3673  const unsigned int dim = 2;
3674  AssertIndexRange(cell_refinement_case,
3677 
3678  // simple special case
3679  if (cell_refinement_case == RefinementCase<dim>::cut_xy)
3680  return RefinementCase<1>::cut_x;
3681 
3682  const RefinementCase<dim - 1>
3685  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3686  RefinementCase<dim - 1>::no_refinement},
3687 
3688  {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3689 
3690  {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3691 
3692  {RefinementCase<dim - 1>::cut_x, // cut_xy
3693  RefinementCase<dim - 1>::cut_x}};
3694 
3695  return ref_cases[cell_refinement_case][face_no / 2];
3696 }
3697 
3698 
3699 template <>
3700 inline RefinementCase<2>
3702  const RefinementCase<3> &cell_refinement_case,
3703  const unsigned int face_no,
3704  const bool face_orientation,
3705  const bool /*face_flip*/,
3706  const bool face_rotation)
3707 {
3708  const unsigned int dim = 3;
3709  AssertIndexRange(cell_refinement_case,
3712 
3713  // simple special case
3714  if (cell_refinement_case == RefinementCase<dim>::cut_xyz)
3715  return RefinementCase<dim - 1>::cut_xy;
3716 
3717  const RefinementCase<dim - 1>
3720  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3721  RefinementCase<dim - 1>::no_refinement,
3722  RefinementCase<dim - 1>::no_refinement},
3723 
3724  {RefinementCase<dim - 1>::no_refinement, // cut_x
3725  RefinementCase<dim - 1>::cut_y,
3726  RefinementCase<dim - 1>::cut_x},
3727 
3728  {RefinementCase<dim - 1>::cut_x, // cut_y
3729  RefinementCase<dim - 1>::no_refinement,
3730  RefinementCase<dim - 1>::cut_y},
3731 
3732  {RefinementCase<dim - 1>::cut_x, // cut_xy
3733  RefinementCase<dim - 1>::cut_y,
3734  RefinementCase<dim - 1>::cut_xy},
3735 
3736  {RefinementCase<dim - 1>::cut_y, // cut_z
3737  RefinementCase<dim - 1>::cut_x,
3738  RefinementCase<dim - 1>::no_refinement},
3739 
3740  {RefinementCase<dim - 1>::cut_y, // cut_xz
3741  RefinementCase<dim - 1>::cut_xy,
3742  RefinementCase<dim - 1>::cut_x},
3743 
3744  {RefinementCase<dim - 1>::cut_xy, // cut_yz
3745  RefinementCase<dim - 1>::cut_x,
3746  RefinementCase<dim - 1>::cut_y},
3747 
3748  {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3749  RefinementCase<dim - 1>::cut_xy,
3750  RefinementCase<dim - 1>::cut_xy},
3751  };
3752 
3753  const RefinementCase<dim - 1> ref_case =
3754  ref_cases[cell_refinement_case][face_no / 2];
3755 
3756  const RefinementCase<dim - 1> flip[4] = {
3757  RefinementCase<dim - 1>::no_refinement,
3758  RefinementCase<dim - 1>::cut_y,
3759  RefinementCase<dim - 1>::cut_x,
3760  RefinementCase<dim - 1>::cut_xy};
3761 
3762  // correct the ref_case for face_orientation
3763  // and face_rotation. for face_orientation,
3764  // 'true' is the default value whereas for
3765  // face_rotation, 'false' is standard. If
3766  // <tt>face_rotation==face_orientation</tt>,
3767  // then one of them is non-standard and we
3768  // have to swap cut_x and cut_y, otherwise no
3769  // change is necessary. face_flip has no
3770  // influence. however, in order to keep the
3771  // interface consistent with other functions,
3772  // we still include it as an argument to this
3773  // function
3774  return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3775 }
3776 
3777 
3778 
3779 template <int dim>
3780 inline RefinementCase<1>
3782  const unsigned int)
3783 {
3784  Assert(false, ExcNotImplemented());
3786 }
3787 
3788 template <>
3789 inline RefinementCase<1>
3791  const RefinementCase<1> &cell_refinement_case,
3792  const unsigned int line_no)
3793 {
3794  (void)line_no;
3795  const unsigned int dim = 1;
3796  (void)dim;
3797  AssertIndexRange(cell_refinement_case,
3800 
3801  return cell_refinement_case;
3802 }
3803 
3804 
3805 template <>
3806 inline RefinementCase<1>
3808  const RefinementCase<2> &cell_refinement_case,
3809  const unsigned int line_no)
3810 {
3811  // Assertions are in face_refinement_case()
3812  return face_refinement_case(cell_refinement_case, line_no);
3813 }
3814 
3815 
3816 template <>
3817 inline RefinementCase<1>
3819  const RefinementCase<3> &cell_refinement_case,
3820  const unsigned int line_no)
3821 {
3822  const unsigned int dim = 3;
3823  AssertIndexRange(cell_refinement_case,
3826 
3827  // simple special case
3828  if (cell_refinement_case == RefinementCase<dim>::cut_xyz)
3829  return RefinementCase<1>::cut_x;
3830 
3831  // array indicating, which simple refine
3832  // case cuts a line in direction x, y or
3833  // z. For example, cut_y and everything
3834  // containing cut_y (cut_xy, cut_yz,
3835  // cut_xyz) cuts lines, which are in y
3836  // direction.
3837  const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3840 
3841  // order the direction of lines
3842  // 0->x, 1->y, 2->z
3843  const unsigned int direction[lines_per_cell] = {
3844  1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3845 
3846  return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3849 }
3850 
3851 
3852 
3853 template <int dim>
3854 inline RefinementCase<dim>
3856  const RefinementCase<dim - 1> &,
3857  const unsigned int,
3858  const bool,
3859  const bool,
3860  const bool)
3861 {
3862  Assert(false, ExcNotImplemented());
3863 
3865 }
3866 
3867 template <>
3868 inline RefinementCase<1>
3870  const RefinementCase<0> &,
3871  const unsigned int,
3872  const bool,
3873  const bool,
3874  const bool)
3875 {
3876  const unsigned int dim = 1;
3877  Assert(false, ExcImpossibleInDim(dim));
3878 
3880 }
3881 
3882 
3883 template <>
3884 inline RefinementCase<2>
3886  const RefinementCase<1> &face_refinement_case,
3887  const unsigned int face_no,
3888  const bool,
3889  const bool,
3890  const bool)
3891 {
3892  const unsigned int dim = 2;
3893  AssertIndexRange(face_refinement_case,
3896 
3897  if (face_refinement_case == RefinementCase<dim>::cut_x)
3898  return (face_no / 2) != 0u ? RefinementCase<dim>::cut_x :
3900  else
3902 }
3903 
3904 
3905 template <>
3906 inline RefinementCase<3>
3908  const RefinementCase<2> &face_refinement_case,
3909  const unsigned int face_no,
3910  const bool face_orientation,
3911  const bool /*face_flip*/,
3912  const bool face_rotation)
3913 {
3914  const unsigned int dim = 3;
3915  AssertIndexRange(face_refinement_case,
3918 
3923 
3924  // correct the face_refinement_case for
3925  // face_orientation and face_rotation. for
3926  // face_orientation, 'true' is the default
3927  // value whereas for face_rotation, 'false'
3928  // is standard. If
3929  // <tt>face_rotation==face_orientation</tt>,
3930  // then one of them is non-standard and we
3931  // have to swap cut_x and cut_y, otherwise no
3932  // change is necessary. face_flip has no
3933  // influence. however, in order to keep the
3934  // interface consistent with other functions,
3935  // we still include it as an argument to this
3936  // function
3937  const RefinementCase<dim - 1> std_face_ref =
3938  (face_orientation == face_rotation) ? flip[face_refinement_case] :
3939  face_refinement_case;
3940 
3941  const RefinementCase<dim> face_to_cell[3][4] = {
3942  {RefinementCase<dim>::no_refinement, // faces 0 and 1
3943  RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3946 
3947  {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3948  // "exchanged on faces 2 and 3")
3952 
3953  {RefinementCase<dim>::no_refinement, // faces 4 and 5
3957 
3958  return face_to_cell[face_no / 2][std_face_ref];
3959 }
3960 
3961 
3962 
3963 template <int dim>
3964 inline RefinementCase<dim>
3966  const unsigned int)
3967 {
3968  Assert(false, ExcNotImplemented());
3969 
3971 }
3972 
3973 template <>
3974 inline RefinementCase<1>
3976  const unsigned int line_no)
3977 {
3978  (void)line_no;
3979  AssertIndexRange(line_no, 1);
3980 
3981  return RefinementCase<1>::cut_x;
3982 }
3983 
3984 
3985 template <>
3986 inline RefinementCase<2>
3988  const unsigned int line_no)
3989 {
3990  const unsigned int dim = 2;
3991  (void)dim;
3993 
3994  return (line_no / 2) != 0u ? RefinementCase<2>::cut_x :
3996 }
3997 
3998 
3999 template <>
4000 inline RefinementCase<3>
4002  const unsigned int line_no)
4003 {
4004  const unsigned int dim = 3;
4006 
4007  const RefinementCase<dim> ref_cases[6] = {
4008  RefinementCase<dim>::cut_y, // lines 0 and 1
4009  RefinementCase<dim>::cut_x, // lines 2 and 3
4010  RefinementCase<dim>::cut_y, // lines 4 and 5
4011  RefinementCase<dim>::cut_x, // lines 6 and 7
4012  RefinementCase<dim>::cut_z, // lines 8 and 9
4013  RefinementCase<dim>::cut_z}; // lines 10 and 11
4014 
4015  return ref_cases[line_no / 2];
4016 }
4017 
4018 
4019 
4020 template <>
4021 inline unsigned int
4022 GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
4023  const bool face_orientation,
4024  const bool face_flip,
4025  const bool face_rotation)
4026 {
4028 
4029  // set up a table to make sure that
4030  // we handle non-standard faces correctly
4031  //
4032  // so set up a table that for each vertex (of
4033  // a quad in standard position) describes
4034  // which vertex to take
4035  //
4036  // first index: four vertices 0...3
4037  //
4038  // second index: face_orientation; 0:
4039  // opposite normal, 1: standard
4040  //
4041  // third index: face_flip; 0: standard, 1:
4042  // face rotated by 180 degrees
4043  //
4044  // forth index: face_rotation: 0: standard,
4045  // 1: face rotated by 90 degrees
4046 
4047  const unsigned int vertex_translation[4][2][2][2] = {
4048  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
4049  // face_rotation=false and true
4050  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
4051  // face_rotation=false and true
4052  {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
4053  // face_rotation=false and true
4054  {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
4055  // face_rotation=false and true
4056 
4057  {{{2, 3}, // vertex 1 ...
4058  {1, 0}},
4059  {{1, 3}, {2, 0}}},
4060 
4061  {{{1, 0}, // vertex 2 ...
4062  {2, 3}},
4063  {{2, 0}, {1, 3}}},
4064 
4065  {{{3, 1}, // vertex 3 ...
4066  {0, 2}},
4067  {{3, 2}, {0, 1}}}};
4068 
4069  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
4070 }
4071 
4072 
4073 
4074 template <int dim>
4075 inline unsigned int
4076 GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
4077  const bool,
4078  const bool,
4079  const bool)
4080 {
4081  Assert(dim > 1, ExcImpossibleInDim(dim));
4083  return vertex;
4084 }
4085 
4086 
4087 
4088 template <>
4089 inline unsigned int
4090 GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
4091  const bool face_orientation,
4092  const bool face_flip,
4093  const bool face_rotation)
4094 {
4096 
4097 
4098  // make sure we handle
4099  // non-standard faces correctly
4100  //
4101  // so set up a table that for each line (of a
4102  // quad) describes which line to take
4103  //
4104  // first index: four lines 0...3
4105  //
4106  // second index: face_orientation; 0:
4107  // opposite normal, 1: standard
4108  //
4109  // third index: face_flip; 0: standard, 1:
4110  // face rotated by 180 degrees
4111  //
4112  // forth index: face_rotation: 0: standard,
4113  // 1: face rotated by 90 degrees
4114 
4115  const unsigned int line_translation[4][2][2][2] = {
4116  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4117  // face_rotation=false and true
4118  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4119  // face_rotation=false and true
4120  {{0, 3}, // line 0, face_orientation=true, face_flip=false,
4121  // face_rotation=false and true
4122  {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
4123  // face_rotation=false and true
4124 
4125  {{{3, 1}, // line 1 ...
4126  {2, 0}},
4127  {{1, 2}, {0, 3}}},
4128 
4129  {{{0, 3}, // line 2 ...
4130  {1, 2}},
4131  {{2, 0}, {3, 1}}},
4132 
4133  {{{1, 2}, // line 3 ...
4134  {0, 3}},
4135  {{3, 1}, {2, 0}}}};
4136 
4137  return line_translation[line][face_orientation][face_flip][face_rotation];
4138 }
4139 
4140 
4141 
4142 template <int dim>
4143 inline unsigned int
4144 GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
4145  const bool,
4146  const bool,
4147  const bool)
4148 {
4149  Assert(false, ExcNotImplemented());
4150  return line;
4151 }
4152 
4153 
4154 
4155 template <>
4156 inline unsigned int
4157 GeometryInfo<2>::standard_to_real_line_vertex(const unsigned int vertex,
4158  const bool line_orientation)
4159 {
4160  return line_orientation ? vertex : (1 - vertex);
4161 }
4162 
4163 
4164 
4165 template <int dim>
4166 inline unsigned int
4167 GeometryInfo<dim>::standard_to_real_line_vertex(const unsigned int vertex,
4168  const bool)
4169 {
4170  Assert(false, ExcNotImplemented());
4171  return vertex;
4172 }
4173 
4174 
4175 
4176 template <>
4177 inline std::array<unsigned int, 2>
4179  const unsigned int vertex)
4180 {
4181  return {{vertex % 2, vertex / 2}};
4182 }
4183 
4184 
4185 
4186 template <int dim>
4187 inline std::array<unsigned int, 2>
4189  const unsigned int vertex)
4190 {
4191  Assert(false, ExcNotImplemented());
4192  (void)vertex;
4193  return {{0, 0}};
4194 }
4195 
4196 
4197 
4198 template <>
4199 inline std::array<unsigned int, 2>
4201 {
4202  // set up a table that for each
4203  // line describes a) from which
4204  // quad to take it, b) which line
4205  // therein it is if the face is
4206  // oriented correctly
4207  static const unsigned int lookup_table[GeometryInfo<3>::lines_per_cell][2] = {
4208  {4, 0}, // take first four lines from bottom face
4209  {4, 1},
4210  {4, 2},
4211  {4, 3},
4212 
4213  {5, 0}, // second four lines from top face
4214  {5, 1},
4215  {5, 2},
4216  {5, 3},
4217 
4218  {0, 0}, // the rest randomly
4219  {1, 0},
4220  {0, 1},
4221  {1, 1}};
4222 
4223  return {{lookup_table[i][0], lookup_table[i][1]}};
4224 }
4225 
4226 
4227 
4228 template <int dim>
4229 inline std::array<unsigned int, 2>
4231 {
4232  Assert(false, ExcNotImplemented());
4233  (void)line;
4234  return {{0, 0}};
4235 }
4236 
4237 
4238 
4239 template <>
4240 inline std::array<unsigned int, 2>
4242  const unsigned int vertex)
4243 {
4244  // get the corner indices by asking either the bottom or the top face for its
4245  // vertices. handle non-standard faces by calling the vertex reordering
4246  // function from GeometryInfo
4247 
4248  // bottom face (4) for first four vertices, top face (5) for the rest
4249  return {{4 + vertex / 4, vertex % 4}};
4250 }
4251 
4252 
4253 
4254 template <int dim>
4255 inline std::array<unsigned int, 2>
4257  const unsigned int vertex)
4258 {
4259  Assert(false, ExcNotImplemented());
4260  (void)vertex;
4261  return {{0, 0}};
4262 }
4263 
4264 
4265 
4266 template <>
4267 inline unsigned int
4268 GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
4269  const bool face_orientation,
4270  const bool face_flip,
4271  const bool face_rotation)
4272 {
4274 
4275 
4276  // make sure we handle
4277  // non-standard faces correctly
4278  //
4279  // so set up a table that for each line (of a
4280  // quad) describes which line to take
4281  //
4282  // first index: four lines 0...3
4283  //
4284  // second index: face_orientation; 0:
4285  // opposite normal, 1: standard
4286  //
4287  // third index: face_flip; 0: standard, 1:
4288  // face rotated by 180 degrees
4289  //
4290  // forth index: face_rotation: 0: standard,
4291  // 1: face rotated by 90 degrees
4292 
4293  const unsigned int line_translation[4][2][2][2] = {
4294  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4295  // face_rotation=false and true
4296  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4297  // face_rotation=false and true
4298  {{0, 2}, // line 0, face_orientation=true, face_flip=false,
4299  // face_rotation=false and true
4300  {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
4301  // face_rotation=false and true
4302 
4303  {{{3, 1}, // line 1 ...
4304  {2, 0}},
4305  {{1, 3}, {0, 2}}},
4306 
4307  {{{0, 3}, // line 2 ...
4308  {1, 2}},
4309  {{2, 1}, {3, 0}}},
4310 
4311  {{{1, 2}, // line 3 ...
4312  {0, 3}},
4313  {{3, 0}, {2, 1}}}};
4314 
4315  return line_translation[line][face_orientation][face_flip][face_rotation];
4316 }
4317 
4318 
4319 
4320 template <int dim>
4321 inline unsigned int
4322 GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
4323  const bool,
4324  const bool,
4325  const bool)
4326 {
4327  Assert(false, ExcNotImplemented());
4328  return line;
4329 }
4330 
4331 
4332 
4333 template <>
4334 inline unsigned int
4336  const unsigned int face,
4337  const unsigned int subface,
4338  const bool,
4339  const bool,
4340  const bool,
4341  const RefinementCase<0> &)
4342 {
4343  (void)subface;
4344  AssertIndexRange(face, faces_per_cell);
4345  AssertIndexRange(subface, max_children_per_face);
4346 
4347  return face;
4348 }
4349 
4350 
4351 
4352 template <>
4353 inline unsigned int
4355  const unsigned int face,
4356  const unsigned int subface,
4357  const bool /*face_orientation*/,
4358  const bool face_flip,
4359  const bool /*face_rotation*/,
4360  const RefinementCase<1> &)
4361 {
4362  AssertIndexRange(face, faces_per_cell);
4363  AssertIndexRange(subface, max_children_per_face);
4364 
4365  // always return the child adjacent to the specified
4366  // subface. if the face of a cell is not refined, don't
4367  // throw an assertion but deliver the child adjacent to
4368  // the face nevertheless, i.e. deliver the child of
4369  // this cell adjacent to the subface of a possibly
4370  // refined neighbor. this simplifies setting neighbor
4371  // information in execute_refinement.
4372  const unsigned int
4373  subcells[2][RefinementCase<2>::isotropic_refinement][faces_per_cell]
4374  [max_children_per_face] = {
4375  {
4376  // Normal orientation (face_flip = false)
4377  {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
4378  {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
4379  {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
4380  },
4381  {
4382  // Flipped orientation (face_flip = true)
4383  {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
4384  {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
4385  {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
4386  }};
4387 
4388  return subcells[face_flip][ref_case - 1][face][subface];
4389 }
4390 
4391 
4392 
4393 template <>
4394 inline unsigned int
4396  const unsigned int face,
4397  const unsigned int subface,
4398  const bool face_orientation,
4399  const bool face_flip,
4400  const bool face_rotation,
4401  const RefinementCase<2> &face_ref_case)
4402 {
4403  const unsigned int dim = 3;
4404 
4406  ExcMessage("Cell has no children."));
4407  AssertIndexRange(face, faces_per_cell);
4408  if (!(subface == 0 &&
4409  face_ref_case == RefinementCase<dim - 1>::no_refinement))
4410  {
4411  AssertIndexRange(subface,
4412  GeometryInfo<dim - 1>::n_children(face_ref_case));
4413  }
4414 
4415  // invalid number used for invalid cases,
4416  // e.g. when the children are more refined at
4417  // a given face than the face itself
4418  const unsigned int e = numbers::invalid_unsigned_int;
4419 
4420  // the whole process of finding a child cell
4421  // at a given subface considering the
4422  // possibly anisotropic refinement cases of
4423  // the cell and the face as well as
4424  // orientation, flip and rotation of the face
4425  // is quite complicated. thus, we break it
4426  // down into several steps.
4427 
4428  // first step: convert the given face refine
4429  // case to a face refine case concerning the
4430  // face in standard orientation (, flip and
4431  // rotation). This only affects cut_x and
4432  // cut_y
4433  const RefinementCase<dim - 1> flip[4] = {
4434  RefinementCase<dim - 1>::no_refinement,
4435  RefinementCase<dim - 1>::cut_y,
4436  RefinementCase<dim - 1>::cut_x,
4437  RefinementCase<dim - 1>::cut_xy};
4438  // for face_orientation, 'true' is the
4439  // default value whereas for face_rotation,
4440  // 'false' is standard. If
4441  // <tt>face_rotation==face_orientation</tt>,
4442  // then one of them is non-standard and we
4443  // have to swap cut_x and cut_y, otherwise no
4444  // change is necessary.
4445  const RefinementCase<dim - 1> std_face_ref =
4446  (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4447 
4448  // second step: convert the given subface
4449  // index to the one for a standard face
4450  // respecting face_orientation, face_flip and
4451  // face_rotation
4452 
4453  // first index: face_ref_case
4454  // second index: face_orientation
4455  // third index: face_flip
4456  // forth index: face_rotation
4457  // fifth index: subface index
4458  const unsigned int subface_exchange[4][2][2][2][4] = {
4459  // no_refinement (subface 0 stays 0,
4460  // all others are invalid)
4461  {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4462  {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4463  // cut_x (here, if the face is only
4464  // rotated OR only falsely oriented,
4465  // then subface 0 of the non-standard
4466  // face does NOT correspond to one of
4467  // the subfaces of a standard
4468  // face. Thus we indicate the subface
4469  // which is located at the lower left
4470  // corner (the origin of the face's
4471  // local coordinate system) with
4472  // '0'. The rest of this issue is
4473  // taken care of using the above
4474  // conversion to a 'standard face
4475  // refine case')
4476  {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4477  {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4478  // cut_y (the same applies as for
4479  // cut_x)
4480  {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4481  {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4482  // cut_xyz: this information is
4483  // identical to the information
4484  // returned by
4485  // GeometryInfo<3>::real_to_standard_face_vertex()
4486  {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4487  // face_rotation=false, subfaces 0,1,2,3
4488  {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4489  // face_rotation=true, subfaces 0,1,2,3
4490  {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4491  // face_rotation=false, subfaces 0,1,2,3
4492  {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4493  // face_rotation=true, subfaces 0,1,2,3
4494  {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4495  // face_rotation=false, subfaces 0,1,2,3
4496  {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4497  // face_rotation=true, subfaces 0,1,2,3
4498  {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4499  // face_rotation=false, subfaces 0,1,2,3
4500  {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4501  // face_rotation=true, subfaces 0,1,2,3
4502 
4503  const unsigned int std_subface =
4504  subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4505  [subface];
4506  Assert(std_subface != e, ExcInternalError());
4507 
4508  // third step: these are the children, which
4509  // can be found at the given subfaces of an
4510  // isotropically refined (standard) face
4511  //
4512  // first index: (refinement_case-1)
4513  // second index: face_index
4514  // third index: subface_index (isotropic refinement)
4515  const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4516  [max_children_per_face] = {
4517  // cut_x
4518  {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4519  {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4520  {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4521  {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4522  {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4523  {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4524  // cut_y
4525  {{0, 1, 0, 1},
4526  {0, 1, 0, 1},
4527  {0, 0, 0, 0},
4528  {1, 1, 1, 1},
4529  {0, 0, 1, 1},
4530  {0, 0, 1, 1}},
4531  // cut_xy
4532  {{0, 2, 0, 2},
4533  {1, 3, 1, 3},
4534  {0, 0, 1, 1},
4535  {2, 2, 3, 3},
4536  {0, 1, 2, 3},
4537  {0, 1, 2, 3}},
4538  // cut_z
4539  {{0, 0, 1, 1},
4540  {0, 0, 1, 1},
4541  {0, 1, 0, 1},
4542  {0, 1, 0, 1},
4543  {0, 0, 0, 0},
4544  {1, 1, 1, 1}},
4545  // cut_xz
4546  {{0, 0, 1, 1},
4547  {2, 2, 3, 3},
4548  {0, 1, 2, 3},
4549  {0, 1, 2, 3},
4550  {0, 2, 0, 2},
4551  {1, 3, 1, 3}},
4552  // cut_yz
4553  {{0, 1, 2, 3},
4554  {0, 1, 2, 3},
4555  {0, 2, 0, 2},
4556  {1, 3, 1, 3},
4557  {0, 0, 1, 1},
4558  {2, 2, 3, 3}},
4559  // cut_xyz
4560  {{0, 2, 4, 6},
4561  {1, 3, 5, 7},
4562  {0, 4, 1, 5},
4563  {2, 6, 3, 7},
4564  {0, 1, 2, 3},
4565  {4, 5, 6, 7}}};
4566 
4567  // forth step: check, whether the given face
4568  // refine case is valid for the given cell
4569  // refine case. this is the case, if the
4570  // given face refine case is at least as
4571  // refined as the face is for the given cell
4572  // refine case
4573 
4574  // note, that we are considering standard
4575  // face refinement cases here and thus must
4576  // not pass the given orientation, flip and
4577  // rotation flags
4578  if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4579  face_refinement_case(ref_case, face))
4580  {
4581  // all is fine. for anisotropic face
4582  // refine cases, select one of the
4583  // isotropic subfaces which neighbors the
4584  // same child
4585 
4586  // first index: (standard) face refine case
4587  // second index: subface index
4588  const unsigned int equivalent_iso_subface[4][4] = {
4589  {0, e, e, e}, // no_refinement
4590  {0, 3, e, e}, // cut_x
4591  {0, 3, e, e}, // cut_y
4592  {0, 1, 2, 3}}; // cut_xy
4593 
4594  const unsigned int equ_std_subface =
4595  equivalent_iso_subface[std_face_ref][std_subface];
4596  Assert(equ_std_subface != e, ExcInternalError());
4597 
4598  return iso_children[ref_case - 1][face][equ_std_subface];
4599  }
4600  else
4601  {
4602  // the face_ref_case was too coarse,
4603  // throw an error
4604  Assert(false,
4605  ExcMessage("The face RefineCase is too coarse "
4606  "for the given cell RefineCase."));
4607  }
4608  // we only get here in case of an error
4609  return e;
4610 }
4611 
4612 
4613 
4614 template <>
4615 inline unsigned int
4617  const unsigned int,
4618  const unsigned int,
4619  const bool,
4620  const bool,
4621  const bool,
4622  const RefinementCase<3> &)
4623 {
4624  Assert(false, ExcNotImplemented());
4626 }
4627 
4628 
4629 
4630 template <>
4631 inline unsigned int
4632 GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4633  const unsigned int line,
4634  const bool,
4635  const bool,
4636  const bool)
4637 {
4638  (void)line;
4639  AssertIndexRange(face, faces_per_cell);
4640  AssertIndexRange(line, lines_per_face);
4641 
4642  // The face is a line itself.
4643  return face;
4644 }
4645 
4646 
4647 
4648 template <>
4649 inline unsigned int
4650 GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4651  const unsigned int line,
4652  const bool face_orientation,
4653  const bool face_flip,
4654  const bool face_rotation)
4655 {
4656  AssertIndexRange(face, faces_per_cell);
4657  AssertIndexRange(line, lines_per_face);
4658 
4659  const unsigned lines[faces_per_cell][lines_per_face] = {
4660  {8, 10, 0, 4}, // left face
4661  {9, 11, 1, 5}, // right face
4662  {2, 6, 8, 9}, // front face
4663  {3, 7, 10, 11}, // back face
4664  {0, 1, 2, 3}, // bottom face
4665  {4, 5, 6, 7}}; // top face
4666  return lines[face][real_to_standard_face_line(
4667  line, face_orientation, face_flip, face_rotation)];
4668 }
4669 
4670 
4671 
4672 inline unsigned int
4673 GeometryInfo<0>::face_to_cell_lines(const unsigned int,
4674  const unsigned int,
4675  const bool,
4676  const bool,
4677  const bool)
4678 {
4679  Assert(false, ExcNotImplemented());
4681 }
4682 
4683 
4684 
4685 template <int dim>
4686 inline unsigned int
4687 GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4688  const unsigned int,
4689  const bool,
4690  const bool,
4691  const bool)
4692 {
4693  Assert(false, ExcNotImplemented());
4695 }
4696 
4697 
4698 
4699 template <int dim>
4700 inline unsigned int
4701 GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4702  const unsigned int vertex,
4703  const bool face_orientation,
4704  const bool face_flip,
4705  const bool face_rotation)
4706 {
4707  return child_cell_on_face(RefinementCase<dim>::isotropic_refinement,
4708  face,
4709  vertex,
4710  face_orientation,
4711  face_flip,
4712  face_rotation);
4713 }
4714 
4715 
4716 
4717 inline unsigned int
4718 GeometryInfo<0>::face_to_cell_vertices(const unsigned int,
4719  const unsigned int,
4720  const bool,
4721  const bool,
4722  const bool)
4723 {
4724  Assert(false, ExcNotImplemented());
4726 }
4727 
4728 
4729 
4730 template <int dim>
4731 template <typename Number>
4732 inline Point<dim, Number>
4734 {
4736  for (unsigned int i = 0; i < dim; ++i)
4737  p[i] = std::min(std::max(q[i], Number(0.)), Number(1.));
4738 
4739  return p;
4740 }
4741 
4742 
4743 
4744 template <int dim>
4745 inline double
4747 {
4748  double result = 0.0;
4749 
4750  for (unsigned int i = 0; i < dim; ++i)
4751  {
4752  result = std::max(result, -p[i]);
4753  result = std::max(result, p[i] - 1.);
4754  }
4755 
4756  return result;
4757 }
4758 
4759 
4760 
4761 template <int dim>
4762 inline double
4764  const unsigned int i)
4765 {
4767 
4768  switch (dim)
4769  {
4770  case 1:
4771  {
4772  const double x = xi[0];
4773  switch (i)
4774  {
4775  case 0:
4776  return 1 - x;
4777  case 1:
4778  return x;
4779  }
4780  break;
4781  }
4782 
4783  case 2:
4784  {
4785  const double x = xi[0];
4786  const double y = xi[1];
4787  switch (i)
4788  {
4789  case 0:
4790  return (1 - x) * (1 - y);
4791  case 1:
4792  return x * (1 - y);
4793  case 2:
4794  return (1 - x) * y;
4795  case 3:
4796  return x * y;
4797  }
4798  break;
4799  }
4800 
4801  case 3:
4802  {
4803  const double x = xi[0];
4804  const double y = xi[1];
4805  const double z = xi[2];
4806  switch (i)
4807  {
4808  case 0:
4809  return (1 - x) * (1 - y) * (1 - z);
4810  case 1:
4811  return x * (1 - y) * (1 - z);
4812  case 2:
4813  return (1 - x) * y * (1 - z);
4814  case 3:
4815  return x * y * (1 - z);
4816  case 4:
4817  return (1 - x) * (1 - y) * z;
4818  case 5:
4819  return x * (1 - y) * z;
4820  case 6:
4821  return (1 - x) * y * z;
4822  case 7:
4823  return x * y * z;
4824  }
4825  break;
4826  }
4827 
4828  default:
4829  Assert(false, ExcNotImplemented());
4830  }
4831  return -1e9;
4832 }
4833 
4834 
4835 
4836 template <>
4838  const Point<1> &,
4839  const unsigned int i)
4840 {
4842 
4843  switch (i)
4844  {
4845  case 0:
4846  return Point<1>(-1.);
4847  case 1:
4848  return Point<1>(1.);
4849  }
4850 
4851  return Point<1>(-1e9);
4852 }
4853 
4854 
4855 
4856 template <>
4858  const Point<2> & xi,
4859  const unsigned int i)
4860 {
4862 
4863  const double x = xi[0];
4864  const double y = xi[1];
4865  switch (i)
4866  {
4867  case 0:
4868  return Point<2>(-(1 - y), -(1 - x));
4869  case 1:
4870  return Point<2>(1 - y, -x);
4871  case 2:
4872  return Point<2>(-y, 1 - x);
4873  case 3:
4874  return Point<2>(y, x);
4875  }
4876  return Point<2>(-1e9, -1e9);
4877 }
4878 
4879 
4880 
4881 template <>
4883  const Point<3> & xi,
4884  const unsigned int i)
4885 {
4887 
4888  const double x = xi[0];
4889  const double y = xi[1];
4890  const double z = xi[2];
4891  switch (i)
4892  {
4893  case 0:
4894  return Point<3>(-(1 - y) * (1 - z),
4895  -(1 - x) * (1 - z),
4896  -(1 - x) * (1 - y));
4897  case 1:
4898  return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4899  case 2:
4900  return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4901  case 3:
4902  return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4903  case 4:
4904  return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4905  case 5:
4906  return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4907  case 6:
4908  return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4909  case 7:
4910  return Point<3>(y * z, x * z, x * y);
4911  }
4912 
4913  return Point<3>(-1e9, -1e9, -1e9);
4914 }
4915 
4916 
4917 
4918 template <int dim>
4919 inline Tensor<1, dim>
4921  const unsigned int)
4922 {
4923  Assert(false, ExcNotImplemented());
4924  return Tensor<1, dim>();
4925 }
4926 
4927 
4928 unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4929 {
4930  return 0;
4931 }
4932 
4933 
4934 namespace internal
4935 {
4936  namespace GeometryInfoHelper
4937  {
4938  // wedge product of a single
4939  // vector in 2d: we just have to
4940  // rotate it by 90 degrees to the
4941  // right
4942  inline Tensor<1, 2>
4943  wedge_product(const Tensor<1, 2> (&derivative)[1])
4944  {
4945  Tensor<1, 2> result;
4946  result[0] = derivative[0][1];
4947  result[1] = -derivative[0][0];
4948 
4949  return result;
4950  }
4951 
4952 
4953  // wedge product of 2 vectors in
4954  // 3d is the cross product
4955  inline Tensor<1, 3>
4956  wedge_product(const Tensor<1, 3> (&derivative)[2])
4957  {
4958  return cross_product_3d(derivative[0], derivative[1]);
4959  }
4960 
4961 
4962  // wedge product of dim vectors
4963  // in dim-d: that's the
4964  // determinant of the matrix
4965  template <int dim>
4966  inline Tensor<0, dim>
4967  wedge_product(const Tensor<1, dim> (&derivative)[dim])
4968  {
4969  Tensor<2, dim> jacobian;
4970  for (unsigned int i = 0; i < dim; ++i)
4971  jacobian[i] = derivative[i];
4972 
4973  return determinant(jacobian);
4974  }
4975  } // namespace GeometryInfoHelper
4976 } // namespace internal
4977 
4978 
4979 template <int dim>
4980 template <int spacedim>
4981 inline void
4983 # ifndef DEAL_II_CXX14_CONSTEXPR_BUG
4984  (const Point<spacedim> (&vertices)[vertices_per_cell],
4985  Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell])
4986 # else
4987  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms)
4988 # endif
4989 {
4990  // for each of the vertices,
4991  // compute the alternating form
4992  // of the mapped unit
4993  // vectors. consider for
4994  // example the case of a quad
4995  // in spacedim==3: to do so, we
4996  // need to see how the
4997  // infinitesimal vectors
4998  // (d\xi_1,0) and (0,d\xi_2)
4999  // are transformed into
5000  // spacedim-dimensional space
5001  // and then form their cross
5002  // product (i.e. the wedge product
5003  // of two vectors). to this end, note
5004  // that
5005  // \vec x = sum_i \vec v_i phi_i(\vec xi)
5006  // so the transformed vectors are
5007  // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
5008  // and
5009  // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
5010  // which boils down to the columns
5011  // of the 3x2 matrix \grad_\xi x(\xi)
5012  //
5013  // a similar reasoning would
5014  // hold for all dim,spacedim
5015  // pairs -- we only have to
5016  // compute the wedge product of
5017  // the columns of the
5018  // derivatives
5019  for (unsigned int i = 0; i < vertices_per_cell; ++i)
5020  {
5021  Tensor<1, spacedim> derivatives[dim];
5022 
5023  for (unsigned int j = 0; j < vertices_per_cell; ++j)
5024  {
5025  const Tensor<1, dim> grad_phi_j =
5026  d_linear_shape_function_gradient(unit_cell_vertex(i), j);
5027  for (unsigned int l = 0; l < dim; ++l)
5028  derivatives[l] += vertices[j] * grad_phi_j[l];
5029  }
5030 
5031  forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
5032  }
5033 }
5034 
5035 #endif // DOXYGEN
5037 
5038 #endif
GeometryPrimitive(const unsigned int object_dimension)
GeometryPrimitive(const Object object)
Definition: point.h:111
RefinementCase operator|(const RefinementCase &r) const
static std::size_t memory_consumption()
RefinementCase operator~() const
RefinementCase operator&(const RefinementCase &r) const
RefinementCase(const typename RefinementPossibilities< dim >::Possibilities refinement_case)
static RefinementCase cut_axis(const unsigned int i)
std::uint8_t value
void serialize(Archive &ar, const unsigned int version)
RefinementCase(const std::uint8_t refinement_case)
Definition: tensor.h:503
static constexpr std::size_t memory_consumption()
SubfaceCase(const typename SubfacePossibilities< dim >::Possibilities subface_possibility)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:457
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:458
Point< 3 > vertices[4]
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInvalidCoordinate(double arg1)
static ::ExceptionBase & ExcInvalidSubface(int arg1, int arg2, int arg3)
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcInvalidSubfaceCase(int arg1)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:555
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
static ::ExceptionBase & ExcInvalidRefinementCase(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
static const char U
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:189
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static const unsigned int invalid_unsigned_int
Definition: types.h:206
boost::integer_range< IncrementableType > iota_view
Definition: iota_view.h:46
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108
static unsigned int n_children(const RefinementCase< 0 > &refinement_case)
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 const std::array< unsigned int, vertices_per_cell > dx_to_deal
static std::array< unsigned int, vertices_per_cell > vertex_indices()
static std::array< unsigned int, 0 > face_indices()
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 const std::array< unsigned int, vertices_per_cell > ucd_to_deal
static unsigned int child_cell_from_point(const Point< dim > &p)
static unsigned int standard_to_real_line_vertex(const unsigned int vertex, const bool line_orientation=true)
static constexpr std::array< unsigned int, vertices_per_cell > ucd_to_deal
static constexpr ndarray< Tensor< 1, dim >, faces_per_cell, dim - 1 > unit_tangential_vectors
static bool is_inside_unit_cell(const Point< dim > &p, const double eps)
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 double distance_to_unit_cell(const Point< dim > &p)
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 constexpr unsigned int quads_per_face
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 std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
static constexpr std::array< unsigned int, faces_per_cell > opposite_face
static constexpr unsigned int max_children_per_face
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 constexpr unsigned int vertices_per_cell
static std::array< unsigned int, 2 > standard_hex_line_to_quad_line_index(const unsigned int line)
static Point< dim > unit_cell_vertex(const unsigned int vertex)
static constexpr std::array< unsigned int, vertices_per_cell > dx_to_deal
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
static constexpr unsigned int lines_per_cell
static RefinementCase< 1 > line_refinement_case(const RefinementCase< dim > &cell_refinement_case, const unsigned int line_no)
static constexpr std::array< unsigned int, faces_per_cell > unit_normal_direction
static RefinementCase< dim > min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
static constexpr unsigned int faces_per_cell
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 double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
static constexpr unsigned int hexes_per_cell
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 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 double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
static constexpr std::array< Tensor< 1, dim >, faces_per_cell > unit_normal_vector
static Point< dim, Number > project_to_unit_cell(const Point< dim, Number > &p)
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
static constexpr unsigned int vertices_per_face
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 unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
static bool is_inside_unit_cell(const Point< dim > &p)
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
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)
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 std::array< unsigned int, 2 > standard_quad_vertex_to_line_vertex_index(const unsigned int vertex)
static constexpr ndarray< unsigned int, vertices_per_cell, dim > vertex_to_face
static void alternating_form_at_vertices(const Point< spacedim >(&vertices)[vertices_per_cell], Tensor< spacedim - dim, spacedim >(&forms)[vertices_per_cell])
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 std::array< int, faces_per_cell > unit_normal_orientation
static constexpr unsigned int quads_per_cell
static constexpr unsigned int max_children_per_cell
static unsigned int n_subfaces(const internal::SubfaceCase< dim > &subface_case)
static constexpr unsigned int lines_per_face
static std::array< unsigned int, 2 > standard_hex_vertex_to_quad_vertex_index(const unsigned int vertex)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)