Reference documentation for deal.II version Git d51799cb54 2020-09-28 09:22:08 +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 
1333  static unsigned int
1334  face_to_cell_vertices(const unsigned int face,
1335  const unsigned int vertex,
1336  const bool face_orientation = true,
1337  const bool face_flip = false,
1338  const bool face_rotation = false);
1339 
1354  static unsigned int
1355  face_to_cell_lines(const unsigned int face,
1356  const unsigned int line,
1357  const bool face_orientation = true,
1358  const bool face_flip = false,
1359  const bool face_rotation = false);
1360 
1367  static constexpr unsigned int vertices_per_face = 0;
1368 
1372  static constexpr unsigned int lines_per_face = 0;
1373 
1377  static constexpr unsigned int quads_per_face = 0;
1378 
1382  static constexpr unsigned int lines_per_cell = 0;
1383 
1387  static constexpr unsigned int quads_per_cell = 0;
1388 
1392  static constexpr unsigned int hexes_per_cell = 0;
1393 
1411  static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1412 
1426  static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1427 };
1428 
1429 
1430 
1954 template <int dim>
1955 struct GeometryInfo
1956 {
1964  static constexpr unsigned int max_children_per_cell = 1 << dim;
1965 
1969  static constexpr unsigned int faces_per_cell = 2 * dim;
1970 
1988  face_indices();
1989 
1997  static constexpr unsigned int max_children_per_face =
1998  GeometryInfo<dim - 1>::max_children_per_cell;
1999 
2003  static constexpr unsigned int vertices_per_cell = 1 << dim;
2004 
2021  vertex_indices();
2022 
2026  static constexpr unsigned int vertices_per_face =
2027  GeometryInfo<dim - 1>::vertices_per_cell;
2028 
2032  static constexpr unsigned int lines_per_face =
2033  GeometryInfo<dim - 1>::lines_per_cell;
2034 
2038  static constexpr unsigned int quads_per_face =
2039  GeometryInfo<dim - 1>::quads_per_cell;
2040 
2050  static constexpr unsigned int lines_per_cell =
2051  (2 * GeometryInfo<dim - 1>::lines_per_cell +
2052  GeometryInfo<dim - 1>::vertices_per_cell);
2053 
2061  static constexpr unsigned int quads_per_cell =
2062  (2 * GeometryInfo<dim - 1>::quads_per_cell +
2063  GeometryInfo<dim - 1>::lines_per_cell);
2064 
2068  static constexpr unsigned int hexes_per_cell =
2069  (2 * GeometryInfo<dim - 1>::hexes_per_cell +
2070  GeometryInfo<dim - 1>::quads_per_cell);
2071 
2089  static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
2090  internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal();
2091 
2105  static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
2106  internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal();
2107 
2118  static constexpr std::array<std::array<unsigned int, dim>, vertices_per_cell>
2119  vertex_to_face =
2120  internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face();
2121 
2126  static unsigned int
2127  n_children(const RefinementCase<dim> &refinement_case);
2128 
2133  static unsigned int
2134  n_subfaces(const internal::SubfaceCase<dim> &subface_case);
2135 
2145  static double
2146  subface_ratio(const internal::SubfaceCase<dim> &subface_case,
2147  const unsigned int subface_no);
2148 
2154  static RefinementCase<dim - 1>
2155  face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2156  const unsigned int face_no,
2157  const bool face_orientation = true,
2158  const bool face_flip = false,
2159  const bool face_rotation = false);
2160 
2166  static RefinementCase<dim>
2167  min_cell_refinement_case_for_face_refinement(
2168  const RefinementCase<dim - 1> &face_refinement_case,
2169  const unsigned int face_no,
2170  const bool face_orientation = true,
2171  const bool face_flip = false,
2172  const bool face_rotation = false);
2173 
2178  static RefinementCase<1>
2179  line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2180  const unsigned int line_no);
2181 
2186  static RefinementCase<dim>
2187  min_cell_refinement_case_for_line_refinement(const unsigned int line_no);
2188 
2235  static unsigned int
2236  child_cell_on_face(const RefinementCase<dim> & ref_case,
2237  const unsigned int face,
2238  const unsigned int subface,
2239  const bool face_orientation = true,
2240  const bool face_flip = false,
2241  const bool face_rotation = false,
2242  const RefinementCase<dim - 1> &face_refinement_case =
2244 
2258  static unsigned int
2259  line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
2260 
2281  static unsigned int
2282  face_to_cell_vertices(const unsigned int face,
2283  const unsigned int vertex,
2284  const bool face_orientation = true,
2285  const bool face_flip = false,
2286  const bool face_rotation = false);
2287 
2299  static unsigned int
2300  face_to_cell_lines(const unsigned int face,
2301  const unsigned int line,
2302  const bool face_orientation = true,
2303  const bool face_flip = false,
2304  const bool face_rotation = false);
2305 
2315  static unsigned int
2316  standard_to_real_face_vertex(const unsigned int vertex,
2317  const bool face_orientation = true,
2318  const bool face_flip = false,
2319  const bool face_rotation = false);
2320 
2330  static unsigned int
2331  real_to_standard_face_vertex(const unsigned int vertex,
2332  const bool face_orientation = true,
2333  const bool face_flip = false,
2334  const bool face_rotation = false);
2335 
2345  static unsigned int
2346  standard_to_real_face_line(const unsigned int line,
2347  const bool face_orientation = true,
2348  const bool face_flip = false,
2349  const bool face_rotation = false);
2350 
2356  static unsigned int
2357  standard_to_real_line_vertex(const unsigned int vertex,
2358  const bool line_orientation = true);
2359 
2367  static std::array<unsigned int, 2>
2368  standard_quad_vertex_to_line_vertex_index(const unsigned int vertex);
2369 
2377  static std::array<unsigned int, 2>
2378  standard_hex_vertex_to_quad_vertex_index(const unsigned int vertex);
2379 
2387  static std::array<unsigned int, 2>
2388  standard_hex_line_to_quad_line_index(const unsigned int line);
2389 
2399  static unsigned int
2400  real_to_standard_face_line(const unsigned int line,
2401  const bool face_orientation = true,
2402  const bool face_flip = false,
2403  const bool face_rotation = false);
2404 
2410  static Point<dim>
2411  unit_cell_vertex(const unsigned int vertex);
2412 
2422  static unsigned int
2423  child_cell_from_point(const Point<dim> &p);
2424 
2432  static Point<dim>
2433  cell_to_child_coordinates(const Point<dim> & p,
2434  const unsigned int child_index,
2435  const RefinementCase<dim> refine_case =
2437 
2443  static Point<dim>
2444  child_to_cell_coordinates(const Point<dim> & p,
2445  const unsigned int child_index,
2446  const RefinementCase<dim> refine_case =
2448 
2453  static bool
2454  is_inside_unit_cell(const Point<dim> &p);
2455 
2467  static bool
2468  is_inside_unit_cell(const Point<dim> &p, const double eps);
2469 
2474  static Point<dim>
2475  project_to_unit_cell(const Point<dim> &p);
2476 
2482  static double
2483  distance_to_unit_cell(const Point<dim> &p);
2484 
2489  static double
2490  d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2491 
2496  static Tensor<1, dim>
2497  d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2498 
2550  template <int spacedim>
2551  static void
2552  alternating_form_at_vertices
2553 #ifndef DEAL_II_CXX14_CONSTEXPR_BUG
2554  (const Point<spacedim> (&vertices)[vertices_per_cell],
2555  Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2556 #else
2557  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2558 #endif
2559 
2569  static constexpr std::array<unsigned int, faces_per_cell>
2570  unit_normal_direction =
2571  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction();
2572 
2589  static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2590  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation();
2591 
2602  static constexpr std::array<Tensor<1, dim>, faces_per_cell>
2603  unit_normal_vector =
2604  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_vector();
2605 
2619  static constexpr std::array<std::array<Tensor<1, dim>, dim - 1>,
2620  faces_per_cell>
2621  unit_tangential_vectors = internal::GeometryInfoHelper::Initializers<
2622  dim>::unit_tangential_vectors();
2623 
2629  static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2630  internal::GeometryInfoHelper::Initializers<dim>::opposite_face();
2631 
2632 
2636  DeclException1(ExcInvalidCoordinate,
2637  double,
2638  << "The coordinates must satisfy 0 <= x_i <= 1, "
2639  << "but here we have x_i=" << arg1);
2640 
2644  DeclException3(ExcInvalidSubface,
2645  int,
2646  int,
2647  int,
2648  << "RefinementCase<dim> " << arg1 << ": face " << arg2
2649  << " has no subface " << arg3);
2650 };
2651 
2652 
2653 
2654 #ifndef DOXYGEN
2655 
2656 
2657 /* -------------- declaration of explicit specializations ------------- */
2658 
2659 
2660 template <>
2663  const unsigned int i);
2664 template <>
2667  const unsigned int i);
2668 template <>
2671  const unsigned int i);
2672 
2673 
2674 
2675 /* -------------- inline functions ------------- */
2676 
2677 
2678 inline GeometryPrimitive::GeometryPrimitive(const Object object)
2679  : object(object)
2680 {}
2681 
2682 
2683 
2684 inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2685  : object(static_cast<Object>(object_dimension))
2686 {}
2687 
2688 
2689 inline GeometryPrimitive::operator unsigned int() const
2690 {
2691  return static_cast<unsigned int>(object);
2692 }
2693 
2694 
2695 
2696 namespace internal
2697 {
2698  template <int dim>
2699  inline SubfaceCase<dim>::SubfaceCase(
2700  const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2701  : value(subface_possibility)
2702  {}
2703 
2704 
2705  template <int dim>
2706  inline SubfaceCase<dim>::operator std::uint8_t() const
2707  {
2708  return value;
2709  }
2710 
2711 
2712 } // namespace internal
2713 
2714 
2715 template <int dim>
2716 inline RefinementCase<dim>
2717 RefinementCase<dim>::cut_axis(const unsigned int)
2718 {
2719  Assert(false, ExcInternalError());
2720  return static_cast<std::uint8_t>(-1);
2721 }
2722 
2723 
2724 template <>
2725 inline RefinementCase<1>
2726 RefinementCase<1>::cut_axis(const unsigned int i)
2727 {
2728  AssertIndexRange(i, 1);
2729 
2731  return options[i];
2732 }
2733 
2734 
2735 
2736 template <>
2737 inline RefinementCase<2>
2738 RefinementCase<2>::cut_axis(const unsigned int i)
2739 {
2740  AssertIndexRange(i, 2);
2741 
2744  return options[i];
2745 }
2746 
2747 
2748 
2749 template <>
2750 inline RefinementCase<3>
2751 RefinementCase<3>::cut_axis(const unsigned int i)
2752 {
2753  AssertIndexRange(i, 3);
2754 
2758  return options[i];
2759 }
2760 
2761 
2762 
2763 template <int dim>
2765  : value(RefinementPossibilities<dim>::no_refinement)
2766 {}
2767 
2768 
2769 
2770 template <int dim>
2772  const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2773  : value(refinement_case)
2774 {
2775  // check that only those bits of
2776  // the given argument are set that
2777  // make sense for a given space
2778  // dimension
2779  Assert((refinement_case &
2781  refinement_case,
2782  ExcInvalidRefinementCase(refinement_case));
2783 }
2784 
2785 
2786 
2787 template <int dim>
2788 inline RefinementCase<dim>::RefinementCase(const std::uint8_t refinement_case)
2789  : value(refinement_case)
2790 {
2791  // check that only those bits of
2792  // the given argument are set that
2793  // make sense for a given space
2794  // dimension
2795  Assert((refinement_case &
2797  refinement_case,
2798  ExcInvalidRefinementCase(refinement_case));
2799 }
2800 
2801 
2802 
2803 template <int dim>
2804 inline RefinementCase<dim>::operator std::uint8_t() const
2805 {
2806  return value;
2807 }
2808 
2809 
2810 
2811 template <int dim>
2812 inline RefinementCase<dim>
2814 {
2815  return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2816 }
2817 
2818 
2819 
2820 template <int dim>
2822  operator&(const RefinementCase<dim> &r) const
2823 {
2824  return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2825 }
2826 
2827 
2828 
2829 template <int dim>
2830 inline RefinementCase<dim>
2832 {
2833  return RefinementCase<dim>(static_cast<std::uint8_t>(
2835 }
2836 
2837 
2838 
2839 template <int dim>
2840 inline std::size_t
2842 {
2843  return sizeof(RefinementCase<dim>);
2844 }
2845 
2846 
2847 
2848 template <int dim>
2849 template <class Archive>
2850 inline void
2851 RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2852 {
2853  // serialization can't deal with bitfields, so copy from/to a full sized
2854  // std::uint8_t
2855  std::uint8_t uchar_value = value;
2856  ar & uchar_value;
2857  value = uchar_value;
2858 }
2859 
2860 
2861 
2862 template <>
2863 inline Point<1>
2864 GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2865 {
2867 
2868  return Point<1>(static_cast<double>(vertex));
2869 }
2870 
2871 
2872 
2873 template <>
2874 inline Point<2>
2875 GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2876 {
2878 
2879  return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2880 }
2881 
2882 
2883 
2884 template <>
2885 inline Point<3>
2886 GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2887 {
2889 
2890  return {static_cast<double>(vertex % 2),
2891  static_cast<double>(vertex / 2 % 2),
2892  static_cast<double>(vertex / 4)};
2893 }
2894 
2895 
2896 
2897 inline std::array<unsigned int, 0>
2899 {
2900  return {{}};
2901 }
2902 
2903 
2904 
2905 inline std::array<unsigned int, 1>
2907 {
2908  return {{0}};
2909 }
2910 
2911 
2912 
2913 template <int dim>
2916 {
2917  return {0U, faces_per_cell};
2918 }
2919 
2920 
2921 
2922 template <int dim>
2925 {
2926  return {0U, vertices_per_cell};
2927 }
2928 
2929 
2930 
2931 template <int dim>
2932 inline Point<dim>
2933 GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2934 {
2935  Assert(false, ExcNotImplemented());
2936 
2937  return Point<dim>();
2938 }
2939 
2940 
2941 
2942 template <>
2943 inline unsigned int
2945 {
2946  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2947 
2948  return (p[0] <= 0.5 ? 0 : 1);
2949 }
2950 
2951 
2952 
2953 template <>
2954 inline unsigned int
2956 {
2957  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2958  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2959 
2960  return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
2961 }
2962 
2963 
2964 
2965 template <>
2966 inline unsigned int
2968 {
2969  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2970  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2971  Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
2972 
2973  return (p[0] <= 0.5 ?
2974  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
2975  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
2976 }
2977 
2978 
2979 template <int dim>
2980 inline unsigned int
2982 {
2983  Assert(false, ExcNotImplemented());
2984 
2985  return 0;
2986 }
2987 
2988 
2989 
2990 template <>
2991 inline Point<1>
2993  const unsigned int child_index,
2994  const RefinementCase<1> refine_case)
2995 
2996 {
2997  AssertIndexRange(child_index, 2);
2999  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3000 
3001  return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
3002 }
3003 
3004 
3005 
3006 template <>
3007 inline Point<2>
3009  const unsigned int child_index,
3010  const RefinementCase<2> refine_case)
3011 
3012 {
3013  AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3014 
3015  Point<2> point = p;
3016  switch (refine_case)
3017  {
3019  point[0] *= 2.0;
3020  if (child_index == 1)
3021  point[0] -= 1.0;
3022  break;
3024  point[1] *= 2.0;
3025  if (child_index == 1)
3026  point[1] -= 1.0;
3027  break;
3029  point *= 2.0;
3030  point -= unit_cell_vertex(child_index);
3031  break;
3032  default:
3033  Assert(false, ExcInternalError());
3034  }
3035 
3036  return point;
3037 }
3038 
3039 
3040 
3041 template <>
3042 inline Point<3>
3044  const unsigned int child_index,
3045  const RefinementCase<3> refine_case)
3046 
3047 {
3048  AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3049 
3050  Point<3> point = p;
3051  // there might be a cleverer way to do
3052  // this, but since this function is called
3053  // in very few places for initialization
3054  // purposes only, I don't care at the
3055  // moment
3056  switch (refine_case)
3057  {
3059  point[0] *= 2.0;
3060  if (child_index == 1)
3061  point[0] -= 1.0;
3062  break;
3064  point[1] *= 2.0;
3065  if (child_index == 1)
3066  point[1] -= 1.0;
3067  break;
3069  point[2] *= 2.0;
3070  if (child_index == 1)
3071  point[2] -= 1.0;
3072  break;
3074  point[0] *= 2.0;
3075  point[1] *= 2.0;
3076  if (child_index % 2 == 1)
3077  point[0] -= 1.0;
3078  if (child_index / 2 == 1)
3079  point[1] -= 1.0;
3080  break;
3082  // careful, this is slightly
3083  // different from xy and yz due to
3084  // different internal numbering of
3085  // children!
3086  point[0] *= 2.0;
3087  point[2] *= 2.0;
3088  if (child_index / 2 == 1)
3089  point[0] -= 1.0;
3090  if (child_index % 2 == 1)
3091  point[2] -= 1.0;
3092  break;
3094  point[1] *= 2.0;
3095  point[2] *= 2.0;
3096  if (child_index % 2 == 1)
3097  point[1] -= 1.0;
3098  if (child_index / 2 == 1)
3099  point[2] -= 1.0;
3100  break;
3102  point *= 2.0;
3103  point -= unit_cell_vertex(child_index);
3104  break;
3105  default:
3106  Assert(false, ExcInternalError());
3107  }
3108 
3109  return point;
3110 }
3111 
3112 
3113 
3114 template <int dim>
3115 inline Point<dim>
3117  const Point<dim> & /*p*/,
3118  const unsigned int /*child_index*/,
3119  const RefinementCase<dim> /*refine_case*/)
3120 
3121 {
3122  Assert(false, ExcNotImplemented());
3123  return Point<dim>();
3124 }
3125 
3126 
3127 
3128 template <>
3129 inline Point<1>
3131  const unsigned int child_index,
3132  const RefinementCase<1> refine_case)
3133 
3134 {
3135  AssertIndexRange(child_index, 2);
3137  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3138 
3139  return (p + unit_cell_vertex(child_index)) * 0.5;
3140 }
3141 
3142 
3143 
3144 template <>
3145 inline Point<3>
3147  const unsigned int child_index,
3148  const RefinementCase<3> refine_case)
3149 
3150 {
3151  AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3152 
3153  Point<3> point = p;
3154  // there might be a cleverer way to do
3155  // this, but since this function is called
3156  // in very few places for initialization
3157  // purposes only, I don't care at the
3158  // moment
3159  switch (refine_case)
3160  {
3162  if (child_index == 1)
3163  point[0] += 1.0;
3164  point[0] *= 0.5;
3165  break;
3167  if (child_index == 1)
3168  point[1] += 1.0;
3169  point[1] *= 0.5;
3170  break;
3172  if (child_index == 1)
3173  point[2] += 1.0;
3174  point[2] *= 0.5;
3175  break;
3177  if (child_index % 2 == 1)
3178  point[0] += 1.0;
3179  if (child_index / 2 == 1)
3180  point[1] += 1.0;
3181  point[0] *= 0.5;
3182  point[1] *= 0.5;
3183  break;
3185  // careful, this is slightly
3186  // different from xy and yz due to
3187  // different internal numbering of
3188  // children!
3189  if (child_index / 2 == 1)
3190  point[0] += 1.0;
3191  if (child_index % 2 == 1)
3192  point[2] += 1.0;
3193  point[0] *= 0.5;
3194  point[2] *= 0.5;
3195  break;
3197  if (child_index % 2 == 1)
3198  point[1] += 1.0;
3199  if (child_index / 2 == 1)
3200  point[2] += 1.0;
3201  point[1] *= 0.5;
3202  point[2] *= 0.5;
3203  break;
3205  point += unit_cell_vertex(child_index);
3206  point *= 0.5;
3207  break;
3208  default:
3209  Assert(false, ExcInternalError());
3210  }
3211 
3212  return point;
3213 }
3214 
3215 
3216 
3217 template <>
3218 inline Point<2>
3220  const unsigned int child_index,
3221  const RefinementCase<2> refine_case)
3222 {
3223  AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3224 
3225  Point<2> point = p;
3226  switch (refine_case)
3227  {
3229  if (child_index == 1)
3230  point[0] += 1.0;
3231  point[0] *= 0.5;
3232  break;
3234  if (child_index == 1)
3235  point[1] += 1.0;
3236  point[1] *= 0.5;
3237  break;
3239  point += unit_cell_vertex(child_index);
3240  point *= 0.5;
3241  break;
3242  default:
3243  Assert(false, ExcInternalError());
3244  }
3245 
3246  return point;
3247 }
3248 
3249 
3250 
3251 template <int dim>
3252 inline Point<dim>
3254  const Point<dim> & /*p*/,
3255  const unsigned int /*child_index*/,
3256  const RefinementCase<dim> /*refine_case*/)
3257 {
3258  Assert(false, ExcNotImplemented());
3259  return Point<dim>();
3260 }
3261 
3262 
3263 
3264 template <int dim>
3265 inline bool
3267 {
3268  Assert(false, ExcNotImplemented());
3269  return false;
3270 }
3271 
3272 template <>
3273 inline bool
3275 {
3276  return (p[0] >= 0.) && (p[0] <= 1.);
3277 }
3278 
3279 
3280 
3281 template <>
3282 inline bool
3284 {
3285  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
3286 }
3287 
3288 
3289 
3290 template <>
3291 inline bool
3293 {
3294  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
3295  (p[2] >= 0.) && (p[2] <= 1.);
3296 }
3297 
3298 
3299 
3300 template <int dim>
3301 inline bool
3303 {
3304  Assert(false, ExcNotImplemented());
3305  return false;
3306 }
3307 
3308 template <>
3309 inline bool
3310 GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
3311 {
3312  return (p[0] >= -eps) && (p[0] <= 1. + eps);
3313 }
3314 
3315 
3316 
3317 template <>
3318 inline bool
3319 GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
3320 {
3321  const double l = -eps, u = 1 + eps;
3322  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
3323 }
3324 
3325 
3326 
3327 template <>
3328 inline bool
3329 GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
3330 {
3331  const double l = -eps, u = 1.0 + eps;
3332  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
3333  (p[2] >= l) && (p[2] <= u);
3334 }
3335 
3336 
3337 
3338 template <>
3339 inline unsigned int
3340 GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
3341  const unsigned int vertex)
3342 {
3343  (void)line;
3345  AssertIndexRange(vertex, 2);
3346 
3347  return vertex;
3348 }
3349 
3350 
3351 template <>
3352 inline unsigned int
3353 GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
3354  const unsigned int vertex)
3355 {
3356  constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
3357  return cell_vertices[line][vertex];
3358 }
3359 
3360 
3361 
3362 template <>
3363 inline unsigned int
3364 GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
3365  const unsigned int vertex)
3366 {
3368  AssertIndexRange(vertex, 2);
3369 
3370  constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3371  {1, 3},
3372  {0, 1},
3373  {2, 3},
3374  {4, 6}, // top face
3375  {5, 7},
3376  {4, 5},
3377  {6, 7},
3378  {0,
3379  4}, // connects of bottom
3380  {1, 5}, // top face
3381  {2, 6},
3382  {3, 7}};
3383  return vertices[line][vertex];
3384 }
3385 
3386 
3387 template <>
3388 inline unsigned int
3389 GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3390 {
3391  Assert(false, ExcNotImplemented());
3393 }
3394 
3395 template <>
3396 inline unsigned int
3397 GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3398  const bool face_orientation,
3399  const bool face_flip,
3400  const bool face_rotation)
3401 {
3403 
3404  // set up a table to make sure that
3405  // we handle non-standard faces correctly
3406  //
3407  // so set up a table that for each vertex (of
3408  // a quad in standard position) describes
3409  // which vertex to take
3410  //
3411  // first index: four vertices 0...3
3412  //
3413  // second index: face_orientation; 0:
3414  // opposite normal, 1: standard
3415  //
3416  // third index: face_flip; 0: standard, 1:
3417  // face rotated by 180 degrees
3418  //
3419  // forth index: face_rotation: 0: standard,
3420  // 1: face rotated by 90 degrees
3421 
3422  constexpr unsigned int vertex_translation[4][2][2][2] = {
3423  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3424  // face_rotation=false and true
3425  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3426  // face_rotation=false and true
3427  {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3428  // face_rotation=false and true
3429  {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3430  // face_rotation=false and true
3431 
3432  {{{2, 3}, // vertex 1 ...
3433  {1, 0}},
3434  {{1, 0}, {2, 3}}},
3435 
3436  {{{1, 0}, // vertex 2 ...
3437  {2, 3}},
3438  {{2, 3}, {1, 0}}},
3439 
3440  {{{3, 1}, // vertex 3 ...
3441  {0, 2}},
3442  {{3, 1}, {0, 2}}}};
3443 
3444  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3445 }
3446 
3447 
3448 
3449 template <int dim>
3450 inline unsigned int
3451 GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3452  const bool,
3453  const bool,
3454  const bool)
3455 {
3456  Assert(dim > 1, ExcImpossibleInDim(dim));
3458  return vertex;
3459 }
3460 
3461 template <int dim>
3462 inline unsigned int
3464 {
3465  constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3466  0, 2, 2, 4, 2, 4, 4, 8};
3467 
3468  return n_children[ref_case];
3469 }
3470 
3471 
3472 
3473 template <int dim>
3474 inline unsigned int
3476 {
3477  Assert(false, ExcNotImplemented());
3478  return 0;
3479 }
3480 
3481 template <>
3482 inline unsigned int
3484 {
3485  Assert(false, ExcImpossibleInDim(1));
3486  return 0;
3487 }
3488 
3489 template <>
3490 inline unsigned int
3492 {
3493  return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3494 }
3495 
3496 
3497 
3498 template <>
3499 inline unsigned int
3501 {
3502  const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3503  0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3504  return nsubs[subface_case];
3505 }
3506 
3507 
3508 
3509 template <int dim>
3510 inline double
3512  const unsigned int)
3513 {
3514  Assert(false, ExcNotImplemented());
3515  return 0.;
3516 }
3517 
3518 template <>
3519 inline double
3521  const unsigned int)
3522 {
3523  return 1;
3524 }
3525 
3526 
3527 template <>
3528 inline double
3530  const unsigned int)
3531 {
3532  double ratio = 1;
3533  switch (subface_case)
3534  {
3536  // Here, an
3537  // Assert(false,ExcInternalError())
3538  // would be the right
3539  // choice, but
3540  // unfortunately the
3541  // current function is
3542  // also called for faces
3543  // without children (see
3544  // tests/fe/mapping.cc).
3545  // Assert(false, ExcMessage("Face has no subfaces."));
3546  // Furthermore, assign
3547  // following value as
3548  // otherwise the
3549  // bits/volume_x tests
3550  // break
3552  break;
3554  ratio = 0.5;
3555  break;
3556  default:
3557  // there should be no
3558  // cases left
3559  Assert(false, ExcInternalError());
3560  break;
3561  }
3562 
3563  return ratio;
3564 }
3565 
3566 
3567 template <>
3568 inline double
3570  const unsigned int subface_no)
3571 {
3572  double ratio = 1;
3573  switch (subface_case)
3574  {
3576  // Here, an
3577  // Assert(false,ExcInternalError())
3578  // would be the right
3579  // choice, but
3580  // unfortunately the
3581  // current function is
3582  // also called for faces
3583  // without children (see
3584  // tests/bits/mesh_3d_16.cc). Add
3585  // following switch to
3586  // avoid diffs in
3587  // tests/bits/mesh_3d_16
3589  break;
3592  ratio = 0.5;
3593  break;
3597  ratio = 0.25;
3598  break;
3601  if (subface_no < 2)
3602  ratio = 0.25;
3603  else
3604  ratio = 0.5;
3605  break;
3608  if (subface_no == 0)
3609  ratio = 0.5;
3610  else
3611  ratio = 0.25;
3612  break;
3613  default:
3614  // there should be no
3615  // cases left
3616  Assert(false, ExcInternalError());
3617  break;
3618  }
3619 
3620  return ratio;
3621 }
3622 
3623 
3624 
3625 template <int dim>
3627  const RefinementCase<dim> &,
3628  const unsigned int,
3629  const bool,
3630  const bool,
3631  const bool)
3632 {
3633  Assert(false, ExcNotImplemented());
3634  return RefinementCase<dim - 1>::no_refinement;
3635 }
3636 
3637 template <>
3639  const RefinementCase<1> &,
3640  const unsigned int,
3641  const bool,
3642  const bool,
3643  const bool)
3644 {
3645  Assert(false, ExcImpossibleInDim(1));
3646 
3648 }
3649 
3650 
3651 template <>
3652 inline RefinementCase<1>
3654  const RefinementCase<2> &cell_refinement_case,
3655  const unsigned int face_no,
3656  const bool,
3657  const bool,
3658  const bool)
3659 {
3660  const unsigned int dim = 2;
3661  AssertIndexRange(cell_refinement_case,
3664 
3665  const RefinementCase<dim - 1>
3668  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3669  RefinementCase<dim - 1>::no_refinement},
3670 
3671  {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3672 
3673  {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3674 
3675  {RefinementCase<dim - 1>::cut_x, // cut_xy
3676  RefinementCase<dim - 1>::cut_x}};
3677 
3678  return ref_cases[cell_refinement_case][face_no / 2];
3679 }
3680 
3681 
3682 template <>
3683 inline RefinementCase<2>
3685  const RefinementCase<3> &cell_refinement_case,
3686  const unsigned int face_no,
3687  const bool face_orientation,
3688  const bool /*face_flip*/,
3689  const bool face_rotation)
3690 {
3691  const unsigned int dim = 3;
3692  AssertIndexRange(cell_refinement_case,
3695 
3696  const RefinementCase<dim - 1>
3699  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3700  RefinementCase<dim - 1>::no_refinement,
3701  RefinementCase<dim - 1>::no_refinement},
3702 
3703  {RefinementCase<dim - 1>::no_refinement, // cut_x
3704  RefinementCase<dim - 1>::cut_y,
3705  RefinementCase<dim - 1>::cut_x},
3706 
3707  {RefinementCase<dim - 1>::cut_x, // cut_y
3708  RefinementCase<dim - 1>::no_refinement,
3709  RefinementCase<dim - 1>::cut_y},
3710 
3711  {RefinementCase<dim - 1>::cut_x, // cut_xy
3712  RefinementCase<dim - 1>::cut_y,
3713  RefinementCase<dim - 1>::cut_xy},
3714 
3715  {RefinementCase<dim - 1>::cut_y, // cut_z
3716  RefinementCase<dim - 1>::cut_x,
3717  RefinementCase<dim - 1>::no_refinement},
3718 
3719  {RefinementCase<dim - 1>::cut_y, // cut_xz
3720  RefinementCase<dim - 1>::cut_xy,
3721  RefinementCase<dim - 1>::cut_x},
3722 
3723  {RefinementCase<dim - 1>::cut_xy, // cut_yz
3724  RefinementCase<dim - 1>::cut_x,
3725  RefinementCase<dim - 1>::cut_y},
3726 
3727  {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3728  RefinementCase<dim - 1>::cut_xy,
3729  RefinementCase<dim - 1>::cut_xy},
3730  };
3731 
3732  const RefinementCase<dim - 1> ref_case =
3733  ref_cases[cell_refinement_case][face_no / 2];
3734 
3735  const RefinementCase<dim - 1> flip[4] = {
3736  RefinementCase<dim - 1>::no_refinement,
3737  RefinementCase<dim - 1>::cut_y,
3738  RefinementCase<dim - 1>::cut_x,
3739  RefinementCase<dim - 1>::cut_xy};
3740 
3741  // correct the ref_case for face_orientation
3742  // and face_rotation. for face_orientation,
3743  // 'true' is the default value whereas for
3744  // face_rotation, 'false' is standard. If
3745  // <tt>face_rotation==face_orientation</tt>,
3746  // then one of them is non-standard and we
3747  // have to swap cut_x and cut_y, otherwise no
3748  // change is necessary. face_flip has no
3749  // influence. however, in order to keep the
3750  // interface consistent with other functions,
3751  // we still include it as an argument to this
3752  // function
3753  return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3754 }
3755 
3756 
3757 
3758 template <int dim>
3759 inline RefinementCase<1>
3761  const unsigned int)
3762 {
3763  Assert(false, ExcNotImplemented());
3765 }
3766 
3767 template <>
3768 inline RefinementCase<1>
3770  const RefinementCase<1> &cell_refinement_case,
3771  const unsigned int line_no)
3772 {
3773  (void)line_no;
3774  const unsigned int dim = 1;
3775  (void)dim;
3776  AssertIndexRange(cell_refinement_case,
3779 
3780  return cell_refinement_case;
3781 }
3782 
3783 
3784 template <>
3785 inline RefinementCase<1>
3787  const RefinementCase<2> &cell_refinement_case,
3788  const unsigned int line_no)
3789 {
3790  // Assertions are in face_refinement_case()
3791  return face_refinement_case(cell_refinement_case, line_no);
3792 }
3793 
3794 
3795 template <>
3796 inline RefinementCase<1>
3798  const RefinementCase<3> &cell_refinement_case,
3799  const unsigned int line_no)
3800 {
3801  const unsigned int dim = 3;
3802  AssertIndexRange(cell_refinement_case,
3805 
3806  // array indicating, which simple refine
3807  // case cuts a line in direction x, y or
3808  // z. For example, cut_y and everything
3809  // containing cut_y (cut_xy, cut_yz,
3810  // cut_xyz) cuts lines, which are in y
3811  // direction.
3812  const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3815 
3816  // order the direction of lines
3817  // 0->x, 1->y, 2->z
3818  const unsigned int direction[lines_per_cell] = {
3819  1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3820 
3821  return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3824 }
3825 
3826 
3827 
3828 template <int dim>
3829 inline RefinementCase<dim>
3831  const RefinementCase<dim - 1> &,
3832  const unsigned int,
3833  const bool,
3834  const bool,
3835  const bool)
3836 {
3837  Assert(false, ExcNotImplemented());
3838 
3840 }
3841 
3842 template <>
3843 inline RefinementCase<1>
3845  const RefinementCase<0> &,
3846  const unsigned int,
3847  const bool,
3848  const bool,
3849  const bool)
3850 {
3851  const unsigned int dim = 1;
3852  Assert(false, ExcImpossibleInDim(dim));
3853 
3855 }
3856 
3857 
3858 template <>
3859 inline RefinementCase<2>
3862  const unsigned int face_no,
3863  const bool,
3864  const bool,
3865  const bool)
3866 {
3867  const unsigned int dim = 2;
3868  AssertIndexRange(face_refinement_case,
3871 
3872  if (face_refinement_case == RefinementCase<dim>::cut_x)
3873  return (face_no / 2) ? RefinementCase<dim>::cut_x :
3875  else
3877 }
3878 
3879 
3880 template <>
3881 inline RefinementCase<3>
3883  const RefinementCase<2> &face_refinement_case,
3884  const unsigned int face_no,
3885  const bool face_orientation,
3886  const bool /*face_flip*/,
3887  const bool face_rotation)
3888 {
3889  const unsigned int dim = 3;
3890  AssertIndexRange(face_refinement_case,
3893 
3898 
3899  // correct the face_refinement_case for
3900  // face_orientation and face_rotation. for
3901  // face_orientation, 'true' is the default
3902  // value whereas for face_rotation, 'false'
3903  // is standard. If
3904  // <tt>face_rotation==face_orientation</tt>,
3905  // then one of them is non-standard and we
3906  // have to swap cut_x and cut_y, otherwise no
3907  // change is necessary. face_flip has no
3908  // influence. however, in order to keep the
3909  // interface consistent with other functions,
3910  // we still include it as an argument to this
3911  // function
3912  const RefinementCase<dim - 1> std_face_ref =
3913  (face_orientation == face_rotation) ? flip[face_refinement_case] :
3914  face_refinement_case;
3915 
3916  const RefinementCase<dim> face_to_cell[3][4] = {
3917  {RefinementCase<dim>::no_refinement, // faces 0 and 1
3918  RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3919  RefinementCase<dim>::cut_z,
3921 
3922  {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3923  // "exchanged on faces 2 and 3")
3924  RefinementCase<dim>::cut_z,
3927 
3928  {RefinementCase<dim>::no_refinement, // faces 4 and 5
3932 
3933  return face_to_cell[face_no / 2][std_face_ref];
3934 }
3935 
3936 
3937 
3938 template <int dim>
3939 inline RefinementCase<dim>
3941  const unsigned int)
3942 {
3943  Assert(false, ExcNotImplemented());
3944 
3946 }
3947 
3948 template <>
3949 inline RefinementCase<1>
3951  const unsigned int line_no)
3952 {
3953  (void)line_no;
3954  AssertIndexRange(line_no, 1);
3955 
3956  return RefinementCase<1>::cut_x;
3957 }
3958 
3959 
3960 template <>
3961 inline RefinementCase<2>
3963  const unsigned int line_no)
3964 {
3965  const unsigned int dim = 2;
3966  (void)dim;
3968 
3969  return (line_no / 2) ? RefinementCase<2>::cut_x : RefinementCase<2>::cut_y;
3970 }
3971 
3972 
3973 template <>
3974 inline RefinementCase<3>
3976  const unsigned int line_no)
3977 {
3978  const unsigned int dim = 3;
3980 
3981  const RefinementCase<dim> ref_cases[6] = {
3982  RefinementCase<dim>::cut_y, // lines 0 and 1
3983  RefinementCase<dim>::cut_x, // lines 2 and 3
3984  RefinementCase<dim>::cut_y, // lines 4 and 5
3985  RefinementCase<dim>::cut_x, // lines 6 and 7
3986  RefinementCase<dim>::cut_z, // lines 8 and 9
3987  RefinementCase<dim>::cut_z}; // lines 10 and 11
3988 
3989  return ref_cases[line_no / 2];
3990 }
3991 
3992 
3993 
3994 template <>
3995 inline unsigned int
3996 GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
3997  const bool face_orientation,
3998  const bool face_flip,
3999  const bool face_rotation)
4000 {
4002 
4003  // set up a table to make sure that
4004  // we handle non-standard faces correctly
4005  //
4006  // so set up a table that for each vertex (of
4007  // a quad in standard position) describes
4008  // which vertex to take
4009  //
4010  // first index: four vertices 0...3
4011  //
4012  // second index: face_orientation; 0:
4013  // opposite normal, 1: standard
4014  //
4015  // third index: face_flip; 0: standard, 1:
4016  // face rotated by 180 degrees
4017  //
4018  // forth index: face_rotation: 0: standard,
4019  // 1: face rotated by 90 degrees
4020 
4021  const unsigned int vertex_translation[4][2][2][2] = {
4022  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
4023  // face_rotation=false and true
4024  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
4025  // face_rotation=false and true
4026  {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
4027  // face_rotation=false and true
4028  {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
4029  // face_rotation=false and true
4030 
4031  {{{2, 3}, // vertex 1 ...
4032  {1, 0}},
4033  {{1, 3}, {2, 0}}},
4034 
4035  {{{1, 0}, // vertex 2 ...
4036  {2, 3}},
4037  {{2, 0}, {1, 3}}},
4038 
4039  {{{3, 1}, // vertex 3 ...
4040  {0, 2}},
4041  {{3, 2}, {0, 1}}}};
4042 
4043  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
4044 }
4045 
4046 
4047 
4048 template <int dim>
4049 inline unsigned int
4050 GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
4051  const bool,
4052  const bool,
4053  const bool)
4054 {
4055  Assert(dim > 1, ExcImpossibleInDim(dim));
4057  return vertex;
4058 }
4059 
4060 
4061 
4062 template <>
4063 inline unsigned int
4064 GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
4065  const bool face_orientation,
4066  const bool face_flip,
4067  const bool face_rotation)
4068 {
4070 
4071 
4072  // make sure we handle
4073  // non-standard faces correctly
4074  //
4075  // so set up a table that for each line (of a
4076  // quad) describes which line to take
4077  //
4078  // first index: four lines 0...3
4079  //
4080  // second index: face_orientation; 0:
4081  // opposite normal, 1: standard
4082  //
4083  // third index: face_flip; 0: standard, 1:
4084  // face rotated by 180 degrees
4085  //
4086  // forth index: face_rotation: 0: standard,
4087  // 1: face rotated by 90 degrees
4088 
4089  const unsigned int line_translation[4][2][2][2] = {
4090  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4091  // face_rotation=false and true
4092  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4093  // face_rotation=false and true
4094  {{0, 3}, // line 0, face_orientation=true, face_flip=false,
4095  // face_rotation=false and true
4096  {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
4097  // face_rotation=false and true
4098 
4099  {{{3, 1}, // line 1 ...
4100  {2, 0}},
4101  {{1, 2}, {0, 3}}},
4102 
4103  {{{0, 3}, // line 2 ...
4104  {1, 2}},
4105  {{2, 0}, {3, 1}}},
4106 
4107  {{{1, 2}, // line 3 ...
4108  {0, 3}},
4109  {{3, 1}, {2, 0}}}};
4110 
4111  return line_translation[line][face_orientation][face_flip][face_rotation];
4112 }
4113 
4114 
4115 
4116 template <int dim>
4117 inline unsigned int
4118 GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
4119  const bool,
4120  const bool,
4121  const bool)
4122 {
4123  Assert(false, ExcNotImplemented());
4124  return line;
4125 }
4126 
4127 
4128 
4129 template <>
4130 inline unsigned int
4131 GeometryInfo<2>::standard_to_real_line_vertex(const unsigned int vertex,
4132  const bool line_orientation)
4133 {
4134  return line_orientation ? vertex : (1 - vertex);
4135 }
4136 
4137 
4138 
4139 template <int dim>
4140 inline unsigned int
4141 GeometryInfo<dim>::standard_to_real_line_vertex(const unsigned int vertex,
4142  const bool)
4143 {
4144  Assert(false, ExcNotImplemented());
4145  return vertex;
4146 }
4147 
4148 
4149 
4150 template <>
4151 inline std::array<unsigned int, 2>
4153  const unsigned int vertex)
4154 {
4155  return {{vertex % 2, vertex / 2}};
4156 }
4157 
4158 
4159 
4160 template <int dim>
4161 inline std::array<unsigned int, 2>
4163  const unsigned int vertex)
4164 {
4165  Assert(false, ExcNotImplemented());
4166  (void)vertex;
4167  return {{0, 0}};
4168 }
4169 
4170 
4171 
4172 template <>
4173 inline std::array<unsigned int, 2>
4175 {
4176  // set up a table that for each
4177  // line describes a) from which
4178  // quad to take it, b) which line
4179  // therein it is if the face is
4180  // oriented correctly
4181  static const unsigned int lookup_table[GeometryInfo<3>::lines_per_cell][2] = {
4182  {4, 0}, // take first four lines from bottom face
4183  {4, 1},
4184  {4, 2},
4185  {4, 3},
4186 
4187  {5, 0}, // second four lines from top face
4188  {5, 1},
4189  {5, 2},
4190  {5, 3},
4191 
4192  {0, 0}, // the rest randomly
4193  {1, 0},
4194  {0, 1},
4195  {1, 1}};
4196 
4197  return {{lookup_table[i][0], lookup_table[i][1]}};
4198 }
4199 
4200 
4201 
4202 template <int dim>
4203 inline std::array<unsigned int, 2>
4205 {
4206  Assert(false, ExcNotImplemented());
4207  (void)line;
4208  return {{0, 0}};
4209 }
4210 
4211 
4212 
4213 template <>
4214 inline std::array<unsigned int, 2>
4216  const unsigned int vertex)
4217 {
4218  // get the corner indices by asking either the bottom or the top face for its
4219  // vertices. handle non-standard faces by calling the vertex reordering
4220  // function from GeometryInfo
4221 
4222  // bottom face (4) for first four vertices, top face (5) for the rest
4223  return {{4 + vertex / 4, vertex % 4}};
4224 }
4225 
4226 
4227 
4228 template <int dim>
4229 inline std::array<unsigned int, 2>
4231  const unsigned int vertex)
4232 {
4233  Assert(false, ExcNotImplemented());
4234  (void)vertex;
4235  return {{0, 0}};
4236 }
4237 
4238 
4239 
4240 template <>
4241 inline unsigned int
4242 GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
4243  const bool face_orientation,
4244  const bool face_flip,
4245  const bool face_rotation)
4246 {
4248 
4249 
4250  // make sure we handle
4251  // non-standard faces correctly
4252  //
4253  // so set up a table that for each line (of a
4254  // quad) describes which line to take
4255  //
4256  // first index: four lines 0...3
4257  //
4258  // second index: face_orientation; 0:
4259  // opposite normal, 1: standard
4260  //
4261  // third index: face_flip; 0: standard, 1:
4262  // face rotated by 180 degrees
4263  //
4264  // forth index: face_rotation: 0: standard,
4265  // 1: face rotated by 90 degrees
4266 
4267  const unsigned int line_translation[4][2][2][2] = {
4268  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4269  // face_rotation=false and true
4270  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4271  // face_rotation=false and true
4272  {{0, 2}, // line 0, face_orientation=true, face_flip=false,
4273  // face_rotation=false and true
4274  {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
4275  // face_rotation=false and true
4276 
4277  {{{3, 1}, // line 1 ...
4278  {2, 0}},
4279  {{1, 3}, {0, 2}}},
4280 
4281  {{{0, 3}, // line 2 ...
4282  {1, 2}},
4283  {{2, 1}, {3, 0}}},
4284 
4285  {{{1, 2}, // line 3 ...
4286  {0, 3}},
4287  {{3, 0}, {2, 1}}}};
4288 
4289  return line_translation[line][face_orientation][face_flip][face_rotation];
4290 }
4291 
4292 
4293 
4294 template <int dim>
4295 inline unsigned int
4296 GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
4297  const bool,
4298  const bool,
4299  const bool)
4300 {
4301  Assert(false, ExcNotImplemented());
4302  return line;
4303 }
4304 
4305 
4306 
4307 template <>
4308 inline unsigned int
4310  const unsigned int face,
4311  const unsigned int subface,
4312  const bool,
4313  const bool,
4314  const bool,
4315  const RefinementCase<0> &)
4316 {
4317  (void)subface;
4320 
4321  return face;
4322 }
4323 
4324 
4325 
4326 template <>
4327 inline unsigned int
4329  const unsigned int face,
4330  const unsigned int subface,
4331  const bool /*face_orientation*/,
4332  const bool face_flip,
4333  const bool /*face_rotation*/,
4334  const RefinementCase<1> &)
4335 {
4338 
4339  // always return the child adjacent to the specified
4340  // subface. if the face of a cell is not refined, don't
4341  // throw an assertion but deliver the child adjacent to
4342  // the face nevertheless, i.e. deliver the child of
4343  // this cell adjacent to the subface of a possibly
4344  // refined neighbor. this simplifies setting neighbor
4345  // information in execute_refinement.
4346  const unsigned int
4348  [max_children_per_face] = {
4349  {
4350  // Normal orientation (face_flip = false)
4351  {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
4352  {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
4353  {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
4354  },
4355  {
4356  // Flipped orientation (face_flip = true)
4357  {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
4358  {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
4359  {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
4360  }};
4361 
4362  return subcells[face_flip][ref_case - 1][face][subface];
4363 }
4364 
4365 
4366 
4367 template <>
4368 inline unsigned int
4370  const unsigned int face,
4371  const unsigned int subface,
4372  const bool face_orientation,
4373  const bool face_flip,
4374  const bool face_rotation,
4375  const RefinementCase<2> &face_ref_case)
4376 {
4377  const unsigned int dim = 3;
4378 
4380  ExcMessage("Cell has no children."));
4382  if (!(subface == 0 &&
4383  face_ref_case == RefinementCase<dim - 1>::no_refinement))
4384  {
4385  AssertIndexRange(subface,
4386  GeometryInfo<dim - 1>::n_children(face_ref_case));
4387  }
4388 
4389  // invalid number used for invalid cases,
4390  // e.g. when the children are more refined at
4391  // a given face than the face itself
4392  const unsigned int e = numbers::invalid_unsigned_int;
4393 
4394  // the whole process of finding a child cell
4395  // at a given subface considering the
4396  // possibly anisotropic refinement cases of
4397  // the cell and the face as well as
4398  // orientation, flip and rotation of the face
4399  // is quite complicated. thus, we break it
4400  // down into several steps.
4401 
4402  // first step: convert the given face refine
4403  // case to a face refine case concerning the
4404  // face in standard orientation (, flip and
4405  // rotation). This only affects cut_x and
4406  // cut_y
4407  const RefinementCase<dim - 1> flip[4] = {
4408  RefinementCase<dim - 1>::no_refinement,
4409  RefinementCase<dim - 1>::cut_y,
4410  RefinementCase<dim - 1>::cut_x,
4411  RefinementCase<dim - 1>::cut_xy};
4412  // for face_orientation, 'true' is the
4413  // default value whereas for face_rotation,
4414  // 'false' is standard. If
4415  // <tt>face_rotation==face_orientation</tt>,
4416  // then one of them is non-standard and we
4417  // have to swap cut_x and cut_y, otherwise no
4418  // change is necessary.
4419  const RefinementCase<dim - 1> std_face_ref =
4420  (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4421 
4422  // second step: convert the given subface
4423  // index to the one for a standard face
4424  // respecting face_orientation, face_flip and
4425  // face_rotation
4426 
4427  // first index: face_ref_case
4428  // second index: face_orientation
4429  // third index: face_flip
4430  // forth index: face_rotation
4431  // fifth index: subface index
4432  const unsigned int subface_exchange[4][2][2][2][4] = {
4433  // no_refinement (subface 0 stays 0,
4434  // all others are invalid)
4435  {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4436  {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4437  // cut_x (here, if the face is only
4438  // rotated OR only falsely oriented,
4439  // then subface 0 of the non-standard
4440  // face does NOT correspond to one of
4441  // the subfaces of a standard
4442  // face. Thus we indicate the subface
4443  // which is located at the lower left
4444  // corner (the origin of the face's
4445  // local coordinate system) with
4446  // '0'. The rest of this issue is
4447  // taken care of using the above
4448  // conversion to a 'standard face
4449  // refine case')
4450  {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4451  {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4452  // cut_y (the same applies as for
4453  // cut_x)
4454  {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4455  {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4456  // cut_xyz: this information is
4457  // identical to the information
4458  // returned by
4459  // GeometryInfo<3>::real_to_standard_face_vertex()
4460  {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4461  // face_rotation=false, subfaces 0,1,2,3
4462  {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4463  // face_rotation=true, subfaces 0,1,2,3
4464  {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4465  // face_rotation=false, subfaces 0,1,2,3
4466  {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4467  // face_rotation=true, subfaces 0,1,2,3
4468  {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4469  // face_rotation=false, subfaces 0,1,2,3
4470  {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4471  // face_rotation=true, subfaces 0,1,2,3
4472  {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4473  // face_rotation=false, subfaces 0,1,2,3
4474  {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4475  // face_rotation=true, subfaces 0,1,2,3
4476 
4477  const unsigned int std_subface =
4478  subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4479  [subface];
4480  Assert(std_subface != e, ExcInternalError());
4481 
4482  // third step: these are the children, which
4483  // can be found at the given subfaces of an
4484  // isotropically refined (standard) face
4485  //
4486  // first index: (refinement_case-1)
4487  // second index: face_index
4488  // third index: subface_index (isotropic refinement)
4489  const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4490  [max_children_per_face] = {
4491  // cut_x
4492  {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4493  {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4494  {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4495  {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4496  {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4497  {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4498  // cut_y
4499  {{0, 1, 0, 1},
4500  {0, 1, 0, 1},
4501  {0, 0, 0, 0},
4502  {1, 1, 1, 1},
4503  {0, 0, 1, 1},
4504  {0, 0, 1, 1}},
4505  // cut_xy
4506  {{0, 2, 0, 2},
4507  {1, 3, 1, 3},
4508  {0, 0, 1, 1},
4509  {2, 2, 3, 3},
4510  {0, 1, 2, 3},
4511  {0, 1, 2, 3}},
4512  // cut_z
4513  {{0, 0, 1, 1},
4514  {0, 0, 1, 1},
4515  {0, 1, 0, 1},
4516  {0, 1, 0, 1},
4517  {0, 0, 0, 0},
4518  {1, 1, 1, 1}},
4519  // cut_xz
4520  {{0, 0, 1, 1},
4521  {2, 2, 3, 3},
4522  {0, 1, 2, 3},
4523  {0, 1, 2, 3},
4524  {0, 2, 0, 2},
4525  {1, 3, 1, 3}},
4526  // cut_yz
4527  {{0, 1, 2, 3},
4528  {0, 1, 2, 3},
4529  {0, 2, 0, 2},
4530  {1, 3, 1, 3},
4531  {0, 0, 1, 1},
4532  {2, 2, 3, 3}},
4533  // cut_xyz
4534  {{0, 2, 4, 6},
4535  {1, 3, 5, 7},
4536  {0, 4, 1, 5},
4537  {2, 6, 3, 7},
4538  {0, 1, 2, 3},
4539  {4, 5, 6, 7}}};
4540 
4541  // forth step: check, whether the given face
4542  // refine case is valid for the given cell
4543  // refine case. this is the case, if the
4544  // given face refine case is at least as
4545  // refined as the face is for the given cell
4546  // refine case
4547 
4548  // note, that we are considering standard
4549  // face refinement cases here and thus must
4550  // not pass the given orientation, flip and
4551  // rotation flags
4552  if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4553  face_refinement_case(ref_case, face))
4554  {
4555  // all is fine. for anisotropic face
4556  // refine cases, select one of the
4557  // isotropic subfaces which neighbors the
4558  // same child
4559 
4560  // first index: (standard) face refine case
4561  // second index: subface index
4562  const unsigned int equivalent_iso_subface[4][4] = {
4563  {0, e, e, e}, // no_refinement
4564  {0, 3, e, e}, // cut_x
4565  {0, 3, e, e}, // cut_y
4566  {0, 1, 2, 3}}; // cut_xy
4567 
4568  const unsigned int equ_std_subface =
4569  equivalent_iso_subface[std_face_ref][std_subface];
4570  Assert(equ_std_subface != e, ExcInternalError());
4571 
4572  return iso_children[ref_case - 1][face][equ_std_subface];
4573  }
4574  else
4575  {
4576  // the face_ref_case was too coarse,
4577  // throw an error
4578  Assert(false,
4579  ExcMessage("The face RefineCase is too coarse "
4580  "for the given cell RefineCase."));
4581  }
4582  // we only get here in case of an error
4583  return e;
4584 }
4585 
4586 
4587 
4588 template <>
4589 inline unsigned int
4591  const unsigned int,
4592  const unsigned int,
4593  const bool,
4594  const bool,
4595  const bool,
4596  const RefinementCase<3> &)
4597 {
4598  Assert(false, ExcNotImplemented());
4600 }
4601 
4602 
4603 
4604 template <>
4605 inline unsigned int
4606 GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4607  const unsigned int line,
4608  const bool,
4609  const bool,
4610  const bool)
4611 {
4612  (void)line;
4615 
4616  // The face is a line itself.
4617  return face;
4618 }
4619 
4620 
4621 
4622 template <>
4623 inline unsigned int
4624 GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4625  const unsigned int line,
4626  const bool face_orientation,
4627  const bool face_flip,
4628  const bool face_rotation)
4629 {
4632 
4633  const unsigned lines[faces_per_cell][lines_per_face] = {
4634  {8, 10, 0, 4}, // left face
4635  {9, 11, 1, 5}, // right face
4636  {2, 6, 8, 9}, // front face
4637  {3, 7, 10, 11}, // back face
4638  {0, 1, 2, 3}, // bottom face
4639  {4, 5, 6, 7}}; // top face
4640  return lines[face][real_to_standard_face_line(
4641  line, face_orientation, face_flip, face_rotation)];
4642 }
4643 
4644 
4645 
4646 inline unsigned int
4647 GeometryInfo<0>::face_to_cell_lines(const unsigned int,
4648  const unsigned int,
4649  const bool,
4650  const bool,
4651  const bool)
4652 {
4653  Assert(false, ExcNotImplemented());
4655 }
4656 
4657 
4658 
4659 template <int dim>
4660 inline unsigned int
4661 GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4662  const unsigned int,
4663  const bool,
4664  const bool,
4665  const bool)
4666 {
4667  Assert(false, ExcNotImplemented());
4669 }
4670 
4671 
4672 
4673 template <int dim>
4674 inline unsigned int
4675 GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4676  const unsigned int vertex,
4677  const bool face_orientation,
4678  const bool face_flip,
4679  const bool face_rotation)
4680 {
4682  face,
4683  vertex,
4684  face_orientation,
4685  face_flip,
4686  face_rotation);
4687 }
4688 
4689 
4690 
4691 inline unsigned int
4692 GeometryInfo<0>::face_to_cell_vertices(const unsigned int,
4693  const unsigned int,
4694  const bool,
4695  const bool,
4696  const bool)
4697 {
4698  Assert(false, ExcNotImplemented());
4700 }
4701 
4702 
4703 
4704 template <int dim>
4705 inline Point<dim>
4707 {
4708  Point<dim> p = q;
4709  for (unsigned int i = 0; i < dim; i++)
4710  if (p[i] < 0.)
4711  p[i] = 0.;
4712  else if (p[i] > 1.)
4713  p[i] = 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  if ((-p[i]) > result)
4728  result = -p[i];
4729  else if ((p[i] - 1.) > result)
4730  result = (p[i] - 1.);
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)
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:1636
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:2436
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:1411
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()