Reference documentation for deal.II version Git 9b341ae816 2021-01-21 15:59:19 -0700
\(\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 
866  template <class Archive>
867  void
868  serialize(Archive &ar, const unsigned int version);
869 
874  ExcInvalidRefinementCase,
875  int,
876  << "The refinement flags given (" << arg1
877  << ") contain set bits that do not "
878  << "make sense for the space dimension of the object to which they are applied.");
879 
880 private:
885  std::uint8_t value : (dim > 0 ? dim : 1);
886 };
887 
888 
889 namespace internal
890 {
910  template <int dim>
912  {
917  {
921  case_none = 0,
922 
926  case_isotropic = static_cast<std::uint8_t>(-1)
927  };
928  };
929 
930 
939  template <>
941  {
948  {
952  case_none = 0,
953 
957  case_isotropic = case_none
958  };
959  };
960 
961 
962 
972  template <>
974  {
981  {
985  case_none = 0,
986 
990  case_isotropic = case_none
991  };
992  };
993 
994 
995 
1006  template <>
1008  {
1016  {
1020  case_none = 0,
1024  case_x = 1,
1028  case_isotropic = case_x
1029  };
1030  };
1031 
1032 
1033 
1125  template <>
1127  {
1135  {
1136  case_none = 0,
1137  case_x = 1,
1138  case_x1y = 2,
1139  case_x2y = 3,
1140  case_x1y2y = 4,
1141  case_y = 5,
1142  case_y1x = 6,
1143  case_y2x = 7,
1144  case_y1x2x = 8,
1145  case_xy = 9,
1146 
1147  case_isotropic = case_xy
1148  };
1149  };
1150 
1151 
1152 
1159  template <int dim>
1160  class SubfaceCase : public SubfacePossibilities<dim>
1161  {
1162  public:
1169  subface_possibility);
1170 
1182  operator std::uint8_t() const;
1183 
1187  static constexpr std::size_t
1189 
1194  ExcInvalidSubfaceCase,
1195  int,
1196  << "The subface case given (" << arg1 << ") does not make sense "
1197  << "for the space dimension of the object to which they are applied.");
1198 
1199  private:
1204  std::uint8_t value : (dim == 3 ? 4 : 1);
1205  };
1206 
1207 } // namespace internal
1208 
1209 
1210 
1211 template <int dim>
1213 
1214 
1215 
1232 template <>
1233 struct GeometryInfo<0>
1234 {
1242  static constexpr unsigned int max_children_per_cell = 1;
1243 
1247  static constexpr unsigned int faces_per_cell = 0;
1248 
1265  static std::array<unsigned int, 0>
1266  face_indices();
1267 
1275  static constexpr unsigned int max_children_per_face = 0;
1276 
1282  static unsigned int
1283  n_children(const RefinementCase<0> &refinement_case);
1284 
1288  static constexpr unsigned int vertices_per_cell = 1;
1289 
1308  static std::array<unsigned int, vertices_per_cell>
1309  vertex_indices();
1310 
1334  static unsigned int
1335  face_to_cell_vertices(const unsigned int face,
1336  const unsigned int vertex,
1337  const bool face_orientation = true,
1338  const bool face_flip = false,
1339  const bool face_rotation = false);
1340 
1355  static unsigned int
1356  face_to_cell_lines(const unsigned int face,
1357  const unsigned int line,
1358  const bool face_orientation = true,
1359  const bool face_flip = false,
1360  const bool face_rotation = false);
1361 
1368  static constexpr unsigned int vertices_per_face = 0;
1369 
1373  static constexpr unsigned int lines_per_face = 0;
1374 
1378  static constexpr unsigned int quads_per_face = 0;
1379 
1383  static constexpr unsigned int lines_per_cell = 0;
1384 
1388  static constexpr unsigned int quads_per_cell = 0;
1389 
1393  static constexpr unsigned int hexes_per_cell = 0;
1394 
1412  static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1413 
1427  static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1428 };
1429 
1430 
1431 
1955 template <int dim>
1956 struct GeometryInfo
1957 {
1965  static constexpr unsigned int max_children_per_cell = 1 << dim;
1966 
1970  static constexpr unsigned int faces_per_cell = 2 * dim;
1971 
1989  face_indices();
1990 
1998  static constexpr unsigned int max_children_per_face =
1999  GeometryInfo<dim - 1>::max_children_per_cell;
2000 
2004  static constexpr unsigned int vertices_per_cell = 1 << dim;
2005 
2022  vertex_indices();
2023 
2027  static constexpr unsigned int vertices_per_face =
2028  GeometryInfo<dim - 1>::vertices_per_cell;
2029 
2033  static constexpr unsigned int lines_per_face =
2034  GeometryInfo<dim - 1>::lines_per_cell;
2035 
2039  static constexpr unsigned int quads_per_face =
2040  GeometryInfo<dim - 1>::quads_per_cell;
2041 
2051  static constexpr unsigned int lines_per_cell =
2052  (2 * GeometryInfo<dim - 1>::lines_per_cell +
2053  GeometryInfo<dim - 1>::vertices_per_cell);
2054 
2062  static constexpr unsigned int quads_per_cell =
2063  (2 * GeometryInfo<dim - 1>::quads_per_cell +
2064  GeometryInfo<dim - 1>::lines_per_cell);
2065 
2069  static constexpr unsigned int hexes_per_cell =
2070  (2 * GeometryInfo<dim - 1>::hexes_per_cell +
2071  GeometryInfo<dim - 1>::quads_per_cell);
2072 
2090  static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
2091  internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal();
2092 
2106  static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
2107  internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal();
2108 
2119  static constexpr std::array<std::array<unsigned int, dim>, vertices_per_cell>
2120  vertex_to_face =
2121  internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face();
2122 
2127  static unsigned int
2128  n_children(const RefinementCase<dim> &refinement_case);
2129 
2134  static unsigned int
2135  n_subfaces(const internal::SubfaceCase<dim> &subface_case);
2136 
2146  static double
2147  subface_ratio(const internal::SubfaceCase<dim> &subface_case,
2148  const unsigned int subface_no);
2149 
2155  static RefinementCase<dim - 1>
2156  face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2157  const unsigned int face_no,
2158  const bool face_orientation = true,
2159  const bool face_flip = false,
2160  const bool face_rotation = false);
2161 
2167  static RefinementCase<dim>
2168  min_cell_refinement_case_for_face_refinement(
2169  const RefinementCase<dim - 1> &face_refinement_case,
2170  const unsigned int face_no,
2171  const bool face_orientation = true,
2172  const bool face_flip = false,
2173  const bool face_rotation = false);
2174 
2179  static RefinementCase<1>
2180  line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2181  const unsigned int line_no);
2182 
2187  static RefinementCase<dim>
2188  min_cell_refinement_case_for_line_refinement(const unsigned int line_no);
2189 
2236  static unsigned int
2237  child_cell_on_face(const RefinementCase<dim> & ref_case,
2238  const unsigned int face,
2239  const unsigned int subface,
2240  const bool face_orientation = true,
2241  const bool face_flip = false,
2242  const bool face_rotation = false,
2243  const RefinementCase<dim - 1> &face_refinement_case =
2245 
2259  static unsigned int
2260  line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
2261 
2282  static unsigned int
2283  face_to_cell_vertices(const unsigned int face,
2284  const unsigned int vertex,
2285  const bool face_orientation = true,
2286  const bool face_flip = false,
2287  const bool face_rotation = false);
2288 
2300  static unsigned int
2301  face_to_cell_lines(const unsigned int face,
2302  const unsigned int line,
2303  const bool face_orientation = true,
2304  const bool face_flip = false,
2305  const bool face_rotation = false);
2306 
2316  static unsigned int
2317  standard_to_real_face_vertex(const unsigned int vertex,
2318  const bool face_orientation = true,
2319  const bool face_flip = false,
2320  const bool face_rotation = false);
2321 
2331  static unsigned int
2332  real_to_standard_face_vertex(const unsigned int vertex,
2333  const bool face_orientation = true,
2334  const bool face_flip = false,
2335  const bool face_rotation = false);
2336 
2346  static unsigned int
2347  standard_to_real_face_line(const unsigned int line,
2348  const bool face_orientation = true,
2349  const bool face_flip = false,
2350  const bool face_rotation = false);
2351 
2357  static unsigned int
2358  standard_to_real_line_vertex(const unsigned int vertex,
2359  const bool line_orientation = true);
2360 
2368  static std::array<unsigned int, 2>
2369  standard_quad_vertex_to_line_vertex_index(const unsigned int vertex);
2370 
2378  static std::array<unsigned int, 2>
2379  standard_hex_vertex_to_quad_vertex_index(const unsigned int vertex);
2380 
2388  static std::array<unsigned int, 2>
2389  standard_hex_line_to_quad_line_index(const unsigned int line);
2390 
2400  static unsigned int
2401  real_to_standard_face_line(const unsigned int line,
2402  const bool face_orientation = true,
2403  const bool face_flip = false,
2404  const bool face_rotation = false);
2405 
2411  static Point<dim>
2412  unit_cell_vertex(const unsigned int vertex);
2413 
2423  static unsigned int
2424  child_cell_from_point(const Point<dim> &p);
2425 
2433  static Point<dim>
2434  cell_to_child_coordinates(const Point<dim> & p,
2435  const unsigned int child_index,
2436  const RefinementCase<dim> refine_case =
2438 
2444  static Point<dim>
2445  child_to_cell_coordinates(const Point<dim> & p,
2446  const unsigned int child_index,
2447  const RefinementCase<dim> refine_case =
2449 
2454  static bool
2455  is_inside_unit_cell(const Point<dim> &p);
2456 
2468  static bool
2469  is_inside_unit_cell(const Point<dim> &p, const double eps);
2470 
2475  template <typename Number = double>
2476  static Point<dim, Number>
2477  project_to_unit_cell(const Point<dim, Number> &p);
2478 
2484  static double
2485  distance_to_unit_cell(const Point<dim> &p);
2486 
2491  static double
2492  d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2493 
2498  static Tensor<1, dim>
2499  d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2500 
2552  template <int spacedim>
2553  static void
2554  alternating_form_at_vertices
2555 #ifndef DEAL_II_CXX14_CONSTEXPR_BUG
2556  (const Point<spacedim> (&vertices)[vertices_per_cell],
2557  Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2558 #else
2559  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2560 #endif
2561 
2571  static constexpr std::array<unsigned int, faces_per_cell>
2572  unit_normal_direction =
2573  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction();
2574 
2591  static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2592  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation();
2593 
2604  static constexpr std::array<Tensor<1, dim>, faces_per_cell>
2605  unit_normal_vector =
2606  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_vector();
2607 
2621  static constexpr std::array<std::array<Tensor<1, dim>, dim - 1>,
2622  faces_per_cell>
2623  unit_tangential_vectors = internal::GeometryInfoHelper::Initializers<
2624  dim>::unit_tangential_vectors();
2625 
2631  static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2632  internal::GeometryInfoHelper::Initializers<dim>::opposite_face();
2633 
2634 
2638  DeclException1(ExcInvalidCoordinate,
2639  double,
2640  << "The coordinates must satisfy 0 <= x_i <= 1, "
2641  << "but here we have x_i=" << arg1);
2642 
2646  DeclException3(ExcInvalidSubface,
2647  int,
2648  int,
2649  int,
2650  << "RefinementCase<dim> " << arg1 << ": face " << arg2
2651  << " has no subface " << arg3);
2652 };
2653 
2654 
2655 
2656 #ifndef DOXYGEN
2657 
2658 
2659 /* -------------- declaration of explicit specializations ------------- */
2660 
2661 
2662 template <>
2665  const unsigned int i);
2666 template <>
2669  const unsigned int i);
2670 template <>
2673  const unsigned int i);
2674 
2675 
2676 
2677 /* -------------- inline functions ------------- */
2678 
2679 
2680 inline GeometryPrimitive::GeometryPrimitive(const Object object)
2681  : object(object)
2682 {}
2683 
2684 
2685 
2686 inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2687  : object(static_cast<Object>(object_dimension))
2688 {}
2689 
2690 
2691 inline GeometryPrimitive::operator unsigned int() const
2692 {
2693  return static_cast<unsigned int>(object);
2694 }
2695 
2696 
2697 
2698 namespace internal
2699 {
2700  template <int dim>
2701  inline SubfaceCase<dim>::SubfaceCase(
2702  const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2703  : value(subface_possibility)
2704  {}
2705 
2706 
2707  template <int dim>
2708  inline SubfaceCase<dim>::operator std::uint8_t() const
2709  {
2710  return value;
2711  }
2712 
2713 
2714 } // namespace internal
2715 
2716 
2717 template <int dim>
2718 inline RefinementCase<dim>
2719 RefinementCase<dim>::cut_axis(const unsigned int)
2720 {
2721  Assert(false, ExcInternalError());
2722  return static_cast<std::uint8_t>(-1);
2723 }
2724 
2725 
2726 template <>
2727 inline RefinementCase<1>
2728 RefinementCase<1>::cut_axis(const unsigned int i)
2729 {
2730  AssertIndexRange(i, 1);
2731 
2733  return options[i];
2734 }
2735 
2736 
2737 
2738 template <>
2739 inline RefinementCase<2>
2740 RefinementCase<2>::cut_axis(const unsigned int i)
2741 {
2742  AssertIndexRange(i, 2);
2743 
2746  return options[i];
2747 }
2748 
2749 
2750 
2751 template <>
2752 inline RefinementCase<3>
2753 RefinementCase<3>::cut_axis(const unsigned int i)
2754 {
2755  AssertIndexRange(i, 3);
2756 
2760  return options[i];
2761 }
2762 
2763 
2764 
2765 template <int dim>
2767  : value(RefinementPossibilities<dim>::no_refinement)
2768 {}
2769 
2770 
2771 
2772 template <int dim>
2774  const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2775  : value(refinement_case)
2776 {
2777  // check that only those bits of
2778  // the given argument are set that
2779  // make sense for a given space
2780  // dimension
2781  Assert((refinement_case &
2783  refinement_case,
2784  ExcInvalidRefinementCase(refinement_case));
2785 }
2786 
2787 
2788 
2789 template <int dim>
2790 inline RefinementCase<dim>::RefinementCase(const std::uint8_t refinement_case)
2791  : value(refinement_case)
2792 {
2793  // check that only those bits of
2794  // the given argument are set that
2795  // make sense for a given space
2796  // dimension
2797  Assert((refinement_case &
2799  refinement_case,
2800  ExcInvalidRefinementCase(refinement_case));
2801 }
2802 
2803 
2804 
2805 template <int dim>
2806 inline RefinementCase<dim>::operator std::uint8_t() const
2807 {
2808  return value;
2809 }
2810 
2811 
2812 
2813 template <int dim>
2814 inline RefinementCase<dim>
2816 {
2817  return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2818 }
2819 
2820 
2821 
2822 template <int dim>
2824  operator&(const RefinementCase<dim> &r) const
2825 {
2826  return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2827 }
2828 
2829 
2830 
2831 template <int dim>
2832 inline RefinementCase<dim>
2834 {
2835  return RefinementCase<dim>(static_cast<std::uint8_t>(
2837 }
2838 
2839 
2840 
2841 template <int dim>
2842 inline std::size_t
2844 {
2845  return sizeof(RefinementCase<dim>);
2846 }
2847 
2848 
2849 
2850 template <int dim>
2851 template <class Archive>
2852 inline void
2853 RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2854 {
2855  // serialization can't deal with bitfields, so copy from/to a full sized
2856  // std::uint8_t
2857  std::uint8_t uchar_value = value;
2858  ar & uchar_value;
2859  value = uchar_value;
2860 }
2861 
2862 
2863 
2864 template <>
2865 inline Point<1>
2866 GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2867 {
2869 
2870  return Point<1>(static_cast<double>(vertex));
2871 }
2872 
2873 
2874 
2875 template <>
2876 inline Point<2>
2877 GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2878 {
2880 
2881  return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2882 }
2883 
2884 
2885 
2886 template <>
2887 inline Point<3>
2888 GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2889 {
2891 
2892  return {static_cast<double>(vertex % 2),
2893  static_cast<double>(vertex / 2 % 2),
2894  static_cast<double>(vertex / 4)};
2895 }
2896 
2897 
2898 
2899 inline std::array<unsigned int, 0>
2901 {
2902  return {{}};
2903 }
2904 
2905 
2906 
2907 inline std::array<unsigned int, 1>
2909 {
2910  return {{0}};
2911 }
2912 
2913 
2914 
2915 template <int dim>
2918 {
2919  return {0U, faces_per_cell};
2920 }
2921 
2922 
2923 
2924 template <int dim>
2927 {
2928  return {0U, vertices_per_cell};
2929 }
2930 
2931 
2932 
2933 template <int dim>
2934 inline Point<dim>
2935 GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2936 {
2937  Assert(false, ExcNotImplemented());
2938 
2939  return {};
2940 }
2941 
2942 
2943 
2944 template <>
2945 inline unsigned int
2947 {
2948  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2949 
2950  return (p[0] <= 0.5 ? 0 : 1);
2951 }
2952 
2953 
2954 
2955 template <>
2956 inline unsigned int
2958 {
2959  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2960  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2961 
2962  return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
2963 }
2964 
2965 
2966 
2967 template <>
2968 inline unsigned int
2970 {
2971  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2972  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2973  Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
2974 
2975  return (p[0] <= 0.5 ?
2976  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
2977  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
2978 }
2979 
2980 
2981 template <int dim>
2982 inline unsigned int
2984 {
2985  Assert(false, ExcNotImplemented());
2986 
2987  return 0;
2988 }
2989 
2990 
2991 
2992 template <>
2993 inline Point<1>
2995  const unsigned int child_index,
2996  const RefinementCase<1> refine_case)
2997 
2998 {
2999  AssertIndexRange(child_index, 2);
3001  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3002 
3003  return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
3004 }
3005 
3006 
3007 
3008 template <>
3009 inline Point<2>
3011  const unsigned int child_index,
3012  const RefinementCase<2> refine_case)
3013 
3014 {
3015  AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3016 
3017  Point<2> point = p;
3018  switch (refine_case)
3019  {
3021  point[0] *= 2.0;
3022  if (child_index == 1)
3023  point[0] -= 1.0;
3024  break;
3026  point[1] *= 2.0;
3027  if (child_index == 1)
3028  point[1] -= 1.0;
3029  break;
3031  point *= 2.0;
3032  point -= unit_cell_vertex(child_index);
3033  break;
3034  default:
3035  Assert(false, ExcInternalError());
3036  }
3037 
3038  return point;
3039 }
3040 
3041 
3042 
3043 template <>
3044 inline Point<3>
3046  const unsigned int child_index,
3047  const RefinementCase<3> refine_case)
3048 
3049 {
3050  AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3051 
3052  Point<3> point = p;
3053  // there might be a cleverer way to do
3054  // this, but since this function is called
3055  // in very few places for initialization
3056  // purposes only, I don't care at the
3057  // moment
3058  switch (refine_case)
3059  {
3061  point[0] *= 2.0;
3062  if (child_index == 1)
3063  point[0] -= 1.0;
3064  break;
3066  point[1] *= 2.0;
3067  if (child_index == 1)
3068  point[1] -= 1.0;
3069  break;
3071  point[2] *= 2.0;
3072  if (child_index == 1)
3073  point[2] -= 1.0;
3074  break;
3076  point[0] *= 2.0;
3077  point[1] *= 2.0;
3078  if (child_index % 2 == 1)
3079  point[0] -= 1.0;
3080  if (child_index / 2 == 1)
3081  point[1] -= 1.0;
3082  break;
3084  // careful, this is slightly
3085  // different from xy and yz due to
3086  // different internal numbering of
3087  // children!
3088  point[0] *= 2.0;
3089  point[2] *= 2.0;
3090  if (child_index / 2 == 1)
3091  point[0] -= 1.0;
3092  if (child_index % 2 == 1)
3093  point[2] -= 1.0;
3094  break;
3096  point[1] *= 2.0;
3097  point[2] *= 2.0;
3098  if (child_index % 2 == 1)
3099  point[1] -= 1.0;
3100  if (child_index / 2 == 1)
3101  point[2] -= 1.0;
3102  break;
3104  point *= 2.0;
3105  point -= unit_cell_vertex(child_index);
3106  break;
3107  default:
3108  Assert(false, ExcInternalError());
3109  }
3110 
3111  return point;
3112 }
3113 
3114 
3115 
3116 template <int dim>
3117 inline Point<dim>
3119  const Point<dim> & /*p*/,
3120  const unsigned int /*child_index*/,
3121  const RefinementCase<dim> /*refine_case*/)
3122 
3123 {
3124  Assert(false, ExcNotImplemented());
3125  return {};
3126 }
3127 
3128 
3129 
3130 template <>
3131 inline Point<1>
3133  const unsigned int child_index,
3134  const RefinementCase<1> refine_case)
3135 
3136 {
3137  AssertIndexRange(child_index, 2);
3139  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3140 
3141  return (p + unit_cell_vertex(child_index)) * 0.5;
3142 }
3143 
3144 
3145 
3146 template <>
3147 inline Point<3>
3149  const unsigned int child_index,
3150  const RefinementCase<3> refine_case)
3151 
3152 {
3153  AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3154 
3155  Point<3> point = p;
3156  // there might be a cleverer way to do
3157  // this, but since this function is called
3158  // in very few places for initialization
3159  // purposes only, I don't care at the
3160  // moment
3161  switch (refine_case)
3162  {
3164  if (child_index == 1)
3165  point[0] += 1.0;
3166  point[0] *= 0.5;
3167  break;
3169  if (child_index == 1)
3170  point[1] += 1.0;
3171  point[1] *= 0.5;
3172  break;
3174  if (child_index == 1)
3175  point[2] += 1.0;
3176  point[2] *= 0.5;
3177  break;
3179  if (child_index % 2 == 1)
3180  point[0] += 1.0;
3181  if (child_index / 2 == 1)
3182  point[1] += 1.0;
3183  point[0] *= 0.5;
3184  point[1] *= 0.5;
3185  break;
3187  // careful, this is slightly
3188  // different from xy and yz due to
3189  // different internal numbering of
3190  // children!
3191  if (child_index / 2 == 1)
3192  point[0] += 1.0;
3193  if (child_index % 2 == 1)
3194  point[2] += 1.0;
3195  point[0] *= 0.5;
3196  point[2] *= 0.5;
3197  break;
3199  if (child_index % 2 == 1)
3200  point[1] += 1.0;
3201  if (child_index / 2 == 1)
3202  point[2] += 1.0;
3203  point[1] *= 0.5;
3204  point[2] *= 0.5;
3205  break;
3207  point += unit_cell_vertex(child_index);
3208  point *= 0.5;
3209  break;
3210  default:
3211  Assert(false, ExcInternalError());
3212  }
3213 
3214  return point;
3215 }
3216 
3217 
3218 
3219 template <>
3220 inline Point<2>
3222  const unsigned int child_index,
3223  const RefinementCase<2> refine_case)
3224 {
3225  AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3226 
3227  Point<2> point = p;
3228  switch (refine_case)
3229  {
3231  if (child_index == 1)
3232  point[0] += 1.0;
3233  point[0] *= 0.5;
3234  break;
3236  if (child_index == 1)
3237  point[1] += 1.0;
3238  point[1] *= 0.5;
3239  break;
3241  point += unit_cell_vertex(child_index);
3242  point *= 0.5;
3243  break;
3244  default:
3245  Assert(false, ExcInternalError());
3246  }
3247 
3248  return point;
3249 }
3250 
3251 
3252 
3253 template <int dim>
3254 inline Point<dim>
3256  const Point<dim> & /*p*/,
3257  const unsigned int /*child_index*/,
3258  const RefinementCase<dim> /*refine_case*/)
3259 {
3260  Assert(false, ExcNotImplemented());
3261  return {};
3262 }
3263 
3264 
3265 
3266 template <int dim>
3267 inline bool
3269 {
3270  Assert(false, ExcNotImplemented());
3271  return false;
3272 }
3273 
3274 template <>
3275 inline bool
3277 {
3278  return (p[0] >= 0.) && (p[0] <= 1.);
3279 }
3280 
3281 
3282 
3283 template <>
3284 inline bool
3286 {
3287  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
3288 }
3289 
3290 
3291 
3292 template <>
3293 inline bool
3295 {
3296  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
3297  (p[2] >= 0.) && (p[2] <= 1.);
3298 }
3299 
3300 
3301 
3302 template <int dim>
3303 inline bool
3305 {
3306  Assert(false, ExcNotImplemented());
3307  return false;
3308 }
3309 
3310 template <>
3311 inline bool
3312 GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
3313 {
3314  return (p[0] >= -eps) && (p[0] <= 1. + eps);
3315 }
3316 
3317 
3318 
3319 template <>
3320 inline bool
3321 GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
3322 {
3323  const double l = -eps, u = 1 + eps;
3324  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
3325 }
3326 
3327 
3328 
3329 template <>
3330 inline bool
3331 GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
3332 {
3333  const double l = -eps, u = 1.0 + eps;
3334  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
3335  (p[2] >= l) && (p[2] <= u);
3336 }
3337 
3338 
3339 
3340 template <>
3341 inline unsigned int
3342 GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
3343  const unsigned int vertex)
3344 {
3345  (void)line;
3347  AssertIndexRange(vertex, 2);
3348 
3349  return vertex;
3350 }
3351 
3352 
3353 template <>
3354 inline unsigned int
3355 GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
3356  const unsigned int vertex)
3357 {
3358  constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
3359  return cell_vertices[line][vertex];
3360 }
3361 
3362 
3363 
3364 template <>
3365 inline unsigned int
3366 GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
3367  const unsigned int vertex)
3368 {
3370  AssertIndexRange(vertex, 2);
3371 
3372  constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3373  {1, 3},
3374  {0, 1},
3375  {2, 3},
3376  {4, 6}, // top face
3377  {5, 7},
3378  {4, 5},
3379  {6, 7},
3380  {0,
3381  4}, // connects of bottom
3382  {1, 5}, // top face
3383  {2, 6},
3384  {3, 7}};
3385  return vertices[line][vertex];
3386 }
3387 
3388 
3389 template <>
3390 inline unsigned int
3391 GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3392 {
3393  Assert(false, ExcNotImplemented());
3395 }
3396 
3397 template <>
3398 inline unsigned int
3399 GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3400  const bool face_orientation,
3401  const bool face_flip,
3402  const bool face_rotation)
3403 {
3405 
3406  // set up a table to make sure that
3407  // we handle non-standard faces correctly
3408  //
3409  // so set up a table that for each vertex (of
3410  // a quad in standard position) describes
3411  // which vertex to take
3412  //
3413  // first index: four vertices 0...3
3414  //
3415  // second index: face_orientation; 0:
3416  // opposite normal, 1: standard
3417  //
3418  // third index: face_flip; 0: standard, 1:
3419  // face rotated by 180 degrees
3420  //
3421  // forth index: face_rotation: 0: standard,
3422  // 1: face rotated by 90 degrees
3423 
3424  constexpr unsigned int vertex_translation[4][2][2][2] = {
3425  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3426  // face_rotation=false and true
3427  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3428  // face_rotation=false and true
3429  {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3430  // face_rotation=false and true
3431  {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3432  // face_rotation=false and true
3433 
3434  {{{2, 3}, // vertex 1 ...
3435  {1, 0}},
3436  {{1, 0}, {2, 3}}},
3437 
3438  {{{1, 0}, // vertex 2 ...
3439  {2, 3}},
3440  {{2, 3}, {1, 0}}},
3441 
3442  {{{3, 1}, // vertex 3 ...
3443  {0, 2}},
3444  {{3, 1}, {0, 2}}}};
3445 
3446  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3447 }
3448 
3449 
3450 
3451 template <int dim>
3452 inline unsigned int
3453 GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3454  const bool,
3455  const bool,
3456  const bool)
3457 {
3458  Assert(dim > 1, ExcImpossibleInDim(dim));
3460  return vertex;
3461 }
3462 
3463 template <int dim>
3464 inline unsigned int
3466 {
3467  constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3468  0, 2, 2, 4, 2, 4, 4, 8};
3469 
3470  return n_children[ref_case];
3471 }
3472 
3473 
3474 
3475 template <int dim>
3476 inline unsigned int
3478 {
3479  Assert(false, ExcNotImplemented());
3480  return 0;
3481 }
3482 
3483 template <>
3484 inline unsigned int
3486 {
3487  Assert(false, ExcImpossibleInDim(1));
3488  return 0;
3489 }
3490 
3491 template <>
3492 inline unsigned int
3494 {
3495  return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3496 }
3497 
3498 
3499 
3500 template <>
3501 inline unsigned int
3503 {
3504  const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3505  0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3506  return nsubs[subface_case];
3507 }
3508 
3509 
3510 
3511 template <int dim>
3512 inline double
3514  const unsigned int)
3515 {
3516  Assert(false, ExcNotImplemented());
3517  return 0.;
3518 }
3519 
3520 template <>
3521 inline double
3523  const unsigned int)
3524 {
3525  return 1;
3526 }
3527 
3528 
3529 template <>
3530 inline double
3532  const unsigned int)
3533 {
3534  double ratio = 1;
3535  switch (subface_case)
3536  {
3538  // Here, an
3539  // Assert(false,ExcInternalError())
3540  // would be the right
3541  // choice, but
3542  // unfortunately the
3543  // current function is
3544  // also called for faces
3545  // without children (see
3546  // tests/fe/mapping.cc).
3547  // Assert(false, ExcMessage("Face has no subfaces."));
3548  // Furthermore, assign
3549  // following value as
3550  // otherwise the
3551  // bits/volume_x tests
3552  // break
3554  break;
3556  ratio = 0.5;
3557  break;
3558  default:
3559  // there should be no
3560  // cases left
3561  Assert(false, ExcInternalError());
3562  break;
3563  }
3564 
3565  return ratio;
3566 }
3567 
3568 
3569 template <>
3570 inline double
3572  const unsigned int subface_no)
3573 {
3574  double ratio = 1;
3575  switch (subface_case)
3576  {
3578  // Here, an
3579  // Assert(false,ExcInternalError())
3580  // would be the right
3581  // choice, but
3582  // unfortunately the
3583  // current function is
3584  // also called for faces
3585  // without children (see
3586  // tests/bits/mesh_3d_16.cc). Add
3587  // following switch to
3588  // avoid diffs in
3589  // tests/bits/mesh_3d_16
3591  break;
3594  ratio = 0.5;
3595  break;
3599  ratio = 0.25;
3600  break;
3603  if (subface_no < 2)
3604  ratio = 0.25;
3605  else
3606  ratio = 0.5;
3607  break;
3610  if (subface_no == 0)
3611  ratio = 0.5;
3612  else
3613  ratio = 0.25;
3614  break;
3615  default:
3616  // there should be no
3617  // cases left
3618  Assert(false, ExcInternalError());
3619  break;
3620  }
3621 
3622  return ratio;
3623 }
3624 
3625 
3626 
3627 template <int dim>
3629  const RefinementCase<dim> &,
3630  const unsigned int,
3631  const bool,
3632  const bool,
3633  const bool)
3634 {
3635  Assert(false, ExcNotImplemented());
3636  return RefinementCase<dim - 1>::no_refinement;
3637 }
3638 
3639 template <>
3641  const RefinementCase<1> &,
3642  const unsigned int,
3643  const bool,
3644  const bool,
3645  const bool)
3646 {
3647  Assert(false, ExcImpossibleInDim(1));
3648 
3650 }
3651 
3652 
3653 template <>
3654 inline RefinementCase<1>
3656  const RefinementCase<2> &cell_refinement_case,
3657  const unsigned int face_no,
3658  const bool,
3659  const bool,
3660  const bool)
3661 {
3662  const unsigned int dim = 2;
3663  AssertIndexRange(cell_refinement_case,
3666 
3667  const RefinementCase<dim - 1>
3670  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3671  RefinementCase<dim - 1>::no_refinement},
3672 
3673  {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3674 
3675  {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3676 
3677  {RefinementCase<dim - 1>::cut_x, // cut_xy
3678  RefinementCase<dim - 1>::cut_x}};
3679 
3680  return ref_cases[cell_refinement_case][face_no / 2];
3681 }
3682 
3683 
3684 template <>
3685 inline RefinementCase<2>
3687  const RefinementCase<3> &cell_refinement_case,
3688  const unsigned int face_no,
3689  const bool face_orientation,
3690  const bool /*face_flip*/,
3691  const bool face_rotation)
3692 {
3693  const unsigned int dim = 3;
3694  AssertIndexRange(cell_refinement_case,
3697 
3698  const RefinementCase<dim - 1>
3701  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3702  RefinementCase<dim - 1>::no_refinement,
3703  RefinementCase<dim - 1>::no_refinement},
3704 
3705  {RefinementCase<dim - 1>::no_refinement, // cut_x
3706  RefinementCase<dim - 1>::cut_y,
3707  RefinementCase<dim - 1>::cut_x},
3708 
3709  {RefinementCase<dim - 1>::cut_x, // cut_y
3710  RefinementCase<dim - 1>::no_refinement,
3711  RefinementCase<dim - 1>::cut_y},
3712 
3713  {RefinementCase<dim - 1>::cut_x, // cut_xy
3714  RefinementCase<dim - 1>::cut_y,
3715  RefinementCase<dim - 1>::cut_xy},
3716 
3717  {RefinementCase<dim - 1>::cut_y, // cut_z
3718  RefinementCase<dim - 1>::cut_x,
3719  RefinementCase<dim - 1>::no_refinement},
3720 
3721  {RefinementCase<dim - 1>::cut_y, // cut_xz
3722  RefinementCase<dim - 1>::cut_xy,
3723  RefinementCase<dim - 1>::cut_x},
3724 
3725  {RefinementCase<dim - 1>::cut_xy, // cut_yz
3726  RefinementCase<dim - 1>::cut_x,
3727  RefinementCase<dim - 1>::cut_y},
3728 
3729  {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3730  RefinementCase<dim - 1>::cut_xy,
3731  RefinementCase<dim - 1>::cut_xy},
3732  };
3733 
3734  const RefinementCase<dim - 1> ref_case =
3735  ref_cases[cell_refinement_case][face_no / 2];
3736 
3737  const RefinementCase<dim - 1> flip[4] = {
3738  RefinementCase<dim - 1>::no_refinement,
3739  RefinementCase<dim - 1>::cut_y,
3740  RefinementCase<dim - 1>::cut_x,
3741  RefinementCase<dim - 1>::cut_xy};
3742 
3743  // correct the ref_case for face_orientation
3744  // and face_rotation. for face_orientation,
3745  // 'true' is the default value whereas for
3746  // face_rotation, 'false' is standard. If
3747  // <tt>face_rotation==face_orientation</tt>,
3748  // then one of them is non-standard and we
3749  // have to swap cut_x and cut_y, otherwise no
3750  // change is necessary. face_flip has no
3751  // influence. however, in order to keep the
3752  // interface consistent with other functions,
3753  // we still include it as an argument to this
3754  // function
3755  return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3756 }
3757 
3758 
3759 
3760 template <int dim>
3761 inline RefinementCase<1>
3763  const unsigned int)
3764 {
3765  Assert(false, ExcNotImplemented());
3767 }
3768 
3769 template <>
3770 inline RefinementCase<1>
3772  const RefinementCase<1> &cell_refinement_case,
3773  const unsigned int line_no)
3774 {
3775  (void)line_no;
3776  const unsigned int dim = 1;
3777  (void)dim;
3778  AssertIndexRange(cell_refinement_case,
3781 
3782  return cell_refinement_case;
3783 }
3784 
3785 
3786 template <>
3787 inline RefinementCase<1>
3789  const RefinementCase<2> &cell_refinement_case,
3790  const unsigned int line_no)
3791 {
3792  // Assertions are in face_refinement_case()
3793  return face_refinement_case(cell_refinement_case, line_no);
3794 }
3795 
3796 
3797 template <>
3798 inline RefinementCase<1>
3800  const RefinementCase<3> &cell_refinement_case,
3801  const unsigned int line_no)
3802 {
3803  const unsigned int dim = 3;
3804  AssertIndexRange(cell_refinement_case,
3807 
3808  // array indicating, which simple refine
3809  // case cuts a line in direction x, y or
3810  // z. For example, cut_y and everything
3811  // containing cut_y (cut_xy, cut_yz,
3812  // cut_xyz) cuts lines, which are in y
3813  // direction.
3814  const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3817 
3818  // order the direction of lines
3819  // 0->x, 1->y, 2->z
3820  const unsigned int direction[lines_per_cell] = {
3821  1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3822 
3823  return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3826 }
3827 
3828 
3829 
3830 template <int dim>
3831 inline RefinementCase<dim>
3833  const RefinementCase<dim - 1> &,
3834  const unsigned int,
3835  const bool,
3836  const bool,
3837  const bool)
3838 {
3839  Assert(false, ExcNotImplemented());
3840 
3842 }
3843 
3844 template <>
3845 inline RefinementCase<1>
3847  const RefinementCase<0> &,
3848  const unsigned int,
3849  const bool,
3850  const bool,
3851  const bool)
3852 {
3853  const unsigned int dim = 1;
3854  Assert(false, ExcImpossibleInDim(dim));
3855 
3857 }
3858 
3859 
3860 template <>
3861 inline RefinementCase<2>
3864  const unsigned int face_no,
3865  const bool,
3866  const bool,
3867  const bool)
3868 {
3869  const unsigned int dim = 2;
3870  AssertIndexRange(face_refinement_case,
3873 
3874  if (face_refinement_case == RefinementCase<dim>::cut_x)
3875  return (face_no / 2) ? RefinementCase<dim>::cut_x :
3877  else
3879 }
3880 
3881 
3882 template <>
3883 inline RefinementCase<3>
3885  const RefinementCase<2> &face_refinement_case,
3886  const unsigned int face_no,
3887  const bool face_orientation,
3888  const bool /*face_flip*/,
3889  const bool face_rotation)
3890 {
3891  const unsigned int dim = 3;
3892  AssertIndexRange(face_refinement_case,
3895 
3900 
3901  // correct the face_refinement_case for
3902  // face_orientation and face_rotation. for
3903  // face_orientation, 'true' is the default
3904  // value whereas for face_rotation, 'false'
3905  // is standard. If
3906  // <tt>face_rotation==face_orientation</tt>,
3907  // then one of them is non-standard and we
3908  // have to swap cut_x and cut_y, otherwise no
3909  // change is necessary. face_flip has no
3910  // influence. however, in order to keep the
3911  // interface consistent with other functions,
3912  // we still include it as an argument to this
3913  // function
3914  const RefinementCase<dim - 1> std_face_ref =
3915  (face_orientation == face_rotation) ? flip[face_refinement_case] :
3916  face_refinement_case;
3917 
3918  const RefinementCase<dim> face_to_cell[3][4] = {
3919  {RefinementCase<dim>::no_refinement, // faces 0 and 1
3920  RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3921  RefinementCase<dim>::cut_z,
3923 
3924  {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3925  // "exchanged on faces 2 and 3")
3926  RefinementCase<dim>::cut_z,
3929 
3930  {RefinementCase<dim>::no_refinement, // faces 4 and 5
3934 
3935  return face_to_cell[face_no / 2][std_face_ref];
3936 }
3937 
3938 
3939 
3940 template <int dim>
3941 inline RefinementCase<dim>
3943  const unsigned int)
3944 {
3945  Assert(false, ExcNotImplemented());
3946 
3948 }
3949 
3950 template <>
3951 inline RefinementCase<1>
3953  const unsigned int line_no)
3954 {
3955  (void)line_no;
3956  AssertIndexRange(line_no, 1);
3957 
3958  return RefinementCase<1>::cut_x;
3959 }
3960 
3961 
3962 template <>
3963 inline RefinementCase<2>
3965  const unsigned int line_no)
3966 {
3967  const unsigned int dim = 2;
3968  (void)dim;
3970 
3971  return (line_no / 2) ? RefinementCase<2>::cut_x : RefinementCase<2>::cut_y;
3972 }
3973 
3974 
3975 template <>
3976 inline RefinementCase<3>
3978  const unsigned int line_no)
3979 {
3980  const unsigned int dim = 3;
3982 
3983  const RefinementCase<dim> ref_cases[6] = {
3984  RefinementCase<dim>::cut_y, // lines 0 and 1
3985  RefinementCase<dim>::cut_x, // lines 2 and 3
3986  RefinementCase<dim>::cut_y, // lines 4 and 5
3987  RefinementCase<dim>::cut_x, // lines 6 and 7
3988  RefinementCase<dim>::cut_z, // lines 8 and 9
3989  RefinementCase<dim>::cut_z}; // lines 10 and 11
3990 
3991  return ref_cases[line_no / 2];
3992 }
3993 
3994 
3995 
3996 template <>
3997 inline unsigned int
3998 GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
3999  const bool face_orientation,
4000  const bool face_flip,
4001  const bool face_rotation)
4002 {
4004 
4005  // set up a table to make sure that
4006  // we handle non-standard faces correctly
4007  //
4008  // so set up a table that for each vertex (of
4009  // a quad in standard position) describes
4010  // which vertex to take
4011  //
4012  // first index: four vertices 0...3
4013  //
4014  // second index: face_orientation; 0:
4015  // opposite normal, 1: standard
4016  //
4017  // third index: face_flip; 0: standard, 1:
4018  // face rotated by 180 degrees
4019  //
4020  // forth index: face_rotation: 0: standard,
4021  // 1: face rotated by 90 degrees
4022 
4023  const unsigned int vertex_translation[4][2][2][2] = {
4024  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
4025  // face_rotation=false and true
4026  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
4027  // face_rotation=false and true
4028  {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
4029  // face_rotation=false and true
4030  {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
4031  // face_rotation=false and true
4032 
4033  {{{2, 3}, // vertex 1 ...
4034  {1, 0}},
4035  {{1, 3}, {2, 0}}},
4036 
4037  {{{1, 0}, // vertex 2 ...
4038  {2, 3}},
4039  {{2, 0}, {1, 3}}},
4040 
4041  {{{3, 1}, // vertex 3 ...
4042  {0, 2}},
4043  {{3, 2}, {0, 1}}}};
4044 
4045  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
4046 }
4047 
4048 
4049 
4050 template <int dim>
4051 inline unsigned int
4052 GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
4053  const bool,
4054  const bool,
4055  const bool)
4056 {
4057  Assert(dim > 1, ExcImpossibleInDim(dim));
4059  return vertex;
4060 }
4061 
4062 
4063 
4064 template <>
4065 inline unsigned int
4066 GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
4067  const bool face_orientation,
4068  const bool face_flip,
4069  const bool face_rotation)
4070 {
4072 
4073 
4074  // make sure we handle
4075  // non-standard faces correctly
4076  //
4077  // so set up a table that for each line (of a
4078  // quad) describes which line to take
4079  //
4080  // first index: four lines 0...3
4081  //
4082  // second index: face_orientation; 0:
4083  // opposite normal, 1: standard
4084  //
4085  // third index: face_flip; 0: standard, 1:
4086  // face rotated by 180 degrees
4087  //
4088  // forth index: face_rotation: 0: standard,
4089  // 1: face rotated by 90 degrees
4090 
4091  const unsigned int line_translation[4][2][2][2] = {
4092  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4093  // face_rotation=false and true
4094  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4095  // face_rotation=false and true
4096  {{0, 3}, // line 0, face_orientation=true, face_flip=false,
4097  // face_rotation=false and true
4098  {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
4099  // face_rotation=false and true
4100 
4101  {{{3, 1}, // line 1 ...
4102  {2, 0}},
4103  {{1, 2}, {0, 3}}},
4104 
4105  {{{0, 3}, // line 2 ...
4106  {1, 2}},
4107  {{2, 0}, {3, 1}}},
4108 
4109  {{{1, 2}, // line 3 ...
4110  {0, 3}},
4111  {{3, 1}, {2, 0}}}};
4112 
4113  return line_translation[line][face_orientation][face_flip][face_rotation];
4114 }
4115 
4116 
4117 
4118 template <int dim>
4119 inline unsigned int
4120 GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
4121  const bool,
4122  const bool,
4123  const bool)
4124 {
4125  Assert(false, ExcNotImplemented());
4126  return line;
4127 }
4128 
4129 
4130 
4131 template <>
4132 inline unsigned int
4133 GeometryInfo<2>::standard_to_real_line_vertex(const unsigned int vertex,
4134  const bool line_orientation)
4135 {
4136  return line_orientation ? vertex : (1 - vertex);
4137 }
4138 
4139 
4140 
4141 template <int dim>
4142 inline unsigned int
4143 GeometryInfo<dim>::standard_to_real_line_vertex(const unsigned int vertex,
4144  const bool)
4145 {
4146  Assert(false, ExcNotImplemented());
4147  return vertex;
4148 }
4149 
4150 
4151 
4152 template <>
4153 inline std::array<unsigned int, 2>
4155  const unsigned int vertex)
4156 {
4157  return {{vertex % 2, vertex / 2}};
4158 }
4159 
4160 
4161 
4162 template <int dim>
4163 inline std::array<unsigned int, 2>
4165  const unsigned int vertex)
4166 {
4167  Assert(false, ExcNotImplemented());
4168  (void)vertex;
4169  return {{0, 0}};
4170 }
4171 
4172 
4173 
4174 template <>
4175 inline std::array<unsigned int, 2>
4177 {
4178  // set up a table that for each
4179  // line describes a) from which
4180  // quad to take it, b) which line
4181  // therein it is if the face is
4182  // oriented correctly
4183  static const unsigned int lookup_table[GeometryInfo<3>::lines_per_cell][2] = {
4184  {4, 0}, // take first four lines from bottom face
4185  {4, 1},
4186  {4, 2},
4187  {4, 3},
4188 
4189  {5, 0}, // second four lines from top face
4190  {5, 1},
4191  {5, 2},
4192  {5, 3},
4193 
4194  {0, 0}, // the rest randomly
4195  {1, 0},
4196  {0, 1},
4197  {1, 1}};
4198 
4199  return {{lookup_table[i][0], lookup_table[i][1]}};
4200 }
4201 
4202 
4203 
4204 template <int dim>
4205 inline std::array<unsigned int, 2>
4207 {
4208  Assert(false, ExcNotImplemented());
4209  (void)line;
4210  return {{0, 0}};
4211 }
4212 
4213 
4214 
4215 template <>
4216 inline std::array<unsigned int, 2>
4218  const unsigned int vertex)
4219 {
4220  // get the corner indices by asking either the bottom or the top face for its
4221  // vertices. handle non-standard faces by calling the vertex reordering
4222  // function from GeometryInfo
4223 
4224  // bottom face (4) for first four vertices, top face (5) for the rest
4225  return {{4 + vertex / 4, vertex % 4}};
4226 }
4227 
4228 
4229 
4230 template <int dim>
4231 inline std::array<unsigned int, 2>
4233  const unsigned int vertex)
4234 {
4235  Assert(false, ExcNotImplemented());
4236  (void)vertex;
4237  return {{0, 0}};
4238 }
4239 
4240 
4241 
4242 template <>
4243 inline unsigned int
4244 GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
4245  const bool face_orientation,
4246  const bool face_flip,
4247  const bool face_rotation)
4248 {
4250 
4251 
4252  // make sure we handle
4253  // non-standard faces correctly
4254  //
4255  // so set up a table that for each line (of a
4256  // quad) describes which line to take
4257  //
4258  // first index: four lines 0...3
4259  //
4260  // second index: face_orientation; 0:
4261  // opposite normal, 1: standard
4262  //
4263  // third index: face_flip; 0: standard, 1:
4264  // face rotated by 180 degrees
4265  //
4266  // forth index: face_rotation: 0: standard,
4267  // 1: face rotated by 90 degrees
4268 
4269  const unsigned int line_translation[4][2][2][2] = {
4270  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4271  // face_rotation=false and true
4272  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4273  // face_rotation=false and true
4274  {{0, 2}, // line 0, face_orientation=true, face_flip=false,
4275  // face_rotation=false and true
4276  {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
4277  // face_rotation=false and true
4278 
4279  {{{3, 1}, // line 1 ...
4280  {2, 0}},
4281  {{1, 3}, {0, 2}}},
4282 
4283  {{{0, 3}, // line 2 ...
4284  {1, 2}},
4285  {{2, 1}, {3, 0}}},
4286 
4287  {{{1, 2}, // line 3 ...
4288  {0, 3}},
4289  {{3, 0}, {2, 1}}}};
4290 
4291  return line_translation[line][face_orientation][face_flip][face_rotation];
4292 }
4293 
4294 
4295 
4296 template <int dim>
4297 inline unsigned int
4298 GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
4299  const bool,
4300  const bool,
4301  const bool)
4302 {
4303  Assert(false, ExcNotImplemented());
4304  return line;
4305 }
4306 
4307 
4308 
4309 template <>
4310 inline unsigned int
4312  const unsigned int face,
4313  const unsigned int subface,
4314  const bool,
4315  const bool,
4316  const bool,
4317  const RefinementCase<0> &)
4318 {
4319  (void)subface;
4322 
4323  return face;
4324 }
4325 
4326 
4327 
4328 template <>
4329 inline unsigned int
4331  const unsigned int face,
4332  const unsigned int subface,
4333  const bool /*face_orientation*/,
4334  const bool face_flip,
4335  const bool /*face_rotation*/,
4336  const RefinementCase<1> &)
4337 {
4340 
4341  // always return the child adjacent to the specified
4342  // subface. if the face of a cell is not refined, don't
4343  // throw an assertion but deliver the child adjacent to
4344  // the face nevertheless, i.e. deliver the child of
4345  // this cell adjacent to the subface of a possibly
4346  // refined neighbor. this simplifies setting neighbor
4347  // information in execute_refinement.
4348  const unsigned int
4350  [max_children_per_face] = {
4351  {
4352  // Normal orientation (face_flip = false)
4353  {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
4354  {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
4355  {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
4356  },
4357  {
4358  // Flipped orientation (face_flip = true)
4359  {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
4360  {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
4361  {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
4362  }};
4363 
4364  return subcells[face_flip][ref_case - 1][face][subface];
4365 }
4366 
4367 
4368 
4369 template <>
4370 inline unsigned int
4372  const unsigned int face,
4373  const unsigned int subface,
4374  const bool face_orientation,
4375  const bool face_flip,
4376  const bool face_rotation,
4377  const RefinementCase<2> &face_ref_case)
4378 {
4379  const unsigned int dim = 3;
4380 
4382  ExcMessage("Cell has no children."));
4384  if (!(subface == 0 &&
4385  face_ref_case == RefinementCase<dim - 1>::no_refinement))
4386  {
4387  AssertIndexRange(subface,
4388  GeometryInfo<dim - 1>::n_children(face_ref_case));
4389  }
4390 
4391  // invalid number used for invalid cases,
4392  // e.g. when the children are more refined at
4393  // a given face than the face itself
4394  const unsigned int e = numbers::invalid_unsigned_int;
4395 
4396  // the whole process of finding a child cell
4397  // at a given subface considering the
4398  // possibly anisotropic refinement cases of
4399  // the cell and the face as well as
4400  // orientation, flip and rotation of the face
4401  // is quite complicated. thus, we break it
4402  // down into several steps.
4403 
4404  // first step: convert the given face refine
4405  // case to a face refine case concerning the
4406  // face in standard orientation (, flip and
4407  // rotation). This only affects cut_x and
4408  // cut_y
4409  const RefinementCase<dim - 1> flip[4] = {
4410  RefinementCase<dim - 1>::no_refinement,
4411  RefinementCase<dim - 1>::cut_y,
4412  RefinementCase<dim - 1>::cut_x,
4413  RefinementCase<dim - 1>::cut_xy};
4414  // for face_orientation, 'true' is the
4415  // default value whereas for face_rotation,
4416  // 'false' is standard. If
4417  // <tt>face_rotation==face_orientation</tt>,
4418  // then one of them is non-standard and we
4419  // have to swap cut_x and cut_y, otherwise no
4420  // change is necessary.
4421  const RefinementCase<dim - 1> std_face_ref =
4422  (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4423 
4424  // second step: convert the given subface
4425  // index to the one for a standard face
4426  // respecting face_orientation, face_flip and
4427  // face_rotation
4428 
4429  // first index: face_ref_case
4430  // second index: face_orientation
4431  // third index: face_flip
4432  // forth index: face_rotation
4433  // fifth index: subface index
4434  const unsigned int subface_exchange[4][2][2][2][4] = {
4435  // no_refinement (subface 0 stays 0,
4436  // all others are invalid)
4437  {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4438  {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4439  // cut_x (here, if the face is only
4440  // rotated OR only falsely oriented,
4441  // then subface 0 of the non-standard
4442  // face does NOT correspond to one of
4443  // the subfaces of a standard
4444  // face. Thus we indicate the subface
4445  // which is located at the lower left
4446  // corner (the origin of the face's
4447  // local coordinate system) with
4448  // '0'. The rest of this issue is
4449  // taken care of using the above
4450  // conversion to a 'standard face
4451  // refine case')
4452  {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4453  {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4454  // cut_y (the same applies as for
4455  // cut_x)
4456  {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4457  {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4458  // cut_xyz: this information is
4459  // identical to the information
4460  // returned by
4461  // GeometryInfo<3>::real_to_standard_face_vertex()
4462  {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4463  // face_rotation=false, subfaces 0,1,2,3
4464  {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4465  // face_rotation=true, subfaces 0,1,2,3
4466  {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4467  // face_rotation=false, subfaces 0,1,2,3
4468  {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4469  // face_rotation=true, subfaces 0,1,2,3
4470  {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4471  // face_rotation=false, subfaces 0,1,2,3
4472  {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4473  // face_rotation=true, subfaces 0,1,2,3
4474  {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4475  // face_rotation=false, subfaces 0,1,2,3
4476  {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4477  // face_rotation=true, subfaces 0,1,2,3
4478 
4479  const unsigned int std_subface =
4480  subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4481  [subface];
4482  Assert(std_subface != e, ExcInternalError());
4483 
4484  // third step: these are the children, which
4485  // can be found at the given subfaces of an
4486  // isotropically refined (standard) face
4487  //
4488  // first index: (refinement_case-1)
4489  // second index: face_index
4490  // third index: subface_index (isotropic refinement)
4491  const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4492  [max_children_per_face] = {
4493  // cut_x
4494  {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4495  {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4496  {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4497  {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4498  {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4499  {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4500  // cut_y
4501  {{0, 1, 0, 1},
4502  {0, 1, 0, 1},
4503  {0, 0, 0, 0},
4504  {1, 1, 1, 1},
4505  {0, 0, 1, 1},
4506  {0, 0, 1, 1}},
4507  // cut_xy
4508  {{0, 2, 0, 2},
4509  {1, 3, 1, 3},
4510  {0, 0, 1, 1},
4511  {2, 2, 3, 3},
4512  {0, 1, 2, 3},
4513  {0, 1, 2, 3}},
4514  // cut_z
4515  {{0, 0, 1, 1},
4516  {0, 0, 1, 1},
4517  {0, 1, 0, 1},
4518  {0, 1, 0, 1},
4519  {0, 0, 0, 0},
4520  {1, 1, 1, 1}},
4521  // cut_xz
4522  {{0, 0, 1, 1},
4523  {2, 2, 3, 3},
4524  {0, 1, 2, 3},
4525  {0, 1, 2, 3},
4526  {0, 2, 0, 2},
4527  {1, 3, 1, 3}},
4528  // cut_yz
4529  {{0, 1, 2, 3},
4530  {0, 1, 2, 3},
4531  {0, 2, 0, 2},
4532  {1, 3, 1, 3},
4533  {0, 0, 1, 1},
4534  {2, 2, 3, 3}},
4535  // cut_xyz
4536  {{0, 2, 4, 6},
4537  {1, 3, 5, 7},
4538  {0, 4, 1, 5},
4539  {2, 6, 3, 7},
4540  {0, 1, 2, 3},
4541  {4, 5, 6, 7}}};
4542 
4543  // forth step: check, whether the given face
4544  // refine case is valid for the given cell
4545  // refine case. this is the case, if the
4546  // given face refine case is at least as
4547  // refined as the face is for the given cell
4548  // refine case
4549 
4550  // note, that we are considering standard
4551  // face refinement cases here and thus must
4552  // not pass the given orientation, flip and
4553  // rotation flags
4554  if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4555  face_refinement_case(ref_case, face))
4556  {
4557  // all is fine. for anisotropic face
4558  // refine cases, select one of the
4559  // isotropic subfaces which neighbors the
4560  // same child
4561 
4562  // first index: (standard) face refine case
4563  // second index: subface index
4564  const unsigned int equivalent_iso_subface[4][4] = {
4565  {0, e, e, e}, // no_refinement
4566  {0, 3, e, e}, // cut_x
4567  {0, 3, e, e}, // cut_y
4568  {0, 1, 2, 3}}; // cut_xy
4569 
4570  const unsigned int equ_std_subface =
4571  equivalent_iso_subface[std_face_ref][std_subface];
4572  Assert(equ_std_subface != e, ExcInternalError());
4573 
4574  return iso_children[ref_case - 1][face][equ_std_subface];
4575  }
4576  else
4577  {
4578  // the face_ref_case was too coarse,
4579  // throw an error
4580  Assert(false,
4581  ExcMessage("The face RefineCase is too coarse "
4582  "for the given cell RefineCase."));
4583  }
4584  // we only get here in case of an error
4585  return e;
4586 }
4587 
4588 
4589 
4590 template <>
4591 inline unsigned int
4593  const unsigned int,
4594  const unsigned int,
4595  const bool,
4596  const bool,
4597  const bool,
4598  const RefinementCase<3> &)
4599 {
4600  Assert(false, ExcNotImplemented());
4602 }
4603 
4604 
4605 
4606 template <>
4607 inline unsigned int
4608 GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4609  const unsigned int line,
4610  const bool,
4611  const bool,
4612  const bool)
4613 {
4614  (void)line;
4617 
4618  // The face is a line itself.
4619  return face;
4620 }
4621 
4622 
4623 
4624 template <>
4625 inline unsigned int
4626 GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4627  const unsigned int line,
4628  const bool face_orientation,
4629  const bool face_flip,
4630  const bool face_rotation)
4631 {
4634 
4635  const unsigned lines[faces_per_cell][lines_per_face] = {
4636  {8, 10, 0, 4}, // left face
4637  {9, 11, 1, 5}, // right face
4638  {2, 6, 8, 9}, // front face
4639  {3, 7, 10, 11}, // back face
4640  {0, 1, 2, 3}, // bottom face
4641  {4, 5, 6, 7}}; // top face
4642  return lines[face][real_to_standard_face_line(
4643  line, face_orientation, face_flip, face_rotation)];
4644 }
4645 
4646 
4647 
4648 inline unsigned int
4649 GeometryInfo<0>::face_to_cell_lines(const unsigned int,
4650  const unsigned int,
4651  const bool,
4652  const bool,
4653  const bool)
4654 {
4655  Assert(false, ExcNotImplemented());
4657 }
4658 
4659 
4660 
4661 template <int dim>
4662 inline unsigned int
4663 GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4664  const unsigned int,
4665  const bool,
4666  const bool,
4667  const bool)
4668 {
4669  Assert(false, ExcNotImplemented());
4671 }
4672 
4673 
4674 
4675 template <int dim>
4676 inline unsigned int
4677 GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4678  const unsigned int vertex,
4679  const bool face_orientation,
4680  const bool face_flip,
4681  const bool face_rotation)
4682 {
4684  face,
4685  vertex,
4686  face_orientation,
4687  face_flip,
4688  face_rotation);
4689 }
4690 
4691 
4692 
4693 inline unsigned int
4694 GeometryInfo<0>::face_to_cell_vertices(const unsigned int,
4695  const unsigned int,
4696  const bool,
4697  const bool,
4698  const bool)
4699 {
4700  Assert(false, ExcNotImplemented());
4702 }
4703 
4704 
4705 
4706 template <int dim>
4707 template <typename Number>
4708 inline Point<dim, Number>
4710 {
4712  for (unsigned int i = 0; i < dim; i++)
4713  p[i] = std::min(std::max(q[i], Number(0.)), Number(1.));
4714 
4715  return p;
4716 }
4717 
4718 
4719 
4720 template <int dim>
4721 inline double
4723 {
4724  double result = 0.0;
4725 
4726  for (unsigned int i = 0; i < dim; i++)
4727  {
4728  result = std::max(result, -p[i]);
4729  result = std::max(result, p[i] - 1.);
4730  }
4731 
4732  return result;
4733 }
4734 
4735 
4736 
4737 template <int dim>
4738 inline double
4740  const unsigned int i)
4741 {
4743 
4744  switch (dim)
4745  {
4746  case 1:
4747  {
4748  const double x = xi[0];
4749  switch (i)
4750  {
4751  case 0:
4752  return 1 - x;
4753  case 1:
4754  return x;
4755  }
4756  break;
4757  }
4758 
4759  case 2:
4760  {
4761  const double x = xi[0];
4762  const double y = xi[1];
4763  switch (i)
4764  {
4765  case 0:
4766  return (1 - x) * (1 - y);
4767  case 1:
4768  return x * (1 - y);
4769  case 2:
4770  return (1 - x) * y;
4771  case 3:
4772  return x * y;
4773  }
4774  break;
4775  }
4776 
4777  case 3:
4778  {
4779  const double x = xi[0];
4780  const double y = xi[1];
4781  const double z = xi[2];
4782  switch (i)
4783  {
4784  case 0:
4785  return (1 - x) * (1 - y) * (1 - z);
4786  case 1:
4787  return x * (1 - y) * (1 - z);
4788  case 2:
4789  return (1 - x) * y * (1 - z);
4790  case 3:
4791  return x * y * (1 - z);
4792  case 4:
4793  return (1 - x) * (1 - y) * z;
4794  case 5:
4795  return x * (1 - y) * z;
4796  case 6:
4797  return (1 - x) * y * z;
4798  case 7:
4799  return x * y * z;
4800  }
4801  break;
4802  }
4803 
4804  default:
4805  Assert(false, ExcNotImplemented());
4806  }
4807  return -1e9;
4808 }
4809 
4810 
4811 
4812 template <>
4814  const Point<1> &,
4815  const unsigned int i)
4816 {
4818 
4819  switch (i)
4820  {
4821  case 0:
4822  return Point<1>(-1.);
4823  case 1:
4824  return Point<1>(1.);
4825  }
4826 
4827  return Point<1>(-1e9);
4828 }
4829 
4830 
4831 
4832 template <>
4834  const Point<2> & xi,
4835  const unsigned int i)
4836 {
4838 
4839  const double x = xi[0];
4840  const double y = xi[1];
4841  switch (i)
4842  {
4843  case 0:
4844  return Point<2>(-(1 - y), -(1 - x));
4845  case 1:
4846  return Point<2>(1 - y, -x);
4847  case 2:
4848  return Point<2>(-y, 1 - x);
4849  case 3:
4850  return Point<2>(y, x);
4851  }
4852  return Point<2>(-1e9, -1e9);
4853 }
4854 
4855 
4856 
4857 template <>
4859  const Point<3> & xi,
4860  const unsigned int i)
4861 {
4863 
4864  const double x = xi[0];
4865  const double y = xi[1];
4866  const double z = xi[2];
4867  switch (i)
4868  {
4869  case 0:
4870  return Point<3>(-(1 - y) * (1 - z),
4871  -(1 - x) * (1 - z),
4872  -(1 - x) * (1 - y));
4873  case 1:
4874  return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4875  case 2:
4876  return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4877  case 3:
4878  return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4879  case 4:
4880  return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4881  case 5:
4882  return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4883  case 6:
4884  return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4885  case 7:
4886  return Point<3>(y * z, x * z, x * y);
4887  }
4888 
4889  return Point<3>(-1e9, -1e9, -1e9);
4890 }
4891 
4892 
4893 
4894 template <int dim>
4895 inline Tensor<1, dim>
4897  const unsigned int)
4898 {
4899  Assert(false, ExcNotImplemented());
4900  return Tensor<1, dim>();
4901 }
4902 
4903 
4904 unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4905 {
4906  return 0;
4907 }
4908 
4909 
4910 namespace internal
4911 {
4912  namespace GeometryInfoHelper
4913  {
4914  // wedge product of a single
4915  // vector in 2d: we just have to
4916  // rotate it by 90 degrees to the
4917  // right
4918  inline Tensor<1, 2>
4919  wedge_product(const Tensor<1, 2> (&derivative)[1])
4920  {
4921  Tensor<1, 2> result;
4922  result[0] = derivative[0][1];
4923  result[1] = -derivative[0][0];
4924 
4925  return result;
4926  }
4927 
4928 
4929  // wedge product of 2 vectors in
4930  // 3d is the cross product
4931  inline Tensor<1, 3>
4932  wedge_product(const Tensor<1, 3> (&derivative)[2])
4933  {
4934  return cross_product_3d(derivative[0], derivative[1]);
4935  }
4936 
4937 
4938  // wedge product of dim vectors
4939  // in dim-d: that's the
4940  // determinant of the matrix
4941  template <int dim>
4942  inline Tensor<0, dim>
4943  wedge_product(const Tensor<1, dim> (&derivative)[dim])
4944  {
4945  Tensor<2, dim> jacobian;
4946  for (unsigned int i = 0; i < dim; ++i)
4947  jacobian[i] = derivative[i];
4948 
4949  return determinant(jacobian);
4950  }
4951  } // namespace GeometryInfoHelper
4952 } // namespace internal
4953 
4954 
4955 template <int dim>
4956 template <int spacedim>
4957 inline void
4959 # ifndef DEAL_II_CXX14_CONSTEXPR_BUG
4962 # else
4964 # endif
4965 {
4966  // for each of the vertices,
4967  // compute the alternating form
4968  // of the mapped unit
4969  // vectors. consider for
4970  // example the case of a quad
4971  // in spacedim==3: to do so, we
4972  // need to see how the
4973  // infinitesimal vectors
4974  // (d\xi_1,0) and (0,d\xi_2)
4975  // are transformed into
4976  // spacedim-dimensional space
4977  // and then form their cross
4978  // product (i.e. the wedge product
4979  // of two vectors). to this end, note
4980  // that
4981  // \vec x = sum_i \vec v_i phi_i(\vec xi)
4982  // so the transformed vectors are
4983  // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
4984  // and
4985  // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
4986  // which boils down to the columns
4987  // of the 3x2 matrix \grad_\xi x(\xi)
4988  //
4989  // a similar reasoning would
4990  // hold for all dim,spacedim
4991  // pairs -- we only have to
4992  // compute the wedge product of
4993  // the columns of the
4994  // derivatives
4995  for (unsigned int i = 0; i < vertices_per_cell; ++i)
4996  {
4997  Tensor<1, spacedim> derivatives[dim];
4998 
4999  for (unsigned int j = 0; j < vertices_per_cell; ++j)
5000  {
5001  const Tensor<1, dim> grad_phi_j =
5003  for (unsigned int l = 0; l < dim; ++l)
5004  derivatives[l] += vertices[j] * grad_phi_j[l];
5005  }
5006 
5007  forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
5008  }
5009 }
5010 
5011 #endif // DOXYGEN
5013 
5014 #endif
static Point< dim > child_to_cell_coordinates(const Point< dim > &p, const unsigned int child_index, const RefinementCase< dim > refine_case=RefinementCase< dim >::isotropic_refinement)
static unsigned int real_to_standard_face_vertex(const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
static const unsigned int invalid_unsigned_int
Definition: types.h:196
static unsigned int standard_to_real_face_line(const unsigned int line, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
static unsigned int child_cell_from_point(const Point< dim > &p)
static unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static RefinementCase< 1 > line_refinement_case(const RefinementCase< dim > &cell_refinement_case, const unsigned int line_no)
static unsigned int face_to_cell_lines(const unsigned int face, const unsigned int line, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static constexpr unsigned int max_children_per_face
static bool is_inside_unit_cell(const Point< dim > &p)
GeometryPrimitive(const Object object)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
void serialize(Archive &ar, const unsigned int version)
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
static Point< dim, Number > project_to_unit_cell(const Point< dim, Number > &p)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1691
RefinementCase operator &(const RefinementCase &r) const
static unsigned int real_to_standard_face_line(const unsigned int line, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static Point< dim > cell_to_child_coordinates(const Point< dim > &p, const unsigned int child_index, const RefinementCase< dim > refine_case=RefinementCase< dim >::isotropic_refinement)
static const char U
static Point< dim > unit_cell_vertex(const unsigned int vertex)
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
static const std::array< unsigned int, vertices_per_cell > ucd_to_deal
static constexpr unsigned int vertices_per_cell
RefinementCase operator~() const
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2443
static double distance_to_unit_cell(const Point< dim > &p)
static std::array< unsigned int, 2 > standard_quad_vertex_to_line_vertex_index(const unsigned int vertex)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
RefinementCase operator|(const RefinementCase &r) const
UpdateFlags operator &(const UpdateFlags f1, const UpdateFlags f2)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
static constexpr unsigned int lines_per_cell
#define Assert(cond, exc)
Definition: exceptions.h:1466
static RefinementCase< dim - 1 > face_refinement_case(const RefinementCase< dim > &cell_refinement_case, const unsigned int face_no, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static constexpr unsigned int faces_per_cell
boost::integer_range< IncrementableType > iota_view
Definition: iota_view.h:46
static ::ExceptionBase & ExcInvalidCoordinate(double arg1)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:380
static unsigned int standard_to_real_line_vertex(const unsigned int vertex, const bool line_orientation=true)
static std::size_t memory_consumption()
Point< 3 > vertices[4]
static unsigned int standard_to_real_face_vertex(const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
static RefinementCase< dim > min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
static const std::array< unsigned int, vertices_per_cell > dx_to_deal
CacheUpdateFlags operator~(const CacheUpdateFlags f1)
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
ParameterHandler::OutputStyle operator|(const ParameterHandler::OutputStyle f1, const ParameterHandler::OutputStyle f2)
Definition: tensor.h:449
static constexpr unsigned int lines_per_face
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:379
T min(const T &t, const MPI_Comm &mpi_communicator)
static unsigned int n_subfaces(const internal::SubfaceCase< dim > &subface_case)
static RefinementCase< dim > min_cell_refinement_case_for_face_refinement(const RefinementCase< dim - 1 > &face_refinement_case, const unsigned int face_no, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static std::array< unsigned int, 2 > standard_hex_vertex_to_quad_vertex_index(const unsigned int vertex)
static ::ExceptionBase & ExcNotImplemented()
std::uint8_t value
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:561
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
static std::array< unsigned int, 2 > standard_hex_line_to_quad_line_index(const unsigned int line)
T max(const T &t, const MPI_Comm &mpi_communicator)
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()