Reference documentation for deal.II version Git 4abc4a1666 2020-07-04 19:58:34 +0200
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
geometry_info.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 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/point.h>
25 
26 #include <array>
27 #include <cstdint>
28 
29 
30 
32 
33 #ifndef DOXYGEN
34 namespace internal
35 {
36  namespace GeometryInfoHelper
37  {
38  // A struct that holds the values for all the arrays we want to initialize
39  // in GeometryInfo
40  template <int dim>
41  struct Initializers;
42 
43  template <>
44  struct Initializers<1>
45  {
46  static constexpr std::array<unsigned int, 2>
47  ucd_to_deal()
48  {
49  return {{0, 1}};
50  }
51 
52  static constexpr std::array<unsigned int, 2>
53  unit_normal_direction()
54  {
55  return {{0, 0}};
56  }
57 
58  static constexpr std::array<int, 2>
59  unit_normal_orientation()
60  {
61  return {{-1, 1}};
62  }
63 
64  static constexpr std::array<Tensor<1, 1>, 2>
65  unit_normal_vector()
66  {
67  return {{Tensor<1, 1>{{-1}}, Tensor<1, 1>{{1}}}};
68  }
69 
70  static constexpr std::array<std::array<Tensor<1, 1>, 0>, 2>
71  unit_tangential_vectors()
72  {
73  return {{{{}}, {{}}}};
74  }
75 
76  static constexpr std::array<unsigned int, 2>
77  opposite_face()
78  {
79  return {{1, 0}};
80  }
81 
82  static constexpr std::array<unsigned int, 2>
83  dx_to_deal()
84  {
85  return {{0, 1}};
86  }
87 
88  static constexpr std::array<std::array<unsigned int, 1>, 2>
89  vertex_to_face()
90  {
91  return {{{{0}}, {{1}}}};
92  }
93  };
94 
95  template <>
96  struct Initializers<2>
97  {
98  static constexpr std::array<unsigned int, 4>
99  ucd_to_deal()
100  {
101  return {{0, 1, 3, 2}};
102  }
103 
104  static constexpr std::array<unsigned int, 4>
105  unit_normal_direction()
106  {
107  return {{0, 0, 1, 1}};
108  }
109 
110  static constexpr std::array<int, 4>
111  unit_normal_orientation()
112  {
113  return {{-1, 1, -1, 1}};
114  }
115 
116  static constexpr std::array<Tensor<1, 2>, 4>
117  unit_normal_vector()
118  {
119  return {{Tensor<1, 2>{{-1., 0.}},
120  Tensor<1, 2>{{1., 0.}},
121  Tensor<1, 2>{{0., -1.}},
122  Tensor<1, 2>{{0., 1.}}}};
123  }
124 
125  static constexpr std::array<std::array<Tensor<1, 2>, 1>, 4>
126  unit_tangential_vectors()
127  {
128  return {{{{Tensor<1, 2>{{0, -1}}}},
129  {{Tensor<1, 2>{{0, 1}}}},
130  {{Tensor<1, 2>{{1, 0}}}},
131  {{Tensor<1, 2>{{-1, 0}}}}}};
132  }
133 
134  static constexpr std::array<unsigned int, 4>
135  opposite_face()
136  {
137  return {{1, 0, 3, 2}};
138  }
139 
140  static constexpr std::array<unsigned int, 4>
141  dx_to_deal()
142  {
143  return {{0, 2, 1, 3}};
144  }
145 
146  static constexpr std::array<std::array<unsigned int, 2>, 4>
147  vertex_to_face()
148  {
149  return {{{{0, 2}}, {{1, 2}}, {{0, 3}}, {{1, 3}}}};
150  }
151  };
152 
153  template <>
154  struct Initializers<3>
155  {
156  static constexpr std::array<unsigned int, 8>
157  ucd_to_deal()
158  {
159  return {{0, 1, 5, 4, 2, 3, 7, 6}};
160  }
161 
162  static constexpr std::array<unsigned int, 6>
163  unit_normal_direction()
164  {
165  return {{0, 0, 1, 1, 2, 2}};
166  }
167 
168  static constexpr std::array<int, 6>
169  unit_normal_orientation()
170  {
171  return {{-1, 1, -1, 1, -1, 1}};
172  }
173 
174  static constexpr std::array<Tensor<1, 3>, 6>
175  unit_normal_vector()
176  {
177  return {{Tensor<1, 3>{{-1, 0, 0}},
178  Tensor<1, 3>{{1, 0, 0}},
179  Tensor<1, 3>{{0, -1, 0}},
180  Tensor<1, 3>{{0, 1, 0}},
181  Tensor<1, 3>{{0, 0, -1}},
182  Tensor<1, 3>{{0, 0, 1}}}};
183  }
184 
185  static constexpr std::array<std::array<Tensor<1, 3>, 2>, 6>
186  unit_tangential_vectors()
187  {
188  return {{{{Tensor<1, 3>{{0, -1, 0}}, Tensor<1, 3>{{0, 0, 1}}}},
189  {{Tensor<1, 3>{{0, 1, 0}}, Tensor<1, 3>{{0, 0, 1}}}},
190  {{Tensor<1, 3>{{0, 0, -1}}, Tensor<1, 3>{{1, 0, 0}}}},
191  {{Tensor<1, 3>{{0, 0, 1}}, Tensor<1, 3>{{1, 0, 0}}}},
192  {{Tensor<1, 3>{{-1, 0, 0}}, Tensor<1, 3>{{0, 1, 0}}}},
193  {{Tensor<1, 3>{{1, 0, 0}}, Tensor<1, 3>{{0, 1, 0}}}}}};
194  }
195 
196  static constexpr std::array<unsigned int, 6>
197  opposite_face()
198  {
199  return {{1, 0, 3, 2, 5, 4}};
200  }
201 
202  static constexpr std::array<unsigned int, 8>
203  dx_to_deal()
204  {
205  return {{0, 4, 2, 6, 1, 5, 3, 7}};
206  }
207 
208  static constexpr std::array<std::array<unsigned int, 3>, 8>
209  vertex_to_face()
210  {
211  return {{{{0, 2, 4}},
212  {{1, 2, 4}},
213  {{0, 3, 4}},
214  {{1, 3, 4}},
215  {{0, 2, 5}},
216  {{1, 2, 5}},
217  {{0, 3, 5}},
218  {{1, 3, 5}}}};
219  }
220  };
221 
222  template <>
223  struct Initializers<4>
224  {
225  static constexpr std::array<unsigned int, 16>
226  ucd_to_deal()
227  {
243  numbers::invalid_unsigned_int}};
244  }
245 
246  static constexpr std::array<unsigned int, 8>
247  unit_normal_direction()
248  {
249  return {{0, 0, 1, 1, 2, 2, 3, 3}};
250  }
251 
252  static constexpr std::array<int, 8>
253  unit_normal_orientation()
254  {
255  return {{-1, 1, -1, 1, -1, 1, -1, 1}};
256  }
257 
258  static constexpr std::array<Tensor<1, 4>, 8>
259  unit_normal_vector()
260  {
261  return {{Tensor<1, 4>{{-1, 0, 0, 0}},
262  Tensor<1, 4>{{1, 0, 0, 0}},
263  Tensor<1, 4>{{0, -1, 0, 0}},
264  Tensor<1, 4>{{0, 1, 0, 0}},
265  Tensor<1, 4>{{0, 0, -1, 0}},
266  Tensor<1, 4>{{0, 0, 1, 0}},
267  Tensor<1, 4>{{0, 0, 0, -1}},
268  Tensor<1, 4>{{0, 0, 0, 1}}}};
269  }
270 
271  static constexpr std::array<std::array<Tensor<1, 4>, 3>, 8>
272  unit_tangential_vectors()
273  {
274  return {{{{Tensor<1, 4>{{0, -1, 0, 0}},
275  Tensor<1, 4>{{0, 0, 1, 0}},
276  Tensor<1, 4>{{0, 0, 0, 1}}}},
277  {{Tensor<1, 4>{{0, 1, 0, 0}},
278  Tensor<1, 4>{{0, 0, 1, 0}},
279  Tensor<1, 4>{{0, 0, 0, 1}}}},
280  {{Tensor<1, 4>{{0, 0, -1, 0}},
281  Tensor<1, 4>{{0, 0, 0, 1}},
282  Tensor<1, 4>{{1, 0, 0, 0}}}},
283  {{Tensor<1, 4>{{0, 0, 1, 0}},
284  Tensor<1, 4>{{0, 0, 0, 1}},
285  Tensor<1, 4>{{1, 0, 0, 0}}}},
286  {{Tensor<1, 4>{{0, 0, 0, -1}},
287  Tensor<1, 4>{{1, 0, 0, 0}},
288  Tensor<1, 4>{{0, 1, 0, 0}}}},
289  {{Tensor<1, 4>{{0, 0, 0, 1}},
290  Tensor<1, 4>{{1, 0, 0, 0}},
291  Tensor<1, 4>{{0, 1, 0, 0}}}},
292  {{Tensor<1, 4>{{-1, 0, 0, 0}},
293  Tensor<1, 4>{{0, 1, 0, 0}},
294  Tensor<1, 4>{{0, 0, 1, 0}}}},
295  {{Tensor<1, 4>{{1, 0, 0, 0}},
296  Tensor<1, 4>{{0, 1, 0, 0}},
297  Tensor<1, 4>{{0, 0, 1, 0}}}}}};
298  }
299 
300  static constexpr std::array<unsigned int, 8>
301  opposite_face()
302  {
303  return {{1, 0, 3, 2, 5, 4, 7, 6}};
304  }
305 
306  static constexpr std::array<unsigned int, 16>
307  dx_to_deal()
308  {
324  numbers::invalid_unsigned_int}};
325  }
326 
327  static constexpr std::array<std::array<unsigned int, 4>, 16>
328  vertex_to_face()
329  {
333  numbers::invalid_unsigned_int}},
337  numbers::invalid_unsigned_int}},
341  numbers::invalid_unsigned_int}},
345  numbers::invalid_unsigned_int}},
349  numbers::invalid_unsigned_int}},
353  numbers::invalid_unsigned_int}},
357  numbers::invalid_unsigned_int}},
361  numbers::invalid_unsigned_int}},
365  numbers::invalid_unsigned_int}},
369  numbers::invalid_unsigned_int}},
373  numbers::invalid_unsigned_int}},
377  numbers::invalid_unsigned_int}},
381  numbers::invalid_unsigned_int}},
385  numbers::invalid_unsigned_int}},
389  numbers::invalid_unsigned_int}},
393  numbers::invalid_unsigned_int}}}};
394  }
395  };
396  } // namespace GeometryInfoHelper
397 } // namespace internal
398 #endif // DOXYGEN
399 
400 
417 {
418 public:
425  enum Object
426  {
430  vertex = 0,
434  line = 1,
438  quad = 2,
442  hex = 3
443  };
444 
449  GeometryPrimitive(const Object object);
450 
456  GeometryPrimitive(const unsigned int object_dimension);
457 
462  operator unsigned int() const;
463 
464 private:
469 };
470 
471 
472 
485 template <int dim>
487 {
523  {
527  no_refinement = 0,
528 
540  isotropic_refinement = 0xFF
541  };
542 };
543 
544 
545 
555 template <>
557 {
593  {
597  no_refinement = 0,
601  cut_x = 1,
605  isotropic_refinement = cut_x
606  };
607 };
608 
609 
610 
621 template <>
623 {
659  {
663  no_refinement = 0,
667  cut_x = 1,
671  cut_y = 2,
675  cut_xy = cut_x | cut_y,
676 
680  isotropic_refinement = cut_xy
681  };
682 };
683 
684 
685 
696 template <>
698 {
734  {
738  no_refinement = 0,
742  cut_x = 1,
746  cut_y = 2,
750  cut_xy = cut_x | cut_y,
754  cut_z = 4,
758  cut_xz = cut_x | cut_z,
762  cut_yz = cut_y | cut_z,
766  cut_xyz = cut_x | cut_y | cut_z,
767 
771  isotropic_refinement = cut_xyz
772  };
773 };
774 
775 
776 
787 template <int dim>
789 {
790 public:
794  RefinementCase();
795 
801  const typename RefinementPossibilities<dim>::Possibilities refinement_case);
802 
808  explicit RefinementCase(const std::uint8_t refinement_case);
809 
821  operator std::uint8_t() const;
822 
828  operator|(const RefinementCase &r) const;
829 
834  RefinementCase operator&(const RefinementCase &r) const;
835 
844  operator~() const;
845 
846 
852  static RefinementCase
853  cut_axis(const unsigned int i);
854 
858  static std::size_t
860 
865  template <class Archive>
866  void
867  serialize(Archive &ar, const unsigned int version);
868 
873  ExcInvalidRefinementCase,
874  int,
875  << "The refinement flags given (" << arg1
876  << ") contain set bits that do not "
877  << "make sense for the space dimension of the object to which they are applied.");
878 
879 private:
884  std::uint8_t value : (dim > 0 ? dim : 1);
885 };
886 
887 
888 namespace internal
889 {
909  template <int dim>
911  {
916  {
920  case_none = 0,
921 
925  case_isotropic = static_cast<std::uint8_t>(-1)
926  };
927  };
928 
929 
938  template <>
940  {
947  {
951  case_none = 0,
952 
956  case_isotropic = case_none
957  };
958  };
959 
960 
961 
971  template <>
973  {
980  {
984  case_none = 0,
985 
989  case_isotropic = case_none
990  };
991  };
992 
993 
994 
1005  template <>
1007  {
1015  {
1019  case_none = 0,
1023  case_x = 1,
1027  case_isotropic = case_x
1028  };
1029  };
1030 
1031 
1032 
1124  template <>
1126  {
1134  {
1135  case_none = 0,
1136  case_x = 1,
1137  case_x1y = 2,
1138  case_x2y = 3,
1139  case_x1y2y = 4,
1140  case_y = 5,
1141  case_y1x = 6,
1142  case_y2x = 7,
1143  case_y1x2x = 8,
1144  case_xy = 9,
1145 
1146  case_isotropic = case_xy
1147  };
1148  };
1149 
1150 
1151 
1158  template <int dim>
1159  class SubfaceCase : public SubfacePossibilities<dim>
1160  {
1161  public:
1168  subface_possibility);
1169 
1181  operator std::uint8_t() const;
1182 
1186  static constexpr std::size_t
1188 
1193  ExcInvalidSubfaceCase,
1194  int,
1195  << "The subface case given (" << arg1 << ") does not make sense "
1196  << "for the space dimension of the object to which they are applied.");
1197 
1198  private:
1203  std::uint8_t value : (dim == 3 ? 4 : 1);
1204  };
1205 
1206 } // namespace internal
1207 
1208 
1209 
1210 template <int dim>
1212 
1213 
1214 
1231 template <>
1232 struct GeometryInfo<0>
1233 {
1241  static constexpr unsigned int max_children_per_cell = 1;
1242 
1246  static constexpr unsigned int faces_per_cell = 0;
1247 
1264  static std::array<unsigned int, 0>
1265  face_indices();
1266 
1274  static constexpr unsigned int max_children_per_face = 0;
1275 
1281  static unsigned int
1282  n_children(const RefinementCase<0> &refinement_case);
1283 
1287  static constexpr unsigned int vertices_per_cell = 1;
1288 
1307  static std::array<unsigned int, vertices_per_cell>
1308  vertex_indices();
1309 
1316  static constexpr unsigned int vertices_per_face = 0;
1317 
1321  static constexpr unsigned int lines_per_face = 0;
1322 
1326  static constexpr unsigned int quads_per_face = 0;
1327 
1331  static constexpr unsigned int lines_per_cell = 0;
1332 
1336  static constexpr unsigned int quads_per_cell = 0;
1337 
1341  static constexpr unsigned int hexes_per_cell = 0;
1342 
1360  static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1361 
1375  static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1376 };
1377 
1378 
1379 
1903 template <int dim>
1904 struct GeometryInfo
1905 {
1913  static constexpr unsigned int max_children_per_cell = 1 << dim;
1914 
1918  static constexpr unsigned int faces_per_cell = 2 * dim;
1919 
1937  face_indices();
1938 
1946  static constexpr unsigned int max_children_per_face =
1947  GeometryInfo<dim - 1>::max_children_per_cell;
1948 
1952  static constexpr unsigned int vertices_per_cell = 1 << dim;
1953 
1970  vertex_indices();
1971 
1975  static constexpr unsigned int vertices_per_face =
1976  GeometryInfo<dim - 1>::vertices_per_cell;
1977 
1981  static constexpr unsigned int lines_per_face =
1982  GeometryInfo<dim - 1>::lines_per_cell;
1983 
1987  static constexpr unsigned int quads_per_face =
1988  GeometryInfo<dim - 1>::quads_per_cell;
1989 
1999  static constexpr unsigned int lines_per_cell =
2000  (2 * GeometryInfo<dim - 1>::lines_per_cell +
2001  GeometryInfo<dim - 1>::vertices_per_cell);
2002 
2010  static constexpr unsigned int quads_per_cell =
2011  (2 * GeometryInfo<dim - 1>::quads_per_cell +
2012  GeometryInfo<dim - 1>::lines_per_cell);
2013 
2017  static constexpr unsigned int hexes_per_cell =
2018  (2 * GeometryInfo<dim - 1>::hexes_per_cell +
2019  GeometryInfo<dim - 1>::quads_per_cell);
2020 
2038  static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
2039  internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal();
2040 
2054  static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
2055  internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal();
2056 
2067  static constexpr std::array<std::array<unsigned int, dim>, vertices_per_cell>
2068  vertex_to_face =
2069  internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face();
2070 
2075  static unsigned int
2076  n_children(const RefinementCase<dim> &refinement_case);
2077 
2082  static unsigned int
2083  n_subfaces(const internal::SubfaceCase<dim> &subface_case);
2084 
2094  static double
2095  subface_ratio(const internal::SubfaceCase<dim> &subface_case,
2096  const unsigned int subface_no);
2097 
2103  static RefinementCase<dim - 1>
2104  face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2105  const unsigned int face_no,
2106  const bool face_orientation = true,
2107  const bool face_flip = false,
2108  const bool face_rotation = false);
2109 
2115  static RefinementCase<dim>
2116  min_cell_refinement_case_for_face_refinement(
2117  const RefinementCase<dim - 1> &face_refinement_case,
2118  const unsigned int face_no,
2119  const bool face_orientation = true,
2120  const bool face_flip = false,
2121  const bool face_rotation = false);
2122 
2127  static RefinementCase<1>
2128  line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2129  const unsigned int line_no);
2130 
2135  static RefinementCase<dim>
2136  min_cell_refinement_case_for_line_refinement(const unsigned int line_no);
2137 
2184  static unsigned int
2185  child_cell_on_face(const RefinementCase<dim> & ref_case,
2186  const unsigned int face,
2187  const unsigned int subface,
2188  const bool face_orientation = true,
2189  const bool face_flip = false,
2190  const bool face_rotation = false,
2191  const RefinementCase<dim - 1> &face_refinement_case =
2193 
2207  static unsigned int
2208  line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
2209 
2230  static unsigned int
2231  face_to_cell_vertices(const unsigned int face,
2232  const unsigned int vertex,
2233  const bool face_orientation = true,
2234  const bool face_flip = false,
2235  const bool face_rotation = false);
2236 
2248  static unsigned int
2249  face_to_cell_lines(const unsigned int face,
2250  const unsigned int line,
2251  const bool face_orientation = true,
2252  const bool face_flip = false,
2253  const bool face_rotation = false);
2254 
2264  static unsigned int
2265  standard_to_real_face_vertex(const unsigned int vertex,
2266  const bool face_orientation = true,
2267  const bool face_flip = false,
2268  const bool face_rotation = false);
2269 
2279  static unsigned int
2280  real_to_standard_face_vertex(const unsigned int vertex,
2281  const bool face_orientation = true,
2282  const bool face_flip = false,
2283  const bool face_rotation = false);
2284 
2294  static unsigned int
2295  standard_to_real_face_line(const unsigned int line,
2296  const bool face_orientation = true,
2297  const bool face_flip = false,
2298  const bool face_rotation = false);
2299 
2305  static unsigned int
2306  standard_to_real_line_vertex(const unsigned int vertex,
2307  const bool line_orientation = true);
2308 
2316  static std::array<unsigned int, 2>
2317  standard_quad_vertex_to_line_vertex_index(const unsigned int vertex);
2318 
2326  static std::array<unsigned int, 2>
2327  standard_hex_vertex_to_quad_vertex_index(const unsigned int vertex);
2328 
2336  static std::array<unsigned int, 2>
2337  standard_hex_line_to_quad_line_index(const unsigned int line);
2338 
2348  static unsigned int
2349  real_to_standard_face_line(const unsigned int line,
2350  const bool face_orientation = true,
2351  const bool face_flip = false,
2352  const bool face_rotation = false);
2353 
2359  static Point<dim>
2360  unit_cell_vertex(const unsigned int vertex);
2361 
2371  static unsigned int
2372  child_cell_from_point(const Point<dim> &p);
2373 
2381  static Point<dim>
2382  cell_to_child_coordinates(const Point<dim> & p,
2383  const unsigned int child_index,
2384  const RefinementCase<dim> refine_case =
2386 
2392  static Point<dim>
2393  child_to_cell_coordinates(const Point<dim> & p,
2394  const unsigned int child_index,
2395  const RefinementCase<dim> refine_case =
2397 
2402  static bool
2403  is_inside_unit_cell(const Point<dim> &p);
2404 
2416  static bool
2417  is_inside_unit_cell(const Point<dim> &p, const double eps);
2418 
2423  static Point<dim>
2424  project_to_unit_cell(const Point<dim> &p);
2425 
2431  static double
2432  distance_to_unit_cell(const Point<dim> &p);
2433 
2438  static double
2439  d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2440 
2445  static Tensor<1, dim>
2446  d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2447 
2499  template <int spacedim>
2500  static void
2501  alternating_form_at_vertices
2502 #ifndef DEAL_II_CXX14_CONSTEXPR_BUG
2503  (const Point<spacedim> (&vertices)[vertices_per_cell],
2504  Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2505 #else
2506  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2507 #endif
2508 
2518  static constexpr std::array<unsigned int, faces_per_cell>
2519  unit_normal_direction =
2520  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction();
2521 
2538  static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2539  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation();
2540 
2551  static constexpr std::array<Tensor<1, dim>, faces_per_cell>
2552  unit_normal_vector =
2553  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_vector();
2554 
2568  static constexpr std::array<std::array<Tensor<1, dim>, dim - 1>,
2569  faces_per_cell>
2570  unit_tangential_vectors = internal::GeometryInfoHelper::Initializers<
2571  dim>::unit_tangential_vectors();
2572 
2578  static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2579  internal::GeometryInfoHelper::Initializers<dim>::opposite_face();
2580 
2581 
2585  DeclException1(ExcInvalidCoordinate,
2586  double,
2587  << "The coordinates must satisfy 0 <= x_i <= 1, "
2588  << "but here we have x_i=" << arg1);
2589 
2593  DeclException3(ExcInvalidSubface,
2594  int,
2595  int,
2596  int,
2597  << "RefinementCase<dim> " << arg1 << ": face " << arg2
2598  << " has no subface " << arg3);
2599 };
2600 
2601 
2602 
2603 #ifndef DOXYGEN
2604 
2605 
2606 /* -------------- declaration of explicit specializations ------------- */
2607 
2608 
2609 template <>
2612  const unsigned int i);
2613 template <>
2616  const unsigned int i);
2617 template <>
2620  const unsigned int i);
2621 
2622 
2623 
2624 /* -------------- inline functions ------------- */
2625 
2626 
2627 inline GeometryPrimitive::GeometryPrimitive(const Object object)
2628  : object(object)
2629 {}
2630 
2631 
2632 
2633 inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2634  : object(static_cast<Object>(object_dimension))
2635 {}
2636 
2637 
2638 inline GeometryPrimitive::operator unsigned int() const
2639 {
2640  return static_cast<unsigned int>(object);
2641 }
2642 
2643 
2644 
2645 namespace internal
2646 {
2647  template <int dim>
2648  inline SubfaceCase<dim>::SubfaceCase(
2649  const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2650  : value(subface_possibility)
2651  {}
2652 
2653 
2654  template <int dim>
2655  inline SubfaceCase<dim>::operator std::uint8_t() const
2656  {
2657  return value;
2658  }
2659 
2660 
2661 } // namespace internal
2662 
2663 
2664 template <int dim>
2665 inline RefinementCase<dim>
2666 RefinementCase<dim>::cut_axis(const unsigned int)
2667 {
2668  Assert(false, ExcInternalError());
2669  return static_cast<std::uint8_t>(-1);
2670 }
2671 
2672 
2673 template <>
2674 inline RefinementCase<1>
2675 RefinementCase<1>::cut_axis(const unsigned int i)
2676 {
2677  AssertIndexRange(i, 1);
2678 
2680  return options[i];
2681 }
2682 
2683 
2684 
2685 template <>
2686 inline RefinementCase<2>
2687 RefinementCase<2>::cut_axis(const unsigned int i)
2688 {
2689  AssertIndexRange(i, 2);
2690 
2693  return options[i];
2694 }
2695 
2696 
2697 
2698 template <>
2699 inline RefinementCase<3>
2700 RefinementCase<3>::cut_axis(const unsigned int i)
2701 {
2702  AssertIndexRange(i, 3);
2703 
2707  return options[i];
2708 }
2709 
2710 
2711 
2712 template <int dim>
2714  : value(RefinementPossibilities<dim>::no_refinement)
2715 {}
2716 
2717 
2718 
2719 template <int dim>
2721  const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2722  : value(refinement_case)
2723 {
2724  // check that only those bits of
2725  // the given argument are set that
2726  // make sense for a given space
2727  // dimension
2728  Assert((refinement_case &
2730  refinement_case,
2731  ExcInvalidRefinementCase(refinement_case));
2732 }
2733 
2734 
2735 
2736 template <int dim>
2737 inline RefinementCase<dim>::RefinementCase(const std::uint8_t refinement_case)
2738  : value(refinement_case)
2739 {
2740  // check that only those bits of
2741  // the given argument are set that
2742  // make sense for a given space
2743  // dimension
2744  Assert((refinement_case &
2746  refinement_case,
2747  ExcInvalidRefinementCase(refinement_case));
2748 }
2749 
2750 
2751 
2752 template <int dim>
2753 inline RefinementCase<dim>::operator std::uint8_t() const
2754 {
2755  return value;
2756 }
2757 
2758 
2759 
2760 template <int dim>
2761 inline RefinementCase<dim>
2763 {
2764  return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2765 }
2766 
2767 
2768 
2769 template <int dim>
2771  operator&(const RefinementCase<dim> &r) const
2772 {
2773  return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2774 }
2775 
2776 
2777 
2778 template <int dim>
2779 inline RefinementCase<dim>
2781 {
2782  return RefinementCase<dim>(static_cast<std::uint8_t>(
2784 }
2785 
2786 
2787 
2788 template <int dim>
2789 inline std::size_t
2791 {
2792  return sizeof(RefinementCase<dim>);
2793 }
2794 
2795 
2796 
2797 template <int dim>
2798 template <class Archive>
2799 inline void
2800 RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2801 {
2802  // serialization can't deal with bitfields, so copy from/to a full sized
2803  // std::uint8_t
2804  std::uint8_t uchar_value = value;
2805  ar & uchar_value;
2806  value = uchar_value;
2807 }
2808 
2809 
2810 
2811 template <>
2812 inline Point<1>
2813 GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2814 {
2816 
2817  return Point<1>(static_cast<double>(vertex));
2818 }
2819 
2820 
2821 
2822 template <>
2823 inline Point<2>
2824 GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2825 {
2827 
2828  return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2829 }
2830 
2831 
2832 
2833 template <>
2834 inline Point<3>
2835 GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2836 {
2838 
2839  return {static_cast<double>(vertex % 2),
2840  static_cast<double>(vertex / 2 % 2),
2841  static_cast<double>(vertex / 4)};
2842 }
2843 
2844 
2845 
2846 inline std::array<unsigned int, 0>
2848 {
2849  return {{}};
2850 }
2851 
2852 
2853 
2854 inline std::array<unsigned int, 1>
2856 {
2857  return {{0}};
2858 }
2859 
2860 
2861 
2862 template <int dim>
2865 {
2866  return {0U, faces_per_cell};
2867 }
2868 
2869 
2870 
2871 template <int dim>
2874 {
2875  return {0U, vertices_per_cell};
2876 }
2877 
2878 
2879 
2880 template <int dim>
2881 inline Point<dim>
2882 GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2883 {
2884  Assert(false, ExcNotImplemented());
2885 
2886  return Point<dim>();
2887 }
2888 
2889 
2890 
2891 template <>
2892 inline unsigned int
2894 {
2895  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2896 
2897  return (p[0] <= 0.5 ? 0 : 1);
2898 }
2899 
2900 
2901 
2902 template <>
2903 inline unsigned int
2905 {
2906  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2907  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2908 
2909  return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
2910 }
2911 
2912 
2913 
2914 template <>
2915 inline unsigned int
2917 {
2918  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2919  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2920  Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
2921 
2922  return (p[0] <= 0.5 ?
2923  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
2924  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
2925 }
2926 
2927 
2928 template <int dim>
2929 inline unsigned int
2931 {
2932  Assert(false, ExcNotImplemented());
2933 
2934  return 0;
2935 }
2936 
2937 
2938 
2939 template <>
2940 inline Point<1>
2942  const unsigned int child_index,
2943  const RefinementCase<1> refine_case)
2944 
2945 {
2946  AssertIndexRange(child_index, 2);
2948  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
2949 
2950  return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
2951 }
2952 
2953 
2954 
2955 template <>
2956 inline Point<2>
2958  const unsigned int child_index,
2959  const RefinementCase<2> refine_case)
2960 
2961 {
2962  AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
2963 
2964  Point<2> point = p;
2965  switch (refine_case)
2966  {
2968  point[0] *= 2.0;
2969  if (child_index == 1)
2970  point[0] -= 1.0;
2971  break;
2973  point[1] *= 2.0;
2974  if (child_index == 1)
2975  point[1] -= 1.0;
2976  break;
2978  point *= 2.0;
2979  point -= unit_cell_vertex(child_index);
2980  break;
2981  default:
2982  Assert(false, ExcInternalError());
2983  }
2984 
2985  return point;
2986 }
2987 
2988 
2989 
2990 template <>
2991 inline Point<3>
2993  const unsigned int child_index,
2994  const RefinementCase<3> refine_case)
2995 
2996 {
2997  AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
2998 
2999  Point<3> point = p;
3000  // there might be a cleverer way to do
3001  // this, but since this function is called
3002  // in very few places for initialization
3003  // purposes only, I don't care at the
3004  // moment
3005  switch (refine_case)
3006  {
3008  point[0] *= 2.0;
3009  if (child_index == 1)
3010  point[0] -= 1.0;
3011  break;
3013  point[1] *= 2.0;
3014  if (child_index == 1)
3015  point[1] -= 1.0;
3016  break;
3018  point[2] *= 2.0;
3019  if (child_index == 1)
3020  point[2] -= 1.0;
3021  break;
3023  point[0] *= 2.0;
3024  point[1] *= 2.0;
3025  if (child_index % 2 == 1)
3026  point[0] -= 1.0;
3027  if (child_index / 2 == 1)
3028  point[1] -= 1.0;
3029  break;
3031  // careful, this is slightly
3032  // different from xy and yz due to
3033  // different internal numbering of
3034  // children!
3035  point[0] *= 2.0;
3036  point[2] *= 2.0;
3037  if (child_index / 2 == 1)
3038  point[0] -= 1.0;
3039  if (child_index % 2 == 1)
3040  point[2] -= 1.0;
3041  break;
3043  point[1] *= 2.0;
3044  point[2] *= 2.0;
3045  if (child_index % 2 == 1)
3046  point[1] -= 1.0;
3047  if (child_index / 2 == 1)
3048  point[2] -= 1.0;
3049  break;
3051  point *= 2.0;
3052  point -= unit_cell_vertex(child_index);
3053  break;
3054  default:
3055  Assert(false, ExcInternalError());
3056  }
3057 
3058  return point;
3059 }
3060 
3061 
3062 
3063 template <int dim>
3064 inline Point<dim>
3066  const Point<dim> & /*p*/,
3067  const unsigned int /*child_index*/,
3068  const RefinementCase<dim> /*refine_case*/)
3069 
3070 {
3071  Assert(false, ExcNotImplemented());
3072  return Point<dim>();
3073 }
3074 
3075 
3076 
3077 template <>
3078 inline Point<1>
3080  const unsigned int child_index,
3081  const RefinementCase<1> refine_case)
3082 
3083 {
3084  AssertIndexRange(child_index, 2);
3086  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3087 
3088  return (p + unit_cell_vertex(child_index)) * 0.5;
3089 }
3090 
3091 
3092 
3093 template <>
3094 inline Point<3>
3096  const unsigned int child_index,
3097  const RefinementCase<3> refine_case)
3098 
3099 {
3100  AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3101 
3102  Point<3> point = p;
3103  // there might be a cleverer way to do
3104  // this, but since this function is called
3105  // in very few places for initialization
3106  // purposes only, I don't care at the
3107  // moment
3108  switch (refine_case)
3109  {
3111  if (child_index == 1)
3112  point[0] += 1.0;
3113  point[0] *= 0.5;
3114  break;
3116  if (child_index == 1)
3117  point[1] += 1.0;
3118  point[1] *= 0.5;
3119  break;
3121  if (child_index == 1)
3122  point[2] += 1.0;
3123  point[2] *= 0.5;
3124  break;
3126  if (child_index % 2 == 1)
3127  point[0] += 1.0;
3128  if (child_index / 2 == 1)
3129  point[1] += 1.0;
3130  point[0] *= 0.5;
3131  point[1] *= 0.5;
3132  break;
3134  // careful, this is slightly
3135  // different from xy and yz due to
3136  // different internal numbering of
3137  // children!
3138  if (child_index / 2 == 1)
3139  point[0] += 1.0;
3140  if (child_index % 2 == 1)
3141  point[2] += 1.0;
3142  point[0] *= 0.5;
3143  point[2] *= 0.5;
3144  break;
3146  if (child_index % 2 == 1)
3147  point[1] += 1.0;
3148  if (child_index / 2 == 1)
3149  point[2] += 1.0;
3150  point[1] *= 0.5;
3151  point[2] *= 0.5;
3152  break;
3154  point += unit_cell_vertex(child_index);
3155  point *= 0.5;
3156  break;
3157  default:
3158  Assert(false, ExcInternalError());
3159  }
3160 
3161  return point;
3162 }
3163 
3164 
3165 
3166 template <>
3167 inline Point<2>
3169  const unsigned int child_index,
3170  const RefinementCase<2> refine_case)
3171 {
3172  AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3173 
3174  Point<2> point = p;
3175  switch (refine_case)
3176  {
3178  if (child_index == 1)
3179  point[0] += 1.0;
3180  point[0] *= 0.5;
3181  break;
3183  if (child_index == 1)
3184  point[1] += 1.0;
3185  point[1] *= 0.5;
3186  break;
3188  point += unit_cell_vertex(child_index);
3189  point *= 0.5;
3190  break;
3191  default:
3192  Assert(false, ExcInternalError());
3193  }
3194 
3195  return point;
3196 }
3197 
3198 
3199 
3200 template <int dim>
3201 inline Point<dim>
3203  const Point<dim> & /*p*/,
3204  const unsigned int /*child_index*/,
3205  const RefinementCase<dim> /*refine_case*/)
3206 {
3207  Assert(false, ExcNotImplemented());
3208  return Point<dim>();
3209 }
3210 
3211 
3212 
3213 template <int dim>
3214 inline bool
3216 {
3217  Assert(false, ExcNotImplemented());
3218  return false;
3219 }
3220 
3221 template <>
3222 inline bool
3224 {
3225  return (p[0] >= 0.) && (p[0] <= 1.);
3226 }
3227 
3228 
3229 
3230 template <>
3231 inline bool
3233 {
3234  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
3235 }
3236 
3237 
3238 
3239 template <>
3240 inline bool
3242 {
3243  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
3244  (p[2] >= 0.) && (p[2] <= 1.);
3245 }
3246 
3247 
3248 
3249 template <int dim>
3250 inline bool
3252 {
3253  Assert(false, ExcNotImplemented());
3254  return false;
3255 }
3256 
3257 template <>
3258 inline bool
3259 GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
3260 {
3261  return (p[0] >= -eps) && (p[0] <= 1. + eps);
3262 }
3263 
3264 
3265 
3266 template <>
3267 inline bool
3268 GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
3269 {
3270  const double l = -eps, u = 1 + eps;
3271  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
3272 }
3273 
3274 
3275 
3276 template <>
3277 inline bool
3278 GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
3279 {
3280  const double l = -eps, u = 1.0 + eps;
3281  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
3282  (p[2] >= l) && (p[2] <= u);
3283 }
3284 
3285 
3286 
3287 template <>
3288 inline unsigned int
3289 GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
3290  const unsigned int vertex)
3291 {
3292  (void)line;
3294  AssertIndexRange(vertex, 2);
3295 
3296  return vertex;
3297 }
3298 
3299 
3300 template <>
3301 inline unsigned int
3302 GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
3303  const unsigned int vertex)
3304 {
3305  constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
3306  return cell_vertices[line][vertex];
3307 }
3308 
3309 
3310 
3311 template <>
3312 inline unsigned int
3313 GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
3314  const unsigned int vertex)
3315 {
3317  AssertIndexRange(vertex, 2);
3318 
3319  constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3320  {1, 3},
3321  {0, 1},
3322  {2, 3},
3323  {4, 6}, // top face
3324  {5, 7},
3325  {4, 5},
3326  {6, 7},
3327  {0,
3328  4}, // connects of bottom
3329  {1, 5}, // top face
3330  {2, 6},
3331  {3, 7}};
3332  return vertices[line][vertex];
3333 }
3334 
3335 
3336 template <>
3337 inline unsigned int
3338 GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3339 {
3340  Assert(false, ExcNotImplemented());
3342 }
3343 
3344 template <>
3345 inline unsigned int
3346 GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3347  const bool face_orientation,
3348  const bool face_flip,
3349  const bool face_rotation)
3350 {
3352 
3353  // set up a table to make sure that
3354  // we handle non-standard faces correctly
3355  //
3356  // so set up a table that for each vertex (of
3357  // a quad in standard position) describes
3358  // which vertex to take
3359  //
3360  // first index: four vertices 0...3
3361  //
3362  // second index: face_orientation; 0:
3363  // opposite normal, 1: standard
3364  //
3365  // third index: face_flip; 0: standard, 1:
3366  // face rotated by 180 degrees
3367  //
3368  // forth index: face_rotation: 0: standard,
3369  // 1: face rotated by 90 degrees
3370 
3371  constexpr unsigned int vertex_translation[4][2][2][2] = {
3372  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3373  // face_rotation=false and true
3374  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3375  // face_rotation=false and true
3376  {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3377  // face_rotation=false and true
3378  {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3379  // face_rotation=false and true
3380 
3381  {{{2, 3}, // vertex 1 ...
3382  {1, 0}},
3383  {{1, 0}, {2, 3}}},
3384 
3385  {{{1, 0}, // vertex 2 ...
3386  {2, 3}},
3387  {{2, 3}, {1, 0}}},
3388 
3389  {{{3, 1}, // vertex 3 ...
3390  {0, 2}},
3391  {{3, 1}, {0, 2}}}};
3392 
3393  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3394 }
3395 
3396 
3397 
3398 template <int dim>
3399 inline unsigned int
3400 GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3401  const bool,
3402  const bool,
3403  const bool)
3404 {
3405  Assert(dim > 1, ExcImpossibleInDim(dim));
3407  return vertex;
3408 }
3409 
3410 template <int dim>
3411 inline unsigned int
3413 {
3414  constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3415  0, 2, 2, 4, 2, 4, 4, 8};
3416 
3417  return n_children[ref_case];
3418 }
3419 
3420 
3421 
3422 template <int dim>
3423 inline unsigned int
3425 {
3426  Assert(false, ExcNotImplemented());
3427  return 0;
3428 }
3429 
3430 template <>
3431 inline unsigned int
3433 {
3434  Assert(false, ExcImpossibleInDim(1));
3435  return 0;
3436 }
3437 
3438 template <>
3439 inline unsigned int
3441 {
3442  return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3443 }
3444 
3445 
3446 
3447 template <>
3448 inline unsigned int
3450 {
3451  const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3452  0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3453  return nsubs[subface_case];
3454 }
3455 
3456 
3457 
3458 template <int dim>
3459 inline double
3461  const unsigned int)
3462 {
3463  Assert(false, ExcNotImplemented());
3464  return 0.;
3465 }
3466 
3467 template <>
3468 inline double
3470  const unsigned int)
3471 {
3472  return 1;
3473 }
3474 
3475 
3476 template <>
3477 inline double
3479  const unsigned int)
3480 {
3481  double ratio = 1;
3482  switch (subface_case)
3483  {
3485  // Here, an
3486  // Assert(false,ExcInternalError())
3487  // would be the right
3488  // choice, but
3489  // unfortunately the
3490  // current function is
3491  // also called for faces
3492  // without children (see
3493  // tests/fe/mapping.cc).
3494  // Assert(false, ExcMessage("Face has no subfaces."));
3495  // Furthermore, assign
3496  // following value as
3497  // otherwise the
3498  // bits/volume_x tests
3499  // break
3501  break;
3503  ratio = 0.5;
3504  break;
3505  default:
3506  // there should be no
3507  // cases left
3508  Assert(false, ExcInternalError());
3509  break;
3510  }
3511 
3512  return ratio;
3513 }
3514 
3515 
3516 template <>
3517 inline double
3519  const unsigned int subface_no)
3520 {
3521  double ratio = 1;
3522  switch (subface_case)
3523  {
3525  // Here, an
3526  // Assert(false,ExcInternalError())
3527  // would be the right
3528  // choice, but
3529  // unfortunately the
3530  // current function is
3531  // also called for faces
3532  // without children (see
3533  // tests/bits/mesh_3d_16.cc). Add
3534  // following switch to
3535  // avoid diffs in
3536  // tests/bits/mesh_3d_16
3538  break;
3541  ratio = 0.5;
3542  break;
3546  ratio = 0.25;
3547  break;
3550  if (subface_no < 2)
3551  ratio = 0.25;
3552  else
3553  ratio = 0.5;
3554  break;
3557  if (subface_no == 0)
3558  ratio = 0.5;
3559  else
3560  ratio = 0.25;
3561  break;
3562  default:
3563  // there should be no
3564  // cases left
3565  Assert(false, ExcInternalError());
3566  break;
3567  }
3568 
3569  return ratio;
3570 }
3571 
3572 
3573 
3574 template <int dim>
3576  const RefinementCase<dim> &,
3577  const unsigned int,
3578  const bool,
3579  const bool,
3580  const bool)
3581 {
3582  Assert(false, ExcNotImplemented());
3583  return RefinementCase<dim - 1>::no_refinement;
3584 }
3585 
3586 template <>
3588  const RefinementCase<1> &,
3589  const unsigned int,
3590  const bool,
3591  const bool,
3592  const bool)
3593 {
3594  Assert(false, ExcImpossibleInDim(1));
3595 
3597 }
3598 
3599 
3600 template <>
3601 inline RefinementCase<1>
3603  const RefinementCase<2> &cell_refinement_case,
3604  const unsigned int face_no,
3605  const bool,
3606  const bool,
3607  const bool)
3608 {
3609  const unsigned int dim = 2;
3610  AssertIndexRange(cell_refinement_case,
3613 
3614  const RefinementCase<dim - 1>
3617  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3618  RefinementCase<dim - 1>::no_refinement},
3619 
3620  {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3621 
3622  {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3623 
3624  {RefinementCase<dim - 1>::cut_x, // cut_xy
3625  RefinementCase<dim - 1>::cut_x}};
3626 
3627  return ref_cases[cell_refinement_case][face_no / 2];
3628 }
3629 
3630 
3631 template <>
3632 inline RefinementCase<2>
3634  const RefinementCase<3> &cell_refinement_case,
3635  const unsigned int face_no,
3636  const bool face_orientation,
3637  const bool /*face_flip*/,
3638  const bool face_rotation)
3639 {
3640  const unsigned int dim = 3;
3641  AssertIndexRange(cell_refinement_case,
3644 
3645  const RefinementCase<dim - 1>
3648  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3649  RefinementCase<dim - 1>::no_refinement,
3650  RefinementCase<dim - 1>::no_refinement},
3651 
3652  {RefinementCase<dim - 1>::no_refinement, // cut_x
3653  RefinementCase<dim - 1>::cut_y,
3654  RefinementCase<dim - 1>::cut_x},
3655 
3656  {RefinementCase<dim - 1>::cut_x, // cut_y
3657  RefinementCase<dim - 1>::no_refinement,
3658  RefinementCase<dim - 1>::cut_y},
3659 
3660  {RefinementCase<dim - 1>::cut_x, // cut_xy
3661  RefinementCase<dim - 1>::cut_y,
3662  RefinementCase<dim - 1>::cut_xy},
3663 
3664  {RefinementCase<dim - 1>::cut_y, // cut_z
3665  RefinementCase<dim - 1>::cut_x,
3666  RefinementCase<dim - 1>::no_refinement},
3667 
3668  {RefinementCase<dim - 1>::cut_y, // cut_xz
3669  RefinementCase<dim - 1>::cut_xy,
3670  RefinementCase<dim - 1>::cut_x},
3671 
3672  {RefinementCase<dim - 1>::cut_xy, // cut_yz
3673  RefinementCase<dim - 1>::cut_x,
3674  RefinementCase<dim - 1>::cut_y},
3675 
3676  {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3677  RefinementCase<dim - 1>::cut_xy,
3678  RefinementCase<dim - 1>::cut_xy},
3679  };
3680 
3681  const RefinementCase<dim - 1> ref_case =
3682  ref_cases[cell_refinement_case][face_no / 2];
3683 
3684  const RefinementCase<dim - 1> flip[4] = {
3685  RefinementCase<dim - 1>::no_refinement,
3686  RefinementCase<dim - 1>::cut_y,
3687  RefinementCase<dim - 1>::cut_x,
3688  RefinementCase<dim - 1>::cut_xy};
3689 
3690  // correct the ref_case for face_orientation
3691  // and face_rotation. for face_orientation,
3692  // 'true' is the default value whereas for
3693  // face_rotation, 'false' is standard. If
3694  // <tt>face_rotation==face_orientation</tt>,
3695  // then one of them is non-standard and we
3696  // have to swap cut_x and cut_y, otherwise no
3697  // change is necessary. face_flip has no
3698  // influence. however, in order to keep the
3699  // interface consistent with other functions,
3700  // we still include it as an argument to this
3701  // function
3702  return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3703 }
3704 
3705 
3706 
3707 template <int dim>
3708 inline RefinementCase<1>
3710  const unsigned int)
3711 {
3712  Assert(false, ExcNotImplemented());
3714 }
3715 
3716 template <>
3717 inline RefinementCase<1>
3719  const RefinementCase<1> &cell_refinement_case,
3720  const unsigned int line_no)
3721 {
3722  (void)line_no;
3723  const unsigned int dim = 1;
3724  (void)dim;
3725  AssertIndexRange(cell_refinement_case,
3728 
3729  return cell_refinement_case;
3730 }
3731 
3732 
3733 template <>
3734 inline RefinementCase<1>
3736  const RefinementCase<2> &cell_refinement_case,
3737  const unsigned int line_no)
3738 {
3739  // Assertions are in face_refinement_case()
3740  return face_refinement_case(cell_refinement_case, line_no);
3741 }
3742 
3743 
3744 template <>
3745 inline RefinementCase<1>
3747  const RefinementCase<3> &cell_refinement_case,
3748  const unsigned int line_no)
3749 {
3750  const unsigned int dim = 3;
3751  AssertIndexRange(cell_refinement_case,
3754 
3755  // array indicating, which simple refine
3756  // case cuts a line in direction x, y or
3757  // z. For example, cut_y and everything
3758  // containing cut_y (cut_xy, cut_yz,
3759  // cut_xyz) cuts lines, which are in y
3760  // direction.
3761  const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3764 
3765  // order the direction of lines
3766  // 0->x, 1->y, 2->z
3767  const unsigned int direction[lines_per_cell] = {
3768  1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3769 
3770  return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3773 }
3774 
3775 
3776 
3777 template <int dim>
3778 inline RefinementCase<dim>
3780  const RefinementCase<dim - 1> &,
3781  const unsigned int,
3782  const bool,
3783  const bool,
3784  const bool)
3785 {
3786  Assert(false, ExcNotImplemented());
3787 
3789 }
3790 
3791 template <>
3792 inline RefinementCase<1>
3794  const RefinementCase<0> &,
3795  const unsigned int,
3796  const bool,
3797  const bool,
3798  const bool)
3799 {
3800  const unsigned int dim = 1;
3801  Assert(false, ExcImpossibleInDim(dim));
3802 
3804 }
3805 
3806 
3807 template <>
3808 inline RefinementCase<2>
3811  const unsigned int face_no,
3812  const bool,
3813  const bool,
3814  const bool)
3815 {
3816  const unsigned int dim = 2;
3817  AssertIndexRange(face_refinement_case,
3820 
3821  if (face_refinement_case == RefinementCase<dim>::cut_x)
3822  return (face_no / 2) ? RefinementCase<dim>::cut_x :
3824  else
3826 }
3827 
3828 
3829 template <>
3830 inline RefinementCase<3>
3832  const RefinementCase<2> &face_refinement_case,
3833  const unsigned int face_no,
3834  const bool face_orientation,
3835  const bool /*face_flip*/,
3836  const bool face_rotation)
3837 {
3838  const unsigned int dim = 3;
3839  AssertIndexRange(face_refinement_case,
3842 
3847 
3848  // correct the face_refinement_case for
3849  // face_orientation and face_rotation. for
3850  // face_orientation, 'true' is the default
3851  // value whereas for face_rotation, 'false'
3852  // is standard. If
3853  // <tt>face_rotation==face_orientation</tt>,
3854  // then one of them is non-standard and we
3855  // have to swap cut_x and cut_y, otherwise no
3856  // change is necessary. face_flip has no
3857  // influence. however, in order to keep the
3858  // interface consistent with other functions,
3859  // we still include it as an argument to this
3860  // function
3861  const RefinementCase<dim - 1> std_face_ref =
3862  (face_orientation == face_rotation) ? flip[face_refinement_case] :
3863  face_refinement_case;
3864 
3865  const RefinementCase<dim> face_to_cell[3][4] = {
3866  {RefinementCase<dim>::no_refinement, // faces 0 and 1
3867  RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3868  RefinementCase<dim>::cut_z,
3870 
3871  {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3872  // "exchanged on faces 2 and 3")
3873  RefinementCase<dim>::cut_z,
3876 
3877  {RefinementCase<dim>::no_refinement, // faces 4 and 5
3881 
3882  return face_to_cell[face_no / 2][std_face_ref];
3883 }
3884 
3885 
3886 
3887 template <int dim>
3888 inline RefinementCase<dim>
3890  const unsigned int)
3891 {
3892  Assert(false, ExcNotImplemented());
3893 
3895 }
3896 
3897 template <>
3898 inline RefinementCase<1>
3900  const unsigned int line_no)
3901 {
3902  (void)line_no;
3903  AssertIndexRange(line_no, 1);
3904 
3905  return RefinementCase<1>::cut_x;
3906 }
3907 
3908 
3909 template <>
3910 inline RefinementCase<2>
3912  const unsigned int line_no)
3913 {
3914  const unsigned int dim = 2;
3915  (void)dim;
3917 
3918  return (line_no / 2) ? RefinementCase<2>::cut_x : RefinementCase<2>::cut_y;
3919 }
3920 
3921 
3922 template <>
3923 inline RefinementCase<3>
3925  const unsigned int line_no)
3926 {
3927  const unsigned int dim = 3;
3929 
3930  const RefinementCase<dim> ref_cases[6] = {
3931  RefinementCase<dim>::cut_y, // lines 0 and 1
3932  RefinementCase<dim>::cut_x, // lines 2 and 3
3933  RefinementCase<dim>::cut_y, // lines 4 and 5
3934  RefinementCase<dim>::cut_x, // lines 6 and 7
3935  RefinementCase<dim>::cut_z, // lines 8 and 9
3936  RefinementCase<dim>::cut_z}; // lines 10 and 11
3937 
3938  return ref_cases[line_no / 2];
3939 }
3940 
3941 
3942 
3943 template <>
3944 inline unsigned int
3945 GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
3946  const bool face_orientation,
3947  const bool face_flip,
3948  const bool face_rotation)
3949 {
3951 
3952  // set up a table to make sure that
3953  // we handle non-standard faces correctly
3954  //
3955  // so set up a table that for each vertex (of
3956  // a quad in standard position) describes
3957  // which vertex to take
3958  //
3959  // first index: four vertices 0...3
3960  //
3961  // second index: face_orientation; 0:
3962  // opposite normal, 1: standard
3963  //
3964  // third index: face_flip; 0: standard, 1:
3965  // face rotated by 180 degrees
3966  //
3967  // forth index: face_rotation: 0: standard,
3968  // 1: face rotated by 90 degrees
3969 
3970  const unsigned int vertex_translation[4][2][2][2] = {
3971  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3972  // face_rotation=false and true
3973  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3974  // face_rotation=false and true
3975  {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
3976  // face_rotation=false and true
3977  {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
3978  // face_rotation=false and true
3979 
3980  {{{2, 3}, // vertex 1 ...
3981  {1, 0}},
3982  {{1, 3}, {2, 0}}},
3983 
3984  {{{1, 0}, // vertex 2 ...
3985  {2, 3}},
3986  {{2, 0}, {1, 3}}},
3987 
3988  {{{3, 1}, // vertex 3 ...
3989  {0, 2}},
3990  {{3, 2}, {0, 1}}}};
3991 
3992  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3993 }
3994 
3995 
3996 
3997 template <int dim>
3998 inline unsigned int
3999 GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
4000  const bool,
4001  const bool,
4002  const bool)
4003 {
4004  Assert(dim > 1, ExcImpossibleInDim(dim));
4006  return vertex;
4007 }
4008 
4009 
4010 
4011 template <>
4012 inline unsigned int
4013 GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
4014  const bool face_orientation,
4015  const bool face_flip,
4016  const bool face_rotation)
4017 {
4019 
4020 
4021  // make sure we handle
4022  // non-standard faces correctly
4023  //
4024  // so set up a table that for each line (of a
4025  // quad) describes which line to take
4026  //
4027  // first index: four lines 0...3
4028  //
4029  // second index: face_orientation; 0:
4030  // opposite normal, 1: standard
4031  //
4032  // third index: face_flip; 0: standard, 1:
4033  // face rotated by 180 degrees
4034  //
4035  // forth index: face_rotation: 0: standard,
4036  // 1: face rotated by 90 degrees
4037 
4038  const unsigned int line_translation[4][2][2][2] = {
4039  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4040  // face_rotation=false and true
4041  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4042  // face_rotation=false and true
4043  {{0, 3}, // line 0, face_orientation=true, face_flip=false,
4044  // face_rotation=false and true
4045  {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
4046  // face_rotation=false and true
4047 
4048  {{{3, 1}, // line 1 ...
4049  {2, 0}},
4050  {{1, 2}, {0, 3}}},
4051 
4052  {{{0, 3}, // line 2 ...
4053  {1, 2}},
4054  {{2, 0}, {3, 1}}},
4055 
4056  {{{1, 2}, // line 3 ...
4057  {0, 3}},
4058  {{3, 1}, {2, 0}}}};
4059 
4060  return line_translation[line][face_orientation][face_flip][face_rotation];
4061 }
4062 
4063 
4064 
4065 template <int dim>
4066 inline unsigned int
4067 GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
4068  const bool,
4069  const bool,
4070  const bool)
4071 {
4072  Assert(false, ExcNotImplemented());
4073  return line;
4074 }
4075 
4076 
4077 
4078 template <>
4079 inline unsigned int
4080 GeometryInfo<2>::standard_to_real_line_vertex(const unsigned int vertex,
4081  const bool line_orientation)
4082 {
4083  return line_orientation ? vertex : (1 - vertex);
4084 }
4085 
4086 
4087 
4088 template <int dim>
4089 inline unsigned int
4090 GeometryInfo<dim>::standard_to_real_line_vertex(const unsigned int vertex,
4091  const bool)
4092 {
4093  Assert(false, ExcNotImplemented());
4094  return vertex;
4095 }
4096 
4097 
4098 
4099 template <>
4100 inline std::array<unsigned int, 2>
4102  const unsigned int vertex)
4103 {
4104  return {{vertex % 2, vertex / 2}};
4105 }
4106 
4107 
4108 
4109 template <int dim>
4110 inline std::array<unsigned int, 2>
4112  const unsigned int vertex)
4113 {
4114  Assert(false, ExcNotImplemented());
4115  (void)vertex;
4116  return {{0, 0}};
4117 }
4118 
4119 
4120 
4121 template <>
4122 inline std::array<unsigned int, 2>
4124 {
4125  // set up a table that for each
4126  // line describes a) from which
4127  // quad to take it, b) which line
4128  // therein it is if the face is
4129  // oriented correctly
4130  static const unsigned int lookup_table[GeometryInfo<3>::lines_per_cell][2] = {
4131  {4, 0}, // take first four lines from bottom face
4132  {4, 1},
4133  {4, 2},
4134  {4, 3},
4135 
4136  {5, 0}, // second four lines from top face
4137  {5, 1},
4138  {5, 2},
4139  {5, 3},
4140 
4141  {0, 0}, // the rest randomly
4142  {1, 0},
4143  {0, 1},
4144  {1, 1}};
4145 
4146  return {{lookup_table[i][0], lookup_table[i][1]}};
4147 }
4148 
4149 
4150 
4151 template <int dim>
4152 inline std::array<unsigned int, 2>
4154 {
4155  Assert(false, ExcNotImplemented());
4156  (void)line;
4157  return {{0, 0}};
4158 }
4159 
4160 
4161 
4162 template <>
4163 inline std::array<unsigned int, 2>
4165  const unsigned int vertex)
4166 {
4167  // get the corner indices by asking either the bottom or the top face for its
4168  // vertices. handle non-standard faces by calling the vertex reordering
4169  // function from GeometryInfo
4170 
4171  // bottom face (4) for first four vertices, top face (5) for the rest
4172  return {{4 + vertex / 4, vertex % 4}};
4173 }
4174 
4175 
4176 
4177 template <int dim>
4178 inline std::array<unsigned int, 2>
4180  const unsigned int vertex)
4181 {
4182  Assert(false, ExcNotImplemented());
4183  (void)vertex;
4184  return {{0, 0}};
4185 }
4186 
4187 
4188 
4189 template <>
4190 inline unsigned int
4191 GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
4192  const bool face_orientation,
4193  const bool face_flip,
4194  const bool face_rotation)
4195 {
4197 
4198 
4199  // make sure we handle
4200  // non-standard faces correctly
4201  //
4202  // so set up a table that for each line (of a
4203  // quad) describes which line to take
4204  //
4205  // first index: four lines 0...3
4206  //
4207  // second index: face_orientation; 0:
4208  // opposite normal, 1: standard
4209  //
4210  // third index: face_flip; 0: standard, 1:
4211  // face rotated by 180 degrees
4212  //
4213  // forth index: face_rotation: 0: standard,
4214  // 1: face rotated by 90 degrees
4215 
4216  const unsigned int line_translation[4][2][2][2] = {
4217  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4218  // face_rotation=false and true
4219  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4220  // face_rotation=false and true
4221  {{0, 2}, // line 0, face_orientation=true, face_flip=false,
4222  // face_rotation=false and true
4223  {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
4224  // face_rotation=false and true
4225 
4226  {{{3, 1}, // line 1 ...
4227  {2, 0}},
4228  {{1, 3}, {0, 2}}},
4229 
4230  {{{0, 3}, // line 2 ...
4231  {1, 2}},
4232  {{2, 1}, {3, 0}}},
4233 
4234  {{{1, 2}, // line 3 ...
4235  {0, 3}},
4236  {{3, 0}, {2, 1}}}};
4237 
4238  return line_translation[line][face_orientation][face_flip][face_rotation];
4239 }
4240 
4241 
4242 
4243 template <int dim>
4244 inline unsigned int
4245 GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
4246  const bool,
4247  const bool,
4248  const bool)
4249 {
4250  Assert(false, ExcNotImplemented());
4251  return line;
4252 }
4253 
4254 
4255 
4256 template <>
4257 inline unsigned int
4259  const unsigned int face,
4260  const unsigned int subface,
4261  const bool,
4262  const bool,
4263  const bool,
4264  const RefinementCase<0> &)
4265 {
4266  (void)subface;
4269 
4270  return face;
4271 }
4272 
4273 
4274 
4275 template <>
4276 inline unsigned int
4278  const unsigned int face,
4279  const unsigned int subface,
4280  const bool /*face_orientation*/,
4281  const bool face_flip,
4282  const bool /*face_rotation*/,
4283  const RefinementCase<1> &)
4284 {
4287 
4288  // always return the child adjacent to the specified
4289  // subface. if the face of a cell is not refined, don't
4290  // throw an assertion but deliver the child adjacent to
4291  // the face nevertheless, i.e. deliver the child of
4292  // this cell adjacent to the subface of a possibly
4293  // refined neighbor. this simplifies setting neighbor
4294  // information in execute_refinement.
4295  const unsigned int
4297  [max_children_per_face] = {
4298  {
4299  // Normal orientation (face_flip = false)
4300  {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
4301  {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
4302  {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
4303  },
4304  {
4305  // Flipped orientation (face_flip = true)
4306  {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
4307  {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
4308  {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
4309  }};
4310 
4311  return subcells[face_flip][ref_case - 1][face][subface];
4312 }
4313 
4314 
4315 
4316 template <>
4317 inline unsigned int
4319  const unsigned int face,
4320  const unsigned int subface,
4321  const bool face_orientation,
4322  const bool face_flip,
4323  const bool face_rotation,
4324  const RefinementCase<2> &face_ref_case)
4325 {
4326  const unsigned int dim = 3;
4327 
4329  ExcMessage("Cell has no children."));
4331  if (!(subface == 0 &&
4332  face_ref_case == RefinementCase<dim - 1>::no_refinement))
4333  {
4334  AssertIndexRange(subface,
4335  GeometryInfo<dim - 1>::n_children(face_ref_case));
4336  }
4337 
4338  // invalid number used for invalid cases,
4339  // e.g. when the children are more refined at
4340  // a given face than the face itself
4341  const unsigned int e = numbers::invalid_unsigned_int;
4342 
4343  // the whole process of finding a child cell
4344  // at a given subface considering the
4345  // possibly anisotropic refinement cases of
4346  // the cell and the face as well as
4347  // orientation, flip and rotation of the face
4348  // is quite complicated. thus, we break it
4349  // down into several steps.
4350 
4351  // first step: convert the given face refine
4352  // case to a face refine case concerning the
4353  // face in standard orientation (, flip and
4354  // rotation). This only affects cut_x and
4355  // cut_y
4356  const RefinementCase<dim - 1> flip[4] = {
4357  RefinementCase<dim - 1>::no_refinement,
4358  RefinementCase<dim - 1>::cut_y,
4359  RefinementCase<dim - 1>::cut_x,
4360  RefinementCase<dim - 1>::cut_xy};
4361  // for face_orientation, 'true' is the
4362  // default value whereas for face_rotation,
4363  // 'false' is standard. If
4364  // <tt>face_rotation==face_orientation</tt>,
4365  // then one of them is non-standard and we
4366  // have to swap cut_x and cut_y, otherwise no
4367  // change is necessary.
4368  const RefinementCase<dim - 1> std_face_ref =
4369  (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4370 
4371  // second step: convert the given subface
4372  // index to the one for a standard face
4373  // respecting face_orientation, face_flip and
4374  // face_rotation
4375 
4376  // first index: face_ref_case
4377  // second index: face_orientation
4378  // third index: face_flip
4379  // forth index: face_rotation
4380  // fifth index: subface index
4381  const unsigned int subface_exchange[4][2][2][2][4] = {
4382  // no_refinement (subface 0 stays 0,
4383  // all others are invalid)
4384  {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4385  {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4386  // cut_x (here, if the face is only
4387  // rotated OR only falsely oriented,
4388  // then subface 0 of the non-standard
4389  // face does NOT correspond to one of
4390  // the subfaces of a standard
4391  // face. Thus we indicate the subface
4392  // which is located at the lower left
4393  // corner (the origin of the face's
4394  // local coordinate system) with
4395  // '0'. The rest of this issue is
4396  // taken care of using the above
4397  // conversion to a 'standard face
4398  // refine case')
4399  {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4400  {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4401  // cut_y (the same applies as for
4402  // cut_x)
4403  {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4404  {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4405  // cut_xyz: this information is
4406  // identical to the information
4407  // returned by
4408  // GeometryInfo<3>::real_to_standard_face_vertex()
4409  {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4410  // face_rotation=false, subfaces 0,1,2,3
4411  {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4412  // face_rotation=true, subfaces 0,1,2,3
4413  {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4414  // face_rotation=false, subfaces 0,1,2,3
4415  {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4416  // face_rotation=true, subfaces 0,1,2,3
4417  {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4418  // face_rotation=false, subfaces 0,1,2,3
4419  {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4420  // face_rotation=true, subfaces 0,1,2,3
4421  {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4422  // face_rotation=false, subfaces 0,1,2,3
4423  {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4424  // face_rotation=true, subfaces 0,1,2,3
4425 
4426  const unsigned int std_subface =
4427  subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4428  [subface];
4429  Assert(std_subface != e, ExcInternalError());
4430 
4431  // third step: these are the children, which
4432  // can be found at the given subfaces of an
4433  // isotropically refined (standard) face
4434  //
4435  // first index: (refinement_case-1)
4436  // second index: face_index
4437  // third index: subface_index (isotropic refinement)
4438  const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4439  [max_children_per_face] = {
4440  // cut_x
4441  {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4442  {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4443  {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4444  {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4445  {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4446  {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4447  // cut_y
4448  {{0, 1, 0, 1},
4449  {0, 1, 0, 1},
4450  {0, 0, 0, 0},
4451  {1, 1, 1, 1},
4452  {0, 0, 1, 1},
4453  {0, 0, 1, 1}},
4454  // cut_xy
4455  {{0, 2, 0, 2},
4456  {1, 3, 1, 3},
4457  {0, 0, 1, 1},
4458  {2, 2, 3, 3},
4459  {0, 1, 2, 3},
4460  {0, 1, 2, 3}},
4461  // cut_z
4462  {{0, 0, 1, 1},
4463  {0, 0, 1, 1},
4464  {0, 1, 0, 1},
4465  {0, 1, 0, 1},
4466  {0, 0, 0, 0},
4467  {1, 1, 1, 1}},
4468  // cut_xz
4469  {{0, 0, 1, 1},
4470  {2, 2, 3, 3},
4471  {0, 1, 2, 3},
4472  {0, 1, 2, 3},
4473  {0, 2, 0, 2},
4474  {1, 3, 1, 3}},
4475  // cut_yz
4476  {{0, 1, 2, 3},
4477  {0, 1, 2, 3},
4478  {0, 2, 0, 2},
4479  {1, 3, 1, 3},
4480  {0, 0, 1, 1},
4481  {2, 2, 3, 3}},
4482  // cut_xyz
4483  {{0, 2, 4, 6},
4484  {1, 3, 5, 7},
4485  {0, 4, 1, 5},
4486  {2, 6, 3, 7},
4487  {0, 1, 2, 3},
4488  {4, 5, 6, 7}}};
4489 
4490  // forth step: check, whether the given face
4491  // refine case is valid for the given cell
4492  // refine case. this is the case, if the
4493  // given face refine case is at least as
4494  // refined as the face is for the given cell
4495  // refine case
4496 
4497  // note, that we are considering standard
4498  // face refinement cases here and thus must
4499  // not pass the given orientation, flip and
4500  // rotation flags
4501  if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4502  face_refinement_case(ref_case, face))
4503  {
4504  // all is fine. for anisotropic face
4505  // refine cases, select one of the
4506  // isotropic subfaces which neighbors the
4507  // same child
4508 
4509  // first index: (standard) face refine case
4510  // second index: subface index
4511  const unsigned int equivalent_iso_subface[4][4] = {
4512  {0, e, e, e}, // no_refinement
4513  {0, 3, e, e}, // cut_x
4514  {0, 3, e, e}, // cut_y
4515  {0, 1, 2, 3}}; // cut_xy
4516 
4517  const unsigned int equ_std_subface =
4518  equivalent_iso_subface[std_face_ref][std_subface];
4519  Assert(equ_std_subface != e, ExcInternalError());
4520 
4521  return iso_children[ref_case - 1][face][equ_std_subface];
4522  }
4523  else
4524  {
4525  // the face_ref_case was too coarse,
4526  // throw an error
4527  Assert(false,
4528  ExcMessage("The face RefineCase is too coarse "
4529  "for the given cell RefineCase."));
4530  }
4531  // we only get here in case of an error
4532  return e;
4533 }
4534 
4535 
4536 
4537 template <>
4538 inline unsigned int
4540  const unsigned int,
4541  const unsigned int,
4542  const bool,
4543  const bool,
4544  const bool,
4545  const RefinementCase<3> &)
4546 {
4547  Assert(false, ExcNotImplemented());
4549 }
4550 
4551 
4552 
4553 template <>
4554 inline unsigned int
4555 GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4556  const unsigned int line,
4557  const bool,
4558  const bool,
4559  const bool)
4560 {
4561  (void)line;
4564 
4565  // The face is a line itself.
4566  return face;
4567 }
4568 
4569 
4570 
4571 template <>
4572 inline unsigned int
4573 GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4574  const unsigned int line,
4575  const bool face_orientation,
4576  const bool face_flip,
4577  const bool face_rotation)
4578 {
4581 
4582  const unsigned lines[faces_per_cell][lines_per_face] = {
4583  {8, 10, 0, 4}, // left face
4584  {9, 11, 1, 5}, // right face
4585  {2, 6, 8, 9}, // front face
4586  {3, 7, 10, 11}, // back face
4587  {0, 1, 2, 3}, // bottom face
4588  {4, 5, 6, 7}}; // top face
4589  return lines[face][real_to_standard_face_line(
4590  line, face_orientation, face_flip, face_rotation)];
4591 }
4592 
4593 
4594 
4595 template <int dim>
4596 inline unsigned int
4597 GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4598  const unsigned int,
4599  const bool,
4600  const bool,
4601  const bool)
4602 {
4603  Assert(false, ExcNotImplemented());
4605 }
4606 
4607 
4608 
4609 template <int dim>
4610 inline unsigned int
4611 GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4612  const unsigned int vertex,
4613  const bool face_orientation,
4614  const bool face_flip,
4615  const bool face_rotation)
4616 {
4618  face,
4619  vertex,
4620  face_orientation,
4621  face_flip,
4622  face_rotation);
4623 }
4624 
4625 
4626 
4627 template <int dim>
4628 inline Point<dim>
4630 {
4631  Point<dim> p = q;
4632  for (unsigned int i = 0; i < dim; i++)
4633  if (p[i] < 0.)
4634  p[i] = 0.;
4635  else if (p[i] > 1.)
4636  p[i] = 1.;
4637 
4638  return p;
4639 }
4640 
4641 
4642 
4643 template <int dim>
4644 inline double
4646 {
4647  double result = 0.0;
4648 
4649  for (unsigned int i = 0; i < dim; i++)
4650  if ((-p[i]) > result)
4651  result = -p[i];
4652  else if ((p[i] - 1.) > result)
4653  result = (p[i] - 1.);
4654 
4655  return result;
4656 }
4657 
4658 
4659 
4660 template <int dim>
4661 inline double
4663  const unsigned int i)
4664 {
4666 
4667  switch (dim)
4668  {
4669  case 1:
4670  {
4671  const double x = xi[0];
4672  switch (i)
4673  {
4674  case 0:
4675  return 1 - x;
4676  case 1:
4677  return x;
4678  }
4679  break;
4680  }
4681 
4682  case 2:
4683  {
4684  const double x = xi[0];
4685  const double y = xi[1];
4686  switch (i)
4687  {
4688  case 0:
4689  return (1 - x) * (1 - y);
4690  case 1:
4691  return x * (1 - y);
4692  case 2:
4693  return (1 - x) * y;
4694  case 3:
4695  return x * y;
4696  }
4697  break;
4698  }
4699 
4700  case 3:
4701  {
4702  const double x = xi[0];
4703  const double y = xi[1];
4704  const double z = xi[2];
4705  switch (i)
4706  {
4707  case 0:
4708  return (1 - x) * (1 - y) * (1 - z);
4709  case 1:
4710  return x * (1 - y) * (1 - z);
4711  case 2:
4712  return (1 - x) * y * (1 - z);
4713  case 3:
4714  return x * y * (1 - z);
4715  case 4:
4716  return (1 - x) * (1 - y) * z;
4717  case 5:
4718  return x * (1 - y) * z;
4719  case 6:
4720  return (1 - x) * y * z;
4721  case 7:
4722  return x * y * z;
4723  }
4724  break;
4725  }
4726 
4727  default:
4728  Assert(false, ExcNotImplemented());
4729  }
4730  return -1e9;
4731 }
4732 
4733 
4734 
4735 template <>
4737  const Point<1> &,
4738  const unsigned int i)
4739 {
4741 
4742  switch (i)
4743  {
4744  case 0:
4745  return Point<1>(-1.);
4746  case 1:
4747  return Point<1>(1.);
4748  }
4749 
4750  return Point<1>(-1e9);
4751 }
4752 
4753 
4754 
4755 template <>
4757  const Point<2> & xi,
4758  const unsigned int i)
4759 {
4761 
4762  const double x = xi[0];
4763  const double y = xi[1];
4764  switch (i)
4765  {
4766  case 0:
4767  return Point<2>(-(1 - y), -(1 - x));
4768  case 1:
4769  return Point<2>(1 - y, -x);
4770  case 2:
4771  return Point<2>(-y, 1 - x);
4772  case 3:
4773  return Point<2>(y, x);
4774  }
4775  return Point<2>(-1e9, -1e9);
4776 }
4777 
4778 
4779 
4780 template <>
4782  const Point<3> & xi,
4783  const unsigned int i)
4784 {
4786 
4787  const double x = xi[0];
4788  const double y = xi[1];
4789  const double z = xi[2];
4790  switch (i)
4791  {
4792  case 0:
4793  return Point<3>(-(1 - y) * (1 - z),
4794  -(1 - x) * (1 - z),
4795  -(1 - x) * (1 - y));
4796  case 1:
4797  return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4798  case 2:
4799  return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4800  case 3:
4801  return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4802  case 4:
4803  return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4804  case 5:
4805  return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4806  case 6:
4807  return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4808  case 7:
4809  return Point<3>(y * z, x * z, x * y);
4810  }
4811 
4812  return Point<3>(-1e9, -1e9, -1e9);
4813 }
4814 
4815 
4816 
4817 template <int dim>
4818 inline Tensor<1, dim>
4820  const unsigned int)
4821 {
4822  Assert(false, ExcNotImplemented());
4823  return Tensor<1, dim>();
4824 }
4825 
4826 
4827 unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4828 {
4829  return 0;
4830 }
4831 
4832 
4833 namespace internal
4834 {
4835  namespace GeometryInfoHelper
4836  {
4837  // wedge product of a single
4838  // vector in 2d: we just have to
4839  // rotate it by 90 degrees to the
4840  // right
4841  inline Tensor<1, 2>
4842  wedge_product(const Tensor<1, 2> (&derivative)[1])
4843  {
4844  Tensor<1, 2> result;
4845  result[0] = derivative[0][1];
4846  result[1] = -derivative[0][0];
4847 
4848  return result;
4849  }
4850 
4851 
4852  // wedge product of 2 vectors in
4853  // 3d is the cross product
4854  inline Tensor<1, 3>
4855  wedge_product(const Tensor<1, 3> (&derivative)[2])
4856  {
4857  return cross_product_3d(derivative[0], derivative[1]);
4858  }
4859 
4860 
4861  // wedge product of dim vectors
4862  // in dim-d: that's the
4863  // determinant of the matrix
4864  template <int dim>
4865  inline Tensor<0, dim>
4866  wedge_product(const Tensor<1, dim> (&derivative)[dim])
4867  {
4868  Tensor<2, dim> jacobian;
4869  for (unsigned int i = 0; i < dim; ++i)
4870  jacobian[i] = derivative[i];
4871 
4872  return determinant(jacobian);
4873  }
4874  } // namespace GeometryInfoHelper
4875 } // namespace internal
4876 
4877 
4878 template <int dim>
4879 template <int spacedim>
4880 inline void
4882 # ifndef DEAL_II_CXX14_CONSTEXPR_BUG
4885 # else
4887 # endif
4888 {
4889  // for each of the vertices,
4890  // compute the alternating form
4891  // of the mapped unit
4892  // vectors. consider for
4893  // example the case of a quad
4894  // in spacedim==3: to do so, we
4895  // need to see how the
4896  // infinitesimal vectors
4897  // (d\xi_1,0) and (0,d\xi_2)
4898  // are transformed into
4899  // spacedim-dimensional space
4900  // and then form their cross
4901  // product (i.e. the wedge product
4902  // of two vectors). to this end, note
4903  // that
4904  // \vec x = sum_i \vec v_i phi_i(\vec xi)
4905  // so the transformed vectors are
4906  // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
4907  // and
4908  // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
4909  // which boils down to the columns
4910  // of the 3x2 matrix \grad_\xi x(\xi)
4911  //
4912  // a similar reasoning would
4913  // hold for all dim,spacedim
4914  // pairs -- we only have to
4915  // compute the wedge product of
4916  // the columns of the
4917  // derivatives
4918  for (unsigned int i = 0; i < vertices_per_cell; ++i)
4919  {
4920  Tensor<1, spacedim> derivatives[dim];
4921 
4922  for (unsigned int j = 0; j < vertices_per_cell; ++j)
4923  {
4924  const Tensor<1, dim> grad_phi_j =
4926  for (unsigned int l = 0; l < dim; ++l)
4927  derivatives[l] += vertices[j] * grad_phi_j[l];
4928  }
4929 
4930  forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
4931  }
4932 }
4933 
4934 #endif // DOXYGEN
4936 
4937 #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:191
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)
std::vector< unsigned int > vertex_indices
Definition: tria.cc:2244
void serialize(Archive &ar, const unsigned int version)
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1628
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:2430
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:1403
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:363
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 Point< dim > project_to_unit_cell(const Point< dim > &p)
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:448
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:362
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)
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()