Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
geometry_info.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2019 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 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/point.h>
24 
25 #include <cstdint>
26 
27 
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 #ifndef DOXYGEN
32 namespace internal
33 {
34  namespace GeometryInfoHelper
35  {
36  // A struct that holds the values for all the arrays we want to initialize
37  // in GeometryInfo
38  template <int dim>
39  struct Initializers;
40 
41  template <>
42  struct Initializers<1>
43  {
44  static constexpr std::array<unsigned int, 2>
45  ucd_to_deal()
46  {
47  return {{0, 1}};
48  }
49 
50  static constexpr std::array<unsigned int, 2>
51  unit_normal_direction()
52  {
53  return {{0, 0}};
54  }
55 
56  static constexpr std::array<int, 2>
57  unit_normal_orientation()
58  {
59  return {{-1, 1}};
60  }
61 
62  static constexpr std::array<Tensor<1, 1>, 2>
63  unit_normal_vector()
64  {
65  return {{Tensor<1, 1>{{-1}}, Tensor<1, 1>{{1}}}};
66  }
67 
68  static constexpr std::array<unsigned int, 2>
69  opposite_face()
70  {
71  return {{1, 0}};
72  }
73 
74  static constexpr std::array<unsigned int, 2>
75  dx_to_deal()
76  {
77  return {{0, 1}};
78  }
79 
80  static constexpr std::array<std::array<unsigned int, 1>, 2>
81  vertex_to_face()
82  {
83  return {{{{0}}, {{1}}}};
84  }
85  };
86 
87  template <>
88  struct Initializers<2>
89  {
90  static constexpr std::array<unsigned int, 4>
91  ucd_to_deal()
92  {
93  return {{0, 1, 3, 2}};
94  }
95 
96  static constexpr std::array<unsigned int, 4>
97  unit_normal_direction()
98  {
99  return {{0, 0, 1, 1}};
100  }
101 
102  static constexpr std::array<int, 4>
103  unit_normal_orientation()
104  {
105  return {{-1, 1, -1, 1}};
106  }
107 
108  static constexpr std::array<Tensor<1, 2>, 4>
109  unit_normal_vector()
110  {
111  return {{Tensor<1, 2>{{-1., 0.}},
112  Tensor<1, 2>{{1., 0.}},
113  Tensor<1, 2>{{0., -1.}},
114  Tensor<1, 2>{{0., 1.}}}};
115  }
116 
117  static constexpr std::array<unsigned int, 4>
118  opposite_face()
119  {
120  return {{1, 0, 3, 2}};
121  }
122 
123  static constexpr std::array<unsigned int, 4>
124  dx_to_deal()
125  {
126  return {{0, 2, 1, 3}};
127  }
128 
129  static constexpr std::array<std::array<unsigned int, 2>, 4>
130  vertex_to_face()
131  {
132  return {{{{0, 2}}, {{1, 2}}, {{0, 3}}, {{1, 3}}}};
133  }
134  };
135 
136  template <>
137  struct Initializers<3>
138  {
139  static constexpr std::array<unsigned int, 8>
140  ucd_to_deal()
141  {
142  return {{0, 1, 5, 4, 2, 3, 7, 6}};
143  }
144 
145  static constexpr std::array<unsigned int, 6>
146  unit_normal_direction()
147  {
148  return {{0, 0, 1, 1, 2, 2}};
149  }
150 
151  static constexpr std::array<int, 6>
152  unit_normal_orientation()
153  {
154  return {{-1, 1, -1, 1, -1, 1}};
155  }
156 
157  static constexpr std::array<Tensor<1, 3>, 6>
158  unit_normal_vector()
159  {
160  return {{Tensor<1, 3>{{-1, 0, 0}},
161  Tensor<1, 3>{{1, 0, 0}},
162  Tensor<1, 3>{{0, -1, 0}},
163  Tensor<1, 3>{{0, 1, 0}},
164  Tensor<1, 3>{{0, 0, -1}},
165  Tensor<1, 3>{{0, 0, 1}}}};
166  }
167 
168  static constexpr std::array<unsigned int, 6>
169  opposite_face()
170  {
171  return {{1, 0, 3, 2, 5, 4}};
172  }
173 
174  static constexpr std::array<unsigned int, 8>
175  dx_to_deal()
176  {
177  return {{0, 4, 2, 6, 1, 5, 3, 7}};
178  }
179 
180  static constexpr std::array<std::array<unsigned int, 3>, 8>
181  vertex_to_face()
182  {
183  return {{{{0, 2, 4}},
184  {{1, 2, 4}},
185  {{0, 3, 4}},
186  {{1, 3, 4}},
187  {{0, 2, 5}},
188  {{1, 2, 5}},
189  {{0, 3, 5}},
190  {{1, 3, 5}}}};
191  }
192  };
193 
194  template <>
195  struct Initializers<4>
196  {
197  static constexpr std::array<unsigned int, 16>
198  ucd_to_deal()
199  {
215  numbers::invalid_unsigned_int}};
216  }
217 
218  static constexpr std::array<unsigned int, 8>
219  unit_normal_direction()
220  {
221  return {{0, 0, 1, 1, 2, 2, 3, 3}};
222  }
223 
224  static constexpr std::array<int, 8>
225  unit_normal_orientation()
226  {
227  return {{-1, 1, -1, 1, -1, 1, -1, 1}};
228  }
229 
230  static constexpr std::array<Tensor<1, 4>, 8>
231  unit_normal_vector()
232  {
233  return {{Tensor<1, 4>{{-1, 0, 0, 0}},
234  Tensor<1, 4>{{1, 0, 0, 0}},
235  Tensor<1, 4>{{0, -1, 0, 0}},
236  Tensor<1, 4>{{0, 1, 0, 0}},
237  Tensor<1, 4>{{0, 0, -1, 0}},
238  Tensor<1, 4>{{0, 0, 1, 0}},
239  Tensor<1, 4>{{0, 0, 0, -1}},
240  Tensor<1, 4>{{0, 0, 0, 1}}}};
241  }
242 
243  static constexpr std::array<unsigned int, 8>
244  opposite_face()
245  {
246  return {{1, 0, 3, 2, 5, 4, 7, 6}};
247  }
248 
249  static constexpr std::array<unsigned int, 16>
250  dx_to_deal()
251  {
267  numbers::invalid_unsigned_int}};
268  }
269 
270  static constexpr std::array<std::array<unsigned int, 4>, 16>
271  vertex_to_face()
272  {
276  numbers::invalid_unsigned_int}},
280  numbers::invalid_unsigned_int}},
284  numbers::invalid_unsigned_int}},
288  numbers::invalid_unsigned_int}},
292  numbers::invalid_unsigned_int}},
296  numbers::invalid_unsigned_int}},
300  numbers::invalid_unsigned_int}},
304  numbers::invalid_unsigned_int}},
308  numbers::invalid_unsigned_int}},
312  numbers::invalid_unsigned_int}},
316  numbers::invalid_unsigned_int}},
320  numbers::invalid_unsigned_int}},
324  numbers::invalid_unsigned_int}},
328  numbers::invalid_unsigned_int}},
332  numbers::invalid_unsigned_int}},
336  numbers::invalid_unsigned_int}}}};
337  }
338  };
339  } // namespace GeometryInfoHelper
340 } // namespace internal
341 #endif // DOXYGEN
342 
343 
362 {
363 public:
370  enum Object
371  {
375  vertex = 0,
379  line = 1,
383  quad = 2,
387  hex = 3
388  };
389 
394  GeometryPrimitive(const Object object);
395 
401  GeometryPrimitive(const unsigned int object_dimension);
402 
407  operator unsigned int() const;
408 
409 private:
414 };
415 
416 
417 
431 template <int dim>
433 {
469  {
473  no_refinement = 0,
474 
486  isotropic_refinement = 0xFF
487  };
488 };
489 
490 
491 
502 template <>
504 {
540  {
544  no_refinement = 0,
548  cut_x = 1,
552  isotropic_refinement = cut_x
553  };
554 };
555 
556 
557 
569 template <>
571 {
607  {
611  no_refinement = 0,
615  cut_x = 1,
619  cut_y = 2,
623  cut_xy = cut_x | cut_y,
624 
628  isotropic_refinement = cut_xy
629  };
630 };
631 
632 
633 
645 template <>
647 {
683  {
687  no_refinement = 0,
691  cut_x = 1,
695  cut_y = 2,
699  cut_xy = cut_x | cut_y,
703  cut_z = 4,
707  cut_xz = cut_x | cut_z,
711  cut_yz = cut_y | cut_z,
715  cut_xyz = cut_x | cut_y | cut_z,
716 
720  isotropic_refinement = cut_xyz
721  };
722 };
723 
724 
725 
737 template <int dim>
739 {
740 public:
744  RefinementCase();
745 
751  const typename RefinementPossibilities<dim>::Possibilities refinement_case);
752 
758  explicit RefinementCase(const std::uint8_t refinement_case);
759 
771  operator std::uint8_t() const;
772 
778  operator|(const RefinementCase &r) const;
779 
784  RefinementCase operator&(const RefinementCase &r) const;
785 
794  operator~() const;
795 
796 
802  static RefinementCase
803  cut_axis(const unsigned int i);
804 
808  static std::size_t
809  memory_consumption();
810 
815  template <class Archive>
816  void
817  serialize(Archive &ar, const unsigned int version);
818 
823  ExcInvalidRefinementCase,
824  int,
825  << "The refinement flags given (" << arg1
826  << ") contain set bits that do not "
827  << "make sense for the space dimension of the object to which they are applied.");
828 
829 private:
834  std::uint8_t value : (dim > 0 ? dim : 1);
835 };
836 
837 
838 namespace internal
839 {
860  template <int dim>
862  {
867  {
871  case_none = 0,
872 
876  case_isotropic = static_cast<std::uint8_t>(-1)
877  };
878  };
879 
880 
890  template <>
892  {
899  {
903  case_none = 0,
904 
908  case_isotropic = case_none
909  };
910  };
911 
912 
913 
924  template <>
926  {
933  {
937  case_none = 0,
938 
942  case_isotropic = case_none
943  };
944  };
945 
946 
947 
959  template <>
961  {
969  {
973  case_none = 0,
977  case_x = 1,
981  case_isotropic = case_x
982  };
983  };
984 
985 
986 
1079  template <>
1081  {
1089  {
1090  case_none = 0,
1091  case_x = 1,
1092  case_x1y = 2,
1093  case_x2y = 3,
1094  case_x1y2y = 4,
1095  case_y = 5,
1096  case_y1x = 6,
1097  case_y2x = 7,
1098  case_y1x2x = 8,
1099  case_xy = 9,
1100 
1101  case_isotropic = case_xy
1102  };
1103  };
1104 
1105 
1106 
1114  template <int dim>
1115  class SubfaceCase : public SubfacePossibilities<dim>
1116  {
1117  public:
1124  subface_possibility);
1125 
1137  operator std::uint8_t() const;
1138 
1142  static constexpr std::size_t
1143  memory_consumption();
1144 
1149  ExcInvalidSubfaceCase,
1150  int,
1151  << "The subface case given (" << arg1 << ") does not make sense "
1152  << "for the space dimension of the object to which they are applied.");
1153 
1154  private:
1159  std::uint8_t value : (dim == 3 ? 4 : 1);
1160  };
1161 
1162 } // namespace internal
1163 
1164 
1165 
1166 template <int dim>
1168 
1169 
1170 
1188 template <>
1189 struct GeometryInfo<0>
1190 {
1198  static constexpr unsigned int max_children_per_cell = 1;
1199 
1203  static constexpr unsigned int faces_per_cell = 0;
1204 
1212  static constexpr unsigned int max_children_per_face = 0;
1213 
1219  static unsigned int
1220  n_children(const RefinementCase<0> &refinement_case);
1221 
1225  static constexpr unsigned int vertices_per_cell = 1;
1226 
1233  static constexpr unsigned int vertices_per_face = 0;
1234 
1238  static constexpr unsigned int lines_per_face = 0;
1239 
1243  static constexpr unsigned int quads_per_face = 0;
1244 
1248  static constexpr unsigned int lines_per_cell = 0;
1249 
1253  static constexpr unsigned int quads_per_cell = 0;
1254 
1258  static constexpr unsigned int hexes_per_cell = 0;
1259 
1277  static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1278 
1292  static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1293 };
1294 
1295 
1296 
1821 template <int dim>
1822 struct GeometryInfo
1823 {
1831  static constexpr unsigned int max_children_per_cell = 1 << dim;
1832 
1836  static constexpr unsigned int faces_per_cell = 2 * dim;
1837 
1845  static constexpr unsigned int max_children_per_face =
1846  GeometryInfo<dim - 1>::max_children_per_cell;
1847 
1851  static constexpr unsigned int vertices_per_cell = 1 << dim;
1852 
1856  static constexpr unsigned int vertices_per_face =
1857  GeometryInfo<dim - 1>::vertices_per_cell;
1858 
1862  static constexpr unsigned int lines_per_face =
1863  GeometryInfo<dim - 1>::lines_per_cell;
1864 
1868  static constexpr unsigned int quads_per_face =
1869  GeometryInfo<dim - 1>::quads_per_cell;
1870 
1880  static constexpr unsigned int lines_per_cell =
1881  (2 * GeometryInfo<dim - 1>::lines_per_cell +
1882  GeometryInfo<dim - 1>::vertices_per_cell);
1883 
1891  static constexpr unsigned int quads_per_cell =
1892  (2 * GeometryInfo<dim - 1>::quads_per_cell +
1893  GeometryInfo<dim - 1>::lines_per_cell);
1894 
1898  static constexpr unsigned int hexes_per_cell =
1899  (2 * GeometryInfo<dim - 1>::hexes_per_cell +
1900  GeometryInfo<dim - 1>::quads_per_cell);
1901 
1919  static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
1920  internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal();
1921 
1935  static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
1936  internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal();
1937 
1948  static constexpr std::array<std::array<unsigned int, dim>, vertices_per_cell>
1949  vertex_to_face =
1950  internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face();
1951 
1956  static unsigned int
1957  n_children(const RefinementCase<dim> &refinement_case);
1958 
1963  static unsigned int
1964  n_subfaces(const internal::SubfaceCase<dim> &subface_case);
1965 
1975  static double
1976  subface_ratio(const internal::SubfaceCase<dim> &subface_case,
1977  const unsigned int subface_no);
1978 
1984  static RefinementCase<dim - 1>
1985  face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
1986  const unsigned int face_no,
1987  const bool face_orientation = true,
1988  const bool face_flip = false,
1989  const bool face_rotation = false);
1990 
1996  static RefinementCase<dim>
1997  min_cell_refinement_case_for_face_refinement(
1998  const RefinementCase<dim - 1> &face_refinement_case,
1999  const unsigned int face_no,
2000  const bool face_orientation = true,
2001  const bool face_flip = false,
2002  const bool face_rotation = false);
2003 
2008  static RefinementCase<1>
2009  line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2010  const unsigned int line_no);
2011 
2016  static RefinementCase<dim>
2017  min_cell_refinement_case_for_line_refinement(const unsigned int line_no);
2018 
2065  static unsigned int
2066  child_cell_on_face(const RefinementCase<dim> & ref_case,
2067  const unsigned int face,
2068  const unsigned int subface,
2069  const bool face_orientation = true,
2070  const bool face_flip = false,
2071  const bool face_rotation = false,
2072  const RefinementCase<dim - 1> &face_refinement_case =
2074 
2088  static unsigned int
2089  line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
2090 
2111  static unsigned int
2112  face_to_cell_vertices(const unsigned int face,
2113  const unsigned int vertex,
2114  const bool face_orientation = true,
2115  const bool face_flip = false,
2116  const bool face_rotation = false);
2117 
2129  static unsigned int
2130  face_to_cell_lines(const unsigned int face,
2131  const unsigned int line,
2132  const bool face_orientation = true,
2133  const bool face_flip = false,
2134  const bool face_rotation = false);
2135 
2145  static unsigned int
2146  standard_to_real_face_vertex(const unsigned int vertex,
2147  const bool face_orientation = true,
2148  const bool face_flip = false,
2149  const bool face_rotation = false);
2150 
2160  static unsigned int
2161  real_to_standard_face_vertex(const unsigned int vertex,
2162  const bool face_orientation = true,
2163  const bool face_flip = false,
2164  const bool face_rotation = false);
2165 
2175  static unsigned int
2176  standard_to_real_face_line(const unsigned int line,
2177  const bool face_orientation = true,
2178  const bool face_flip = false,
2179  const bool face_rotation = false);
2180 
2190  static unsigned int
2191  real_to_standard_face_line(const unsigned int line,
2192  const bool face_orientation = true,
2193  const bool face_flip = false,
2194  const bool face_rotation = false);
2195 
2201  static Point<dim>
2202  unit_cell_vertex(const unsigned int vertex);
2203 
2213  static unsigned int
2214  child_cell_from_point(const Point<dim> &p);
2215 
2223  static Point<dim>
2224  cell_to_child_coordinates(const Point<dim> & p,
2225  const unsigned int child_index,
2226  const RefinementCase<dim> refine_case =
2228 
2234  static Point<dim>
2235  child_to_cell_coordinates(const Point<dim> & p,
2236  const unsigned int child_index,
2237  const RefinementCase<dim> refine_case =
2239 
2244  static bool
2245  is_inside_unit_cell(const Point<dim> &p);
2246 
2258  static bool
2259  is_inside_unit_cell(const Point<dim> &p, const double eps);
2260 
2265  static Point<dim>
2266  project_to_unit_cell(const Point<dim> &p);
2267 
2273  static double
2274  distance_to_unit_cell(const Point<dim> &p);
2275 
2280  static double
2281  d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2282 
2287  static Tensor<1, dim>
2288  d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2289 
2341  template <int spacedim>
2342  static void
2343  alternating_form_at_vertices
2344 #ifndef DEAL_II_CONSTEXPR_BUG
2345  (const Point<spacedim> (&vertices)[vertices_per_cell],
2346  Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2347 #else
2348  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2349 #endif
2350 
2360  static constexpr std::array<unsigned int, faces_per_cell>
2361  unit_normal_direction =
2362  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction();
2363 
2380  static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2381  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation();
2382 
2393  static constexpr std::array<Tensor<1, dim>, faces_per_cell>
2394  unit_normal_vector =
2395  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_vector();
2396 
2402  static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2403  internal::GeometryInfoHelper::Initializers<dim>::opposite_face();
2404 
2405 
2409  DeclException1(ExcInvalidCoordinate,
2410  double,
2411  << "The coordinates must satisfy 0 <= x_i <= 1, "
2412  << "but here we have x_i=" << arg1);
2413 
2417  DeclException3(ExcInvalidSubface,
2418  int,
2419  int,
2420  int,
2421  << "RefinementCase<dim> " << arg1 << ": face " << arg2
2422  << " has no subface " << arg3);
2423 };
2424 
2425 
2426 
2427 #ifndef DOXYGEN
2428 
2429 
2430 /* -------------- declaration of explicit specializations ------------- */
2431 
2432 
2433 template <>
2436  const unsigned int i);
2437 template <>
2440  const unsigned int i);
2441 template <>
2444  const unsigned int i);
2445 
2446 
2447 
2448 /* -------------- inline functions ------------- */
2449 
2450 
2451 inline GeometryPrimitive::GeometryPrimitive(const Object object)
2452  : object(object)
2453 {}
2454 
2455 
2456 
2457 inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2458  : object(static_cast<Object>(object_dimension))
2459 {}
2460 
2461 
2462 inline GeometryPrimitive::operator unsigned int() const
2463 {
2464  return static_cast<unsigned int>(object);
2465 }
2466 
2467 
2468 
2469 namespace internal
2470 {
2471  template <int dim>
2472  inline SubfaceCase<dim>::SubfaceCase(
2473  const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2474  : value(subface_possibility)
2475  {}
2476 
2477 
2478  template <int dim>
2479  inline SubfaceCase<dim>::operator std::uint8_t() const
2480  {
2481  return value;
2482  }
2483 
2484 
2485 } // namespace internal
2486 
2487 
2488 template <int dim>
2489 inline RefinementCase<dim>
2490 RefinementCase<dim>::cut_axis(const unsigned int)
2491 {
2492  Assert(false, ExcInternalError());
2493  return static_cast<std::uint8_t>(-1);
2494 }
2495 
2496 
2497 template <>
2498 inline RefinementCase<1>
2499 RefinementCase<1>::cut_axis(const unsigned int i)
2500 {
2501  Assert(i < 1, ExcIndexRange(i, 0, 1));
2502 
2504  return options[i];
2505 }
2506 
2507 
2508 
2509 template <>
2510 inline RefinementCase<2>
2511 RefinementCase<2>::cut_axis(const unsigned int i)
2512 {
2513  Assert(i < 2, ExcIndexRange(i, 0, 2));
2514 
2517  return options[i];
2518 }
2519 
2520 
2521 
2522 template <>
2523 inline RefinementCase<3>
2524 RefinementCase<3>::cut_axis(const unsigned int i)
2525 {
2526  Assert(i < 3, ExcIndexRange(i, 0, 3));
2527 
2531  return options[i];
2532 }
2533 
2534 
2535 
2536 template <int dim>
2538  : value(RefinementPossibilities<dim>::no_refinement)
2539 {}
2540 
2541 
2542 
2543 template <int dim>
2545  const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2546  : value(refinement_case)
2547 {
2548  // check that only those bits of
2549  // the given argument are set that
2550  // make sense for a given space
2551  // dimension
2552  Assert((refinement_case &
2554  refinement_case,
2555  ExcInvalidRefinementCase(refinement_case));
2556 }
2557 
2558 
2559 
2560 template <int dim>
2561 inline RefinementCase<dim>::RefinementCase(const std::uint8_t refinement_case)
2562  : value(refinement_case)
2563 {
2564  // check that only those bits of
2565  // the given argument are set that
2566  // make sense for a given space
2567  // dimension
2568  Assert((refinement_case &
2570  refinement_case,
2571  ExcInvalidRefinementCase(refinement_case));
2572 }
2573 
2574 
2575 
2576 template <int dim>
2577 inline RefinementCase<dim>::operator std::uint8_t() const
2578 {
2579  return value;
2580 }
2581 
2582 
2583 
2584 template <int dim>
2585 inline RefinementCase<dim>
2587 {
2588  return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2589 }
2590 
2591 
2592 
2593 template <int dim>
2595  operator&(const RefinementCase<dim> &r) const
2596 {
2597  return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2598 }
2599 
2600 
2601 
2602 template <int dim>
2603 inline RefinementCase<dim>
2605 {
2606  return RefinementCase<dim>(static_cast<std::uint8_t>(
2608 }
2609 
2610 
2611 
2612 template <int dim>
2613 inline std::size_t
2615 {
2616  return sizeof(RefinementCase<dim>);
2617 }
2618 
2619 
2620 
2621 template <int dim>
2622 template <class Archive>
2623 inline void
2624 RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2625 {
2626  // serialization can't deal with bitfields, so copy from/to a full sized
2627  // std::uint8_t
2628  std::uint8_t uchar_value = value;
2629  ar & uchar_value;
2630  value = uchar_value;
2631 }
2632 
2633 
2634 
2635 template <>
2636 inline Point<1>
2637 GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2638 {
2639  Assert(vertex < vertices_per_cell,
2640  ExcIndexRange(vertex, 0, vertices_per_cell));
2641 
2642  return Point<1>(static_cast<double>(vertex));
2643 }
2644 
2645 
2646 
2647 template <>
2648 inline Point<2>
2649 GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2650 {
2651  Assert(vertex < vertices_per_cell,
2652  ExcIndexRange(vertex, 0, vertices_per_cell));
2653 
2654  return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2655 }
2656 
2657 
2658 
2659 template <>
2660 inline Point<3>
2661 GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2662 {
2663  Assert(vertex < vertices_per_cell,
2664  ExcIndexRange(vertex, 0, vertices_per_cell));
2665 
2666  return {static_cast<double>(vertex % 2),
2667  static_cast<double>(vertex / 2 % 2),
2668  static_cast<double>(vertex / 4)};
2669 }
2670 
2671 
2672 
2673 template <int dim>
2674 inline Point<dim>
2675 GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2676 {
2677  Assert(false, ExcNotImplemented());
2678 
2679  return Point<dim>();
2680 }
2681 
2682 
2683 
2684 template <>
2685 inline unsigned int
2687 {
2688  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2689 
2690  return (p[0] <= 0.5 ? 0 : 1);
2691 }
2692 
2693 
2694 
2695 template <>
2696 inline unsigned int
2698 {
2699  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2700  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2701 
2702  return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
2703 }
2704 
2705 
2706 
2707 template <>
2708 inline unsigned int
2710 {
2711  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2712  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2713  Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
2714 
2715  return (p[0] <= 0.5 ?
2716  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
2717  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
2718 }
2719 
2720 
2721 template <int dim>
2722 inline unsigned int
2724 {
2725  Assert(false, ExcNotImplemented());
2726 
2727  return 0;
2728 }
2729 
2730 
2731 
2732 template <>
2733 inline Point<1>
2735  const unsigned int child_index,
2736  const RefinementCase<1> refine_case)
2737 
2738 {
2739  Assert(child_index < 2, ExcIndexRange(child_index, 0, 2));
2741  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
2742 
2743  return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
2744 }
2745 
2746 
2747 
2748 template <>
2749 inline Point<2>
2751  const unsigned int child_index,
2752  const RefinementCase<2> refine_case)
2753 
2754 {
2755  Assert(child_index < GeometryInfo<2>::n_children(refine_case),
2756  ExcIndexRange(child_index,
2757  0,
2758  GeometryInfo<2>::n_children(refine_case)));
2759 
2760  Point<2> point = p;
2761  switch (refine_case)
2762  {
2764  point[0] *= 2.0;
2765  if (child_index == 1)
2766  point[0] -= 1.0;
2767  break;
2769  point[1] *= 2.0;
2770  if (child_index == 1)
2771  point[1] -= 1.0;
2772  break;
2774  point *= 2.0;
2775  point -= unit_cell_vertex(child_index);
2776  break;
2777  default:
2778  Assert(false, ExcInternalError());
2779  }
2780 
2781  return point;
2782 }
2783 
2784 
2785 
2786 template <>
2787 inline Point<3>
2789  const unsigned int child_index,
2790  const RefinementCase<3> refine_case)
2791 
2792 {
2793  Assert(child_index < GeometryInfo<3>::n_children(refine_case),
2794  ExcIndexRange(child_index,
2795  0,
2796  GeometryInfo<3>::n_children(refine_case)));
2797 
2798  Point<3> point = p;
2799  // there might be a cleverer way to do
2800  // this, but since this function is called
2801  // in very few places for initialization
2802  // purposes only, I don't care at the
2803  // moment
2804  switch (refine_case)
2805  {
2807  point[0] *= 2.0;
2808  if (child_index == 1)
2809  point[0] -= 1.0;
2810  break;
2812  point[1] *= 2.0;
2813  if (child_index == 1)
2814  point[1] -= 1.0;
2815  break;
2817  point[2] *= 2.0;
2818  if (child_index == 1)
2819  point[2] -= 1.0;
2820  break;
2822  point[0] *= 2.0;
2823  point[1] *= 2.0;
2824  if (child_index % 2 == 1)
2825  point[0] -= 1.0;
2826  if (child_index / 2 == 1)
2827  point[1] -= 1.0;
2828  break;
2830  // careful, this is slightly
2831  // different from xy and yz due to
2832  // different internal numbering of
2833  // children!
2834  point[0] *= 2.0;
2835  point[2] *= 2.0;
2836  if (child_index / 2 == 1)
2837  point[0] -= 1.0;
2838  if (child_index % 2 == 1)
2839  point[2] -= 1.0;
2840  break;
2842  point[1] *= 2.0;
2843  point[2] *= 2.0;
2844  if (child_index % 2 == 1)
2845  point[1] -= 1.0;
2846  if (child_index / 2 == 1)
2847  point[2] -= 1.0;
2848  break;
2850  point *= 2.0;
2851  point -= unit_cell_vertex(child_index);
2852  break;
2853  default:
2854  Assert(false, ExcInternalError());
2855  }
2856 
2857  return point;
2858 }
2859 
2860 
2861 
2862 template <int dim>
2863 inline Point<dim>
2865  const Point<dim> & /*p*/,
2866  const unsigned int /*child_index*/,
2867  const RefinementCase<dim> /*refine_case*/)
2868 
2869 {
2870  Assert(false, ExcNotImplemented());
2871  return Point<dim>();
2872 }
2873 
2874 
2875 
2876 template <>
2877 inline Point<1>
2879  const unsigned int child_index,
2880  const RefinementCase<1> refine_case)
2881 
2882 {
2883  Assert(child_index < 2, ExcIndexRange(child_index, 0, 2));
2885  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
2886 
2887  return (p + unit_cell_vertex(child_index)) * 0.5;
2888 }
2889 
2890 
2891 
2892 template <>
2893 inline Point<3>
2895  const unsigned int child_index,
2896  const RefinementCase<3> refine_case)
2897 
2898 {
2899  Assert(child_index < GeometryInfo<3>::n_children(refine_case),
2900  ExcIndexRange(child_index,
2901  0,
2902  GeometryInfo<3>::n_children(refine_case)));
2903 
2904  Point<3> point = p;
2905  // there might be a cleverer way to do
2906  // this, but since this function is called
2907  // in very few places for initialization
2908  // purposes only, I don't care at the
2909  // moment
2910  switch (refine_case)
2911  {
2913  if (child_index == 1)
2914  point[0] += 1.0;
2915  point[0] *= 0.5;
2916  break;
2918  if (child_index == 1)
2919  point[1] += 1.0;
2920  point[1] *= 0.5;
2921  break;
2923  if (child_index == 1)
2924  point[2] += 1.0;
2925  point[2] *= 0.5;
2926  break;
2928  if (child_index % 2 == 1)
2929  point[0] += 1.0;
2930  if (child_index / 2 == 1)
2931  point[1] += 1.0;
2932  point[0] *= 0.5;
2933  point[1] *= 0.5;
2934  break;
2936  // careful, this is slightly
2937  // different from xy and yz due to
2938  // different internal numbering of
2939  // children!
2940  if (child_index / 2 == 1)
2941  point[0] += 1.0;
2942  if (child_index % 2 == 1)
2943  point[2] += 1.0;
2944  point[0] *= 0.5;
2945  point[2] *= 0.5;
2946  break;
2948  if (child_index % 2 == 1)
2949  point[1] += 1.0;
2950  if (child_index / 2 == 1)
2951  point[2] += 1.0;
2952  point[1] *= 0.5;
2953  point[2] *= 0.5;
2954  break;
2956  point += unit_cell_vertex(child_index);
2957  point *= 0.5;
2958  break;
2959  default:
2960  Assert(false, ExcInternalError());
2961  }
2962 
2963  return point;
2964 }
2965 
2966 
2967 
2968 template <>
2969 inline Point<2>
2971  const unsigned int child_index,
2972  const RefinementCase<2> refine_case)
2973 {
2974  Assert(child_index < GeometryInfo<2>::n_children(refine_case),
2975  ExcIndexRange(child_index,
2976  0,
2977  GeometryInfo<2>::n_children(refine_case)));
2978 
2979  Point<2> point = p;
2980  switch (refine_case)
2981  {
2983  if (child_index == 1)
2984  point[0] += 1.0;
2985  point[0] *= 0.5;
2986  break;
2988  if (child_index == 1)
2989  point[1] += 1.0;
2990  point[1] *= 0.5;
2991  break;
2993  point += unit_cell_vertex(child_index);
2994  point *= 0.5;
2995  break;
2996  default:
2997  Assert(false, ExcInternalError());
2998  }
2999 
3000  return point;
3001 }
3002 
3003 
3004 
3005 template <int dim>
3006 inline Point<dim>
3008  const Point<dim> & /*p*/,
3009  const unsigned int /*child_index*/,
3010  const RefinementCase<dim> /*refine_case*/)
3011 {
3012  Assert(false, ExcNotImplemented());
3013  return Point<dim>();
3014 }
3015 
3016 
3017 
3018 template <int dim>
3019 inline bool
3021 {
3022  Assert(false, ExcNotImplemented());
3023  return false;
3024 }
3025 
3026 template <>
3027 inline bool
3029 {
3030  return (p[0] >= 0.) && (p[0] <= 1.);
3031 }
3032 
3033 
3034 
3035 template <>
3036 inline bool
3038 {
3039  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
3040 }
3041 
3042 
3043 
3044 template <>
3045 inline bool
3047 {
3048  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
3049  (p[2] >= 0.) && (p[2] <= 1.);
3050 }
3051 
3052 
3053 
3054 template <int dim>
3055 inline bool
3057 {
3058  Assert(false, ExcNotImplemented());
3059  return false;
3060 }
3061 
3062 template <>
3063 inline bool
3064 GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
3065 {
3066  return (p[0] >= -eps) && (p[0] <= 1. + eps);
3067 }
3068 
3069 
3070 
3071 template <>
3072 inline bool
3073 GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
3074 {
3075  const double l = -eps, u = 1 + eps;
3076  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
3077 }
3078 
3079 
3080 
3081 template <>
3082 inline bool
3083 GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
3084 {
3085  const double l = -eps, u = 1.0 + eps;
3086  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
3087  (p[2] >= l) && (p[2] <= u);
3088 }
3089 
3090 
3091 
3092 template <>
3093 inline unsigned int
3094 GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
3095  const unsigned int vertex)
3096 {
3097  (void)line;
3099  Assert(vertex < 2, ExcIndexRange(vertex, 0, 2));
3100 
3101  return vertex;
3102 }
3103 
3104 
3105 template <>
3106 inline unsigned int
3107 GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
3108  const unsigned int vertex)
3109 {
3110  constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
3111  return cell_vertices[line][vertex];
3112 }
3113 
3114 
3115 
3116 template <>
3117 inline unsigned int
3118 GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
3119  const unsigned int vertex)
3120 {
3122  Assert(vertex < 2, ExcIndexRange(vertex, 0, 2));
3123 
3124  constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3125  {1, 3},
3126  {0, 1},
3127  {2, 3},
3128  {4, 6}, // top face
3129  {5, 7},
3130  {4, 5},
3131  {6, 7},
3132  {0,
3133  4}, // connects of bottom
3134  {1, 5}, // top face
3135  {2, 6},
3136  {3, 7}};
3137  return vertices[line][vertex];
3138 }
3139 
3140 
3141 template <>
3142 inline unsigned int
3143 GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3144 {
3145  Assert(false, ExcNotImplemented());
3147 }
3148 
3149 template <>
3150 inline unsigned int
3151 GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3152  const bool face_orientation,
3153  const bool face_flip,
3154  const bool face_rotation)
3155 {
3158 
3159  // set up a table to make sure that
3160  // we handle non-standard faces correctly
3161  //
3162  // so set up a table that for each vertex (of
3163  // a quad in standard position) describes
3164  // which vertex to take
3165  //
3166  // first index: four vertices 0...3
3167  //
3168  // second index: face_orientation; 0:
3169  // opposite normal, 1: standard
3170  //
3171  // third index: face_flip; 0: standard, 1:
3172  // face rotated by 180 degrees
3173  //
3174  // forth index: face_rotation: 0: standard,
3175  // 1: face rotated by 90 degrees
3176 
3177  constexpr unsigned int vertex_translation[4][2][2][2] = {
3178  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3179  // face_rotation=false and true
3180  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3181  // face_rotation=false and true
3182  {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3183  // face_rotation=false and true
3184  {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3185  // face_rotation=false and true
3186 
3187  {{{2, 3}, // vertex 1 ...
3188  {1, 0}},
3189  {{1, 0}, {2, 3}}},
3190 
3191  {{{1, 0}, // vertex 2 ...
3192  {2, 3}},
3193  {{2, 3}, {1, 0}}},
3194 
3195  {{{3, 1}, // vertex 3 ...
3196  {0, 2}},
3197  {{3, 1}, {0, 2}}}};
3198 
3199  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3200 }
3201 
3202 
3203 
3204 template <int dim>
3205 inline unsigned int
3206 GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3207  const bool,
3208  const bool,
3209  const bool)
3210 {
3211  Assert(dim > 1, ExcImpossibleInDim(dim));
3214  return vertex;
3215 }
3216 
3217 template <int dim>
3218 inline unsigned int
3220 {
3221  constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3222  0, 2, 2, 4, 2, 4, 4, 8};
3223 
3224  return n_children[ref_case];
3225 }
3226 
3227 
3228 
3229 template <int dim>
3230 inline unsigned int
3232 {
3233  Assert(false, ExcNotImplemented());
3234  return 0;
3235 }
3236 
3237 template <>
3238 inline unsigned int
3240 {
3241  Assert(false, ExcImpossibleInDim(1));
3242  return 0;
3243 }
3244 
3245 template <>
3246 inline unsigned int
3248 {
3249  return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3250 }
3251 
3252 
3253 
3254 template <>
3255 inline unsigned int
3257 {
3258  const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3259  0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3260  return nsubs[subface_case];
3261 }
3262 
3263 
3264 
3265 template <int dim>
3266 inline double
3268  const unsigned int)
3269 {
3270  Assert(false, ExcNotImplemented());
3271  return 0.;
3272 }
3273 
3274 template <>
3275 inline double
3277  const unsigned int)
3278 {
3279  return 1;
3280 }
3281 
3282 
3283 template <>
3284 inline double
3286  const unsigned int)
3287 {
3288  double ratio = 1;
3289  switch (subface_case)
3290  {
3292  // Here, an
3293  // Assert(false,ExcInternalError())
3294  // would be the right
3295  // choice, but
3296  // unfortunately the
3297  // current function is
3298  // also called for faces
3299  // without children (see
3300  // tests/fe/mapping.cc).
3301  // Assert(false, ExcMessage("Face has no subfaces."));
3302  // Furthermore, assign
3303  // following value as
3304  // otherwise the
3305  // bits/volume_x tests
3306  // break
3308  break;
3310  ratio = 0.5;
3311  break;
3312  default:
3313  // there should be no
3314  // cases left
3315  Assert(false, ExcInternalError());
3316  break;
3317  }
3318 
3319  return ratio;
3320 }
3321 
3322 
3323 template <>
3324 inline double
3326  const unsigned int subface_no)
3327 {
3328  double ratio = 1;
3329  switch (subface_case)
3330  {
3332  // Here, an
3333  // Assert(false,ExcInternalError())
3334  // would be the right
3335  // choice, but
3336  // unfortunately the
3337  // current function is
3338  // also called for faces
3339  // without children (see
3340  // tests/bits/mesh_3d_16.cc). Add
3341  // following switch to
3342  // avoid diffs in
3343  // tests/bits/mesh_3d_16
3345  break;
3348  ratio = 0.5;
3349  break;
3353  ratio = 0.25;
3354  break;
3357  if (subface_no < 2)
3358  ratio = 0.25;
3359  else
3360  ratio = 0.5;
3361  break;
3364  if (subface_no == 0)
3365  ratio = 0.5;
3366  else
3367  ratio = 0.25;
3368  break;
3369  default:
3370  // there should be no
3371  // cases left
3372  Assert(false, ExcInternalError());
3373  break;
3374  }
3375 
3376  return ratio;
3377 }
3378 
3379 
3380 
3381 template <int dim>
3383  const RefinementCase<dim> &,
3384  const unsigned int,
3385  const bool,
3386  const bool,
3387  const bool)
3388 {
3389  Assert(false, ExcNotImplemented());
3390  return RefinementCase<dim - 1>::no_refinement;
3391 }
3392 
3393 template <>
3395  const RefinementCase<1> &,
3396  const unsigned int,
3397  const bool,
3398  const bool,
3399  const bool)
3400 {
3401  Assert(false, ExcImpossibleInDim(1));
3402 
3404 }
3405 
3406 
3407 template <>
3408 inline RefinementCase<1>
3410  const RefinementCase<2> &cell_refinement_case,
3411  const unsigned int face_no,
3412  const bool,
3413  const bool,
3414  const bool)
3415 {
3416  const unsigned int dim = 2;
3417  Assert(cell_refinement_case < RefinementCase<dim>::isotropic_refinement + 1,
3418  ExcIndexRange(cell_refinement_case,
3419  0,
3423 
3424  const RefinementCase<dim - 1>
3427  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3428  RefinementCase<dim - 1>::no_refinement},
3429 
3430  {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3431 
3432  {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3433 
3434  {RefinementCase<dim - 1>::cut_x, // cut_xy
3435  RefinementCase<dim - 1>::cut_x}};
3436 
3437  return ref_cases[cell_refinement_case][face_no / 2];
3438 }
3439 
3440 
3441 template <>
3442 inline RefinementCase<2>
3444  const RefinementCase<3> &cell_refinement_case,
3445  const unsigned int face_no,
3446  const bool face_orientation,
3447  const bool /*face_flip*/,
3448  const bool face_rotation)
3449 {
3450  const unsigned int dim = 3;
3451  Assert(cell_refinement_case < RefinementCase<dim>::isotropic_refinement + 1,
3452  ExcIndexRange(cell_refinement_case,
3453  0,
3457 
3458  const RefinementCase<dim - 1>
3461  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3462  RefinementCase<dim - 1>::no_refinement,
3463  RefinementCase<dim - 1>::no_refinement},
3464 
3465  {RefinementCase<dim - 1>::no_refinement, // cut_x
3466  RefinementCase<dim - 1>::cut_y,
3467  RefinementCase<dim - 1>::cut_x},
3468 
3469  {RefinementCase<dim - 1>::cut_x, // cut_y
3470  RefinementCase<dim - 1>::no_refinement,
3471  RefinementCase<dim - 1>::cut_y},
3472 
3473  {RefinementCase<dim - 1>::cut_x, // cut_xy
3474  RefinementCase<dim - 1>::cut_y,
3475  RefinementCase<dim - 1>::cut_xy},
3476 
3477  {RefinementCase<dim - 1>::cut_y, // cut_z
3478  RefinementCase<dim - 1>::cut_x,
3479  RefinementCase<dim - 1>::no_refinement},
3480 
3481  {RefinementCase<dim - 1>::cut_y, // cut_xz
3482  RefinementCase<dim - 1>::cut_xy,
3483  RefinementCase<dim - 1>::cut_x},
3484 
3485  {RefinementCase<dim - 1>::cut_xy, // cut_yz
3486  RefinementCase<dim - 1>::cut_x,
3487  RefinementCase<dim - 1>::cut_y},
3488 
3489  {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3490  RefinementCase<dim - 1>::cut_xy,
3491  RefinementCase<dim - 1>::cut_xy},
3492  };
3493 
3494  const RefinementCase<dim - 1> ref_case =
3495  ref_cases[cell_refinement_case][face_no / 2];
3496 
3497  const RefinementCase<dim - 1> flip[4] = {
3498  RefinementCase<dim - 1>::no_refinement,
3499  RefinementCase<dim - 1>::cut_y,
3500  RefinementCase<dim - 1>::cut_x,
3501  RefinementCase<dim - 1>::cut_xy};
3502 
3503  // correct the ref_case for face_orientation
3504  // and face_rotation. for face_orientation,
3505  // 'true' is the default value whereas for
3506  // face_rotation, 'false' is standard. If
3507  // <tt>face_rotation==face_orientation</tt>,
3508  // then one of them is non-standard and we
3509  // have to swap cut_x and cut_y, otherwise no
3510  // change is necessary. face_flip has no
3511  // influence. however, in order to keep the
3512  // interface consistent with other functions,
3513  // we still include it as an argument to this
3514  // function
3515  return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3516 }
3517 
3518 
3519 
3520 template <int dim>
3521 inline RefinementCase<1>
3523  const unsigned int)
3524 {
3525  Assert(false, ExcNotImplemented());
3527 }
3528 
3529 template <>
3530 inline RefinementCase<1>
3532  const RefinementCase<1> &cell_refinement_case,
3533  const unsigned int line_no)
3534 {
3535  (void)line_no;
3536  const unsigned int dim = 1;
3537  (void)dim;
3538  Assert(cell_refinement_case < RefinementCase<dim>::isotropic_refinement + 1,
3539  ExcIndexRange(cell_refinement_case,
3540  0,
3544 
3545  return cell_refinement_case;
3546 }
3547 
3548 
3549 template <>
3550 inline RefinementCase<1>
3552  const RefinementCase<2> &cell_refinement_case,
3553  const unsigned int line_no)
3554 {
3555  // Assertions are in face_refinement_case()
3556  return face_refinement_case(cell_refinement_case, line_no);
3557 }
3558 
3559 
3560 template <>
3561 inline RefinementCase<1>
3563  const RefinementCase<3> &cell_refinement_case,
3564  const unsigned int line_no)
3565 {
3566  const unsigned int dim = 3;
3567  Assert(cell_refinement_case < RefinementCase<dim>::isotropic_refinement + 1,
3568  ExcIndexRange(cell_refinement_case,
3569  0,
3573 
3574  // array indicating, which simple refine
3575  // case cuts a line in direction x, y or
3576  // z. For example, cut_y and everything
3577  // containing cut_y (cut_xy, cut_yz,
3578  // cut_xyz) cuts lines, which are in y
3579  // direction.
3580  const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3583 
3584  // order the direction of lines
3585  // 0->x, 1->y, 2->z
3586  const unsigned int direction[lines_per_cell] = {
3587  1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3588 
3589  return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3592 }
3593 
3594 
3595 
3596 template <int dim>
3597 inline RefinementCase<dim>
3599  const RefinementCase<dim - 1> &,
3600  const unsigned int,
3601  const bool,
3602  const bool,
3603  const bool)
3604 {
3605  Assert(false, ExcNotImplemented());
3606 
3608 }
3609 
3610 template <>
3611 inline RefinementCase<1>
3613  const RefinementCase<0> &,
3614  const unsigned int,
3615  const bool,
3616  const bool,
3617  const bool)
3618 {
3619  const unsigned int dim = 1;
3620  Assert(false, ExcImpossibleInDim(dim));
3621 
3623 }
3624 
3625 
3626 template <>
3627 inline RefinementCase<2>
3630  const unsigned int face_no,
3631  const bool,
3632  const bool,
3633  const bool)
3634 {
3635  const unsigned int dim = 2;
3636  Assert(face_refinement_case <
3638  ExcIndexRange(face_refinement_case,
3639  0,
3643 
3644  if (face_refinement_case == RefinementCase<dim>::cut_x)
3645  return (face_no / 2) ? RefinementCase<dim>::cut_x :
3647  else
3649 }
3650 
3651 
3652 template <>
3653 inline RefinementCase<3>
3655  const RefinementCase<2> &face_refinement_case,
3656  const unsigned int face_no,
3657  const bool face_orientation,
3658  const bool /*face_flip*/,
3659  const bool face_rotation)
3660 {
3661  const unsigned int dim = 3;
3662  Assert(face_refinement_case <
3664  ExcIndexRange(face_refinement_case,
3665  0,
3669 
3674 
3675  // correct the face_refinement_case for
3676  // face_orientation and face_rotation. for
3677  // face_orientation, 'true' is the default
3678  // value whereas for face_rotation, 'false'
3679  // is standard. If
3680  // <tt>face_rotation==face_orientation</tt>,
3681  // then one of them is non-standard and we
3682  // have to swap cut_x and cut_y, otherwise no
3683  // change is necessary. face_flip has no
3684  // influence. however, in order to keep the
3685  // interface consistent with other functions,
3686  // we still include it as an argument to this
3687  // function
3688  const RefinementCase<dim - 1> std_face_ref =
3689  (face_orientation == face_rotation) ? flip[face_refinement_case] :
3690  face_refinement_case;
3691 
3692  const RefinementCase<dim> face_to_cell[3][4] = {
3693  {RefinementCase<dim>::no_refinement, // faces 0 and 1
3694  RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3695  RefinementCase<dim>::cut_z,
3697 
3698  {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3699  // "exchanged on faces 2 and 3")
3700  RefinementCase<dim>::cut_z,
3703 
3704  {RefinementCase<dim>::no_refinement, // faces 4 and 5
3708 
3709  return face_to_cell[face_no / 2][std_face_ref];
3710 }
3711 
3712 
3713 
3714 template <int dim>
3715 inline RefinementCase<dim>
3717  const unsigned int)
3718 {
3719  Assert(false, ExcNotImplemented());
3720 
3722 }
3723 
3724 template <>
3725 inline RefinementCase<1>
3727  const unsigned int line_no)
3728 {
3729  (void)line_no;
3730  Assert(line_no == 0, ExcIndexRange(line_no, 0, 1));
3731 
3732  return RefinementCase<1>::cut_x;
3733 }
3734 
3735 
3736 template <>
3737 inline RefinementCase<2>
3739  const unsigned int line_no)
3740 {
3741  const unsigned int dim = 2;
3742  (void)dim;
3745 
3746  return (line_no / 2) ? RefinementCase<2>::cut_x : RefinementCase<2>::cut_y;
3747 }
3748 
3749 
3750 template <>
3751 inline RefinementCase<3>
3753  const unsigned int line_no)
3754 {
3755  const unsigned int dim = 3;
3758 
3759  const RefinementCase<dim> ref_cases[6] = {
3760  RefinementCase<dim>::cut_y, // lines 0 and 1
3761  RefinementCase<dim>::cut_x, // lines 2 and 3
3762  RefinementCase<dim>::cut_y, // lines 4 and 5
3763  RefinementCase<dim>::cut_x, // lines 6 and 7
3764  RefinementCase<dim>::cut_z, // lines 8 and 9
3765  RefinementCase<dim>::cut_z}; // lines 10 and 11
3766 
3767  return ref_cases[line_no / 2];
3768 }
3769 
3770 
3771 
3772 template <>
3773 inline unsigned int
3774 GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
3775  const bool face_orientation,
3776  const bool face_flip,
3777  const bool face_rotation)
3778 {
3781 
3782  // set up a table to make sure that
3783  // we handle non-standard faces correctly
3784  //
3785  // so set up a table that for each vertex (of
3786  // a quad in standard position) describes
3787  // which vertex to take
3788  //
3789  // first index: four vertices 0...3
3790  //
3791  // second index: face_orientation; 0:
3792  // opposite normal, 1: standard
3793  //
3794  // third index: face_flip; 0: standard, 1:
3795  // face rotated by 180 degrees
3796  //
3797  // forth index: face_rotation: 0: standard,
3798  // 1: face rotated by 90 degrees
3799 
3800  const unsigned int vertex_translation[4][2][2][2] = {
3801  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3802  // face_rotation=false and true
3803  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3804  // face_rotation=false and true
3805  {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
3806  // face_rotation=false and true
3807  {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
3808  // face_rotation=false and true
3809 
3810  {{{2, 3}, // vertex 1 ...
3811  {1, 0}},
3812  {{1, 3}, {2, 0}}},
3813 
3814  {{{1, 0}, // vertex 2 ...
3815  {2, 3}},
3816  {{2, 0}, {1, 3}}},
3817 
3818  {{{3, 1}, // vertex 3 ...
3819  {0, 2}},
3820  {{3, 2}, {0, 1}}}};
3821 
3822  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3823 }
3824 
3825 
3826 
3827 template <int dim>
3828 inline unsigned int
3829 GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
3830  const bool,
3831  const bool,
3832  const bool)
3833 {
3834  Assert(dim > 1, ExcImpossibleInDim(dim));
3837  return vertex;
3838 }
3839 
3840 
3841 
3842 template <>
3843 inline unsigned int
3844 GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
3845  const bool face_orientation,
3846  const bool face_flip,
3847  const bool face_rotation)
3848 {
3851 
3852 
3853  // make sure we handle
3854  // non-standard faces correctly
3855  //
3856  // so set up a table that for each line (of a
3857  // quad) describes which line to take
3858  //
3859  // first index: four lines 0...3
3860  //
3861  // second index: face_orientation; 0:
3862  // opposite normal, 1: standard
3863  //
3864  // third index: face_flip; 0: standard, 1:
3865  // face rotated by 180 degrees
3866  //
3867  // forth index: face_rotation: 0: standard,
3868  // 1: face rotated by 90 degrees
3869 
3870  const unsigned int line_translation[4][2][2][2] = {
3871  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
3872  // face_rotation=false and true
3873  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
3874  // face_rotation=false and true
3875  {{0, 3}, // line 0, face_orientation=true, face_flip=false,
3876  // face_rotation=false and true
3877  {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
3878  // face_rotation=false and true
3879 
3880  {{{3, 1}, // line 1 ...
3881  {2, 0}},
3882  {{1, 2}, {0, 3}}},
3883 
3884  {{{0, 3}, // line 2 ...
3885  {1, 2}},
3886  {{2, 0}, {3, 1}}},
3887 
3888  {{{1, 2}, // line 3 ...
3889  {0, 3}},
3890  {{3, 1}, {2, 0}}}};
3891 
3892  return line_translation[line][face_orientation][face_flip][face_rotation];
3893 }
3894 
3895 
3896 
3897 template <int dim>
3898 inline unsigned int
3899 GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
3900  const bool,
3901  const bool,
3902  const bool)
3903 {
3904  Assert(false, ExcNotImplemented());
3905  return line;
3906 }
3907 
3908 
3909 
3910 template <>
3911 inline unsigned int
3912 GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
3913  const bool face_orientation,
3914  const bool face_flip,
3915  const bool face_rotation)
3916 {
3919 
3920 
3921  // make sure we handle
3922  // non-standard faces correctly
3923  //
3924  // so set up a table that for each line (of a
3925  // quad) describes which line to take
3926  //
3927  // first index: four lines 0...3
3928  //
3929  // second index: face_orientation; 0:
3930  // opposite normal, 1: standard
3931  //
3932  // third index: face_flip; 0: standard, 1:
3933  // face rotated by 180 degrees
3934  //
3935  // forth index: face_rotation: 0: standard,
3936  // 1: face rotated by 90 degrees
3937 
3938  const unsigned int line_translation[4][2][2][2] = {
3939  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
3940  // face_rotation=false and true
3941  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
3942  // face_rotation=false and true
3943  {{0, 2}, // line 0, face_orientation=true, face_flip=false,
3944  // face_rotation=false and true
3945  {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
3946  // face_rotation=false and true
3947 
3948  {{{3, 1}, // line 1 ...
3949  {2, 0}},
3950  {{1, 3}, {0, 2}}},
3951 
3952  {{{0, 3}, // line 2 ...
3953  {1, 2}},
3954  {{2, 1}, {3, 0}}},
3955 
3956  {{{1, 2}, // line 3 ...
3957  {0, 3}},
3958  {{3, 0}, {2, 1}}}};
3959 
3960  return line_translation[line][face_orientation][face_flip][face_rotation];
3961 }
3962 
3963 
3964 
3965 template <int dim>
3966 inline unsigned int
3967 GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
3968  const bool,
3969  const bool,
3970  const bool)
3971 {
3972  Assert(false, ExcNotImplemented());
3973  return line;
3974 }
3975 
3976 
3977 
3978 template <>
3979 inline unsigned int
3981  const unsigned int face,
3982  const unsigned int subface,
3983  const bool,
3984  const bool,
3985  const bool,
3986  const RefinementCase<0> &)
3987 {
3988  (void)subface;
3990  Assert(subface < max_children_per_face,
3991  ExcIndexRange(subface, 0, max_children_per_face));
3992 
3993  return face;
3994 }
3995 
3996 
3997 
3998 template <>
3999 inline unsigned int
4001  const unsigned int face,
4002  const unsigned int subface,
4003  const bool /*face_orientation*/,
4004  const bool face_flip,
4005  const bool /*face_rotation*/,
4006  const RefinementCase<1> &)
4007 {
4009  Assert(subface < max_children_per_face,
4010  ExcIndexRange(subface, 0, max_children_per_face));
4011 
4012  // always return the child adjacent to the specified
4013  // subface. if the face of a cell is not refined, don't
4014  // throw an assertion but deliver the child adjacent to
4015  // the face nevertheless, i.e. deliver the child of
4016  // this cell adjacent to the subface of a possibly
4017  // refined neighbor. this simplifies setting neighbor
4018  // information in execute_refinement.
4019  const unsigned int
4021  [max_children_per_face] = {
4022  {
4023  // Normal orientation (face_flip = false)
4024  {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
4025  {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
4026  {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
4027  },
4028  {
4029  // Flipped orientation (face_flip = true)
4030  {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
4031  {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
4032  {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
4033  }};
4034 
4035  return subcells[face_flip][ref_case - 1][face][subface];
4036 }
4037 
4038 
4039 
4040 template <>
4041 inline unsigned int
4043  const unsigned int face,
4044  const unsigned int subface,
4045  const bool face_orientation,
4046  const bool face_flip,
4047  const bool face_rotation,
4048  const RefinementCase<2> &face_ref_case)
4049 {
4050  const unsigned int dim = 3;
4051 
4053  ExcMessage("Cell has no children."));
4055  Assert(subface < GeometryInfo<dim - 1>::n_children(face_ref_case) ||
4056  (subface == 0 &&
4057  face_ref_case == RefinementCase<dim - 1>::no_refinement),
4058  ExcIndexRange(subface, 0, GeometryInfo<2>::n_children(face_ref_case)));
4059 
4060  // invalid number used for invalid cases,
4061  // e.g. when the children are more refined at
4062  // a given face than the face itself
4063  const unsigned int e = numbers::invalid_unsigned_int;
4064 
4065  // the whole process of finding a child cell
4066  // at a given subface considering the
4067  // possibly anisotropic refinement cases of
4068  // the cell and the face as well as
4069  // orientation, flip and rotation of the face
4070  // is quite complicated. thus, we break it
4071  // down into several steps.
4072 
4073  // first step: convert the given face refine
4074  // case to a face refine case concerning the
4075  // face in standard orientation (, flip and
4076  // rotation). This only affects cut_x and
4077  // cut_y
4078  const RefinementCase<dim - 1> flip[4] = {
4079  RefinementCase<dim - 1>::no_refinement,
4080  RefinementCase<dim - 1>::cut_y,
4081  RefinementCase<dim - 1>::cut_x,
4082  RefinementCase<dim - 1>::cut_xy};
4083  // for face_orientation, 'true' is the
4084  // default value whereas for face_rotation,
4085  // 'false' is standard. If
4086  // <tt>face_rotation==face_orientation</tt>,
4087  // then one of them is non-standard and we
4088  // have to swap cut_x and cut_y, otherwise no
4089  // change is necessary.
4090  const RefinementCase<dim - 1> std_face_ref =
4091  (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4092 
4093  // second step: convert the given subface
4094  // index to the one for a standard face
4095  // respecting face_orientation, face_flip and
4096  // face_rotation
4097 
4098  // first index: face_ref_case
4099  // second index: face_orientation
4100  // third index: face_flip
4101  // forth index: face_rotation
4102  // fifth index: subface index
4103  const unsigned int subface_exchange[4][2][2][2][4] = {
4104  // no_refinement (subface 0 stays 0,
4105  // all others are invalid)
4106  {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4107  {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4108  // cut_x (here, if the face is only
4109  // rotated OR only falsely oriented,
4110  // then subface 0 of the non-standard
4111  // face does NOT correspond to one of
4112  // the subfaces of a standard
4113  // face. Thus we indicate the subface
4114  // which is located at the lower left
4115  // corner (the origin of the face's
4116  // local coordinate system) with
4117  // '0'. The rest of this issue is
4118  // taken care of using the above
4119  // conversion to a 'standard face
4120  // refine case')
4121  {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4122  {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4123  // cut_y (the same applies as for
4124  // cut_x)
4125  {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4126  {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4127  // cut_xyz: this information is
4128  // identical to the information
4129  // returned by
4130  // GeometryInfo<3>::real_to_standard_face_vertex()
4131  {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4132  // face_rotation=false, subfaces 0,1,2,3
4133  {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4134  // face_rotation=true, subfaces 0,1,2,3
4135  {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4136  // face_rotation=false, subfaces 0,1,2,3
4137  {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4138  // face_rotation=true, subfaces 0,1,2,3
4139  {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4140  // face_rotation=false, subfaces 0,1,2,3
4141  {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4142  // face_rotation=true, subfaces 0,1,2,3
4143  {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4144  // face_rotation=false, subfaces 0,1,2,3
4145  {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4146  // face_rotation=true, subfaces 0,1,2,3
4147 
4148  const unsigned int std_subface =
4149  subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4150  [subface];
4151  Assert(std_subface != e, ExcInternalError());
4152 
4153  // third step: these are the children, which
4154  // can be found at the given subfaces of an
4155  // isotropically refined (standard) face
4156  //
4157  // first index: (refinement_case-1)
4158  // second index: face_index
4159  // third index: subface_index (isotropic refinement)
4160  const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4161  [max_children_per_face] = {
4162  // cut_x
4163  {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4164  {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4165  {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4166  {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4167  {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4168  {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4169  // cut_y
4170  {{0, 1, 0, 1},
4171  {0, 1, 0, 1},
4172  {0, 0, 0, 0},
4173  {1, 1, 1, 1},
4174  {0, 0, 1, 1},
4175  {0, 0, 1, 1}},
4176  // cut_xy
4177  {{0, 2, 0, 2},
4178  {1, 3, 1, 3},
4179  {0, 0, 1, 1},
4180  {2, 2, 3, 3},
4181  {0, 1, 2, 3},
4182  {0, 1, 2, 3}},
4183  // cut_z
4184  {{0, 0, 1, 1},
4185  {0, 0, 1, 1},
4186  {0, 1, 0, 1},
4187  {0, 1, 0, 1},
4188  {0, 0, 0, 0},
4189  {1, 1, 1, 1}},
4190  // cut_xz
4191  {{0, 0, 1, 1},
4192  {2, 2, 3, 3},
4193  {0, 1, 2, 3},
4194  {0, 1, 2, 3},
4195  {0, 2, 0, 2},
4196  {1, 3, 1, 3}},
4197  // cut_yz
4198  {{0, 1, 2, 3},
4199  {0, 1, 2, 3},
4200  {0, 2, 0, 2},
4201  {1, 3, 1, 3},
4202  {0, 0, 1, 1},
4203  {2, 2, 3, 3}},
4204  // cut_xyz
4205  {{0, 2, 4, 6},
4206  {1, 3, 5, 7},
4207  {0, 4, 1, 5},
4208  {2, 6, 3, 7},
4209  {0, 1, 2, 3},
4210  {4, 5, 6, 7}}};
4211 
4212  // forth step: check, whether the given face
4213  // refine case is valid for the given cell
4214  // refine case. this is the case, if the
4215  // given face refine case is at least as
4216  // refined as the face is for the given cell
4217  // refine case
4218 
4219  // note, that we are considering standard
4220  // face refinement cases here and thus must
4221  // not pass the given orientation, flip and
4222  // rotation flags
4223  if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4224  face_refinement_case(ref_case, face))
4225  {
4226  // all is fine. for anisotropic face
4227  // refine cases, select one of the
4228  // isotropic subfaces which neighbors the
4229  // same child
4230 
4231  // first index: (standard) face refine case
4232  // second index: subface index
4233  const unsigned int equivalent_iso_subface[4][4] = {
4234  {0, e, e, e}, // no_refinement
4235  {0, 3, e, e}, // cut_x
4236  {0, 3, e, e}, // cut_y
4237  {0, 1, 2, 3}}; // cut_xy
4238 
4239  const unsigned int equ_std_subface =
4240  equivalent_iso_subface[std_face_ref][std_subface];
4241  Assert(equ_std_subface != e, ExcInternalError());
4242 
4243  return iso_children[ref_case - 1][face][equ_std_subface];
4244  }
4245  else
4246  {
4247  // the face_ref_case was too coarse,
4248  // throw an error
4249  Assert(false,
4250  ExcMessage("The face RefineCase is too coarse "
4251  "for the given cell RefineCase."));
4252  }
4253  // we only get here in case of an error
4254  return e;
4255 }
4256 
4257 
4258 
4259 template <>
4260 inline unsigned int
4262  const unsigned int,
4263  const unsigned int,
4264  const bool,
4265  const bool,
4266  const bool,
4267  const RefinementCase<3> &)
4268 {
4269  Assert(false, ExcNotImplemented());
4271 }
4272 
4273 
4274 
4275 template <>
4276 inline unsigned int
4277 GeometryInfo<1>::face_to_cell_lines(const unsigned int face,
4278  const unsigned int line,
4279  const bool,
4280  const bool,
4281  const bool)
4282 {
4283  (void)face;
4284  (void)line;
4285  Assert(face + 1 < faces_per_cell + 1, ExcIndexRange(face, 0, faces_per_cell));
4286  Assert(line + 1 < lines_per_face + 1, ExcIndexRange(line, 0, lines_per_face));
4287 
4288  // There is only a single line, so
4289  // it must be this.
4290  return 0;
4291 }
4292 
4293 
4294 
4295 template <>
4296 inline unsigned int
4297 GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4298  const unsigned int line,
4299  const bool,
4300  const bool,
4301  const bool)
4302 {
4303  (void)line;
4306 
4307  // The face is a line itself.
4308  return face;
4309 }
4310 
4311 
4312 
4313 template <>
4314 inline unsigned int
4315 GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4316  const unsigned int line,
4317  const bool face_orientation,
4318  const bool face_flip,
4319  const bool face_rotation)
4320 {
4323 
4324  const unsigned lines[faces_per_cell][lines_per_face] = {
4325  {8, 10, 0, 4}, // left face
4326  {9, 11, 1, 5}, // right face
4327  {2, 6, 8, 9}, // front face
4328  {3, 7, 10, 11}, // back face
4329  {0, 1, 2, 3}, // bottom face
4330  {4, 5, 6, 7}}; // top face
4331  return lines[face][real_to_standard_face_line(
4332  line, face_orientation, face_flip, face_rotation)];
4333 }
4334 
4335 
4336 
4337 template <int dim>
4338 inline unsigned int
4339 GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4340  const unsigned int,
4341  const bool,
4342  const bool,
4343  const bool)
4344 {
4345  Assert(false, ExcNotImplemented());
4347 }
4348 
4349 
4350 
4351 template <int dim>
4352 inline unsigned int
4353 GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4354  const unsigned int vertex,
4355  const bool face_orientation,
4356  const bool face_flip,
4357  const bool face_rotation)
4358 {
4360  face,
4361  vertex,
4362  face_orientation,
4363  face_flip,
4364  face_rotation);
4365 }
4366 
4367 
4368 
4369 template <int dim>
4370 inline Point<dim>
4372 {
4373  Point<dim> p = q;
4374  for (unsigned int i = 0; i < dim; i++)
4375  if (p[i] < 0.)
4376  p[i] = 0.;
4377  else if (p[i] > 1.)
4378  p[i] = 1.;
4379 
4380  return p;
4381 }
4382 
4383 
4384 
4385 template <int dim>
4386 inline double
4388 {
4389  double result = 0.0;
4390 
4391  for (unsigned int i = 0; i < dim; i++)
4392  if ((-p[i]) > result)
4393  result = -p[i];
4394  else if ((p[i] - 1.) > result)
4395  result = (p[i] - 1.);
4396 
4397  return result;
4398 }
4399 
4400 
4401 
4402 template <int dim>
4403 inline double
4405  const unsigned int i)
4406 {
4409 
4410  switch (dim)
4411  {
4412  case 1:
4413  {
4414  const double x = xi[0];
4415  switch (i)
4416  {
4417  case 0:
4418  return 1 - x;
4419  case 1:
4420  return x;
4421  }
4422  break;
4423  }
4424 
4425  case 2:
4426  {
4427  const double x = xi[0];
4428  const double y = xi[1];
4429  switch (i)
4430  {
4431  case 0:
4432  return (1 - x) * (1 - y);
4433  case 1:
4434  return x * (1 - y);
4435  case 2:
4436  return (1 - x) * y;
4437  case 3:
4438  return x * y;
4439  }
4440  break;
4441  }
4442 
4443  case 3:
4444  {
4445  const double x = xi[0];
4446  const double y = xi[1];
4447  const double z = xi[2];
4448  switch (i)
4449  {
4450  case 0:
4451  return (1 - x) * (1 - y) * (1 - z);
4452  case 1:
4453  return x * (1 - y) * (1 - z);
4454  case 2:
4455  return (1 - x) * y * (1 - z);
4456  case 3:
4457  return x * y * (1 - z);
4458  case 4:
4459  return (1 - x) * (1 - y) * z;
4460  case 5:
4461  return x * (1 - y) * z;
4462  case 6:
4463  return (1 - x) * y * z;
4464  case 7:
4465  return x * y * z;
4466  }
4467  break;
4468  }
4469 
4470  default:
4471  Assert(false, ExcNotImplemented());
4472  }
4473  return -1e9;
4474 }
4475 
4476 
4477 
4478 template <>
4480  const Point<1> &,
4481  const unsigned int i)
4482 {
4485 
4486  switch (i)
4487  {
4488  case 0:
4489  return Point<1>(-1.);
4490  case 1:
4491  return Point<1>(1.);
4492  }
4493 
4494  return Point<1>(-1e9);
4495 }
4496 
4497 
4498 
4499 template <>
4501  const Point<2> & xi,
4502  const unsigned int i)
4503 {
4506 
4507  const double x = xi[0];
4508  const double y = xi[1];
4509  switch (i)
4510  {
4511  case 0:
4512  return Point<2>(-(1 - y), -(1 - x));
4513  case 1:
4514  return Point<2>(1 - y, -x);
4515  case 2:
4516  return Point<2>(-y, 1 - x);
4517  case 3:
4518  return Point<2>(y, x);
4519  }
4520  return Point<2>(-1e9, -1e9);
4521 }
4522 
4523 
4524 
4525 template <>
4527  const Point<3> & xi,
4528  const unsigned int i)
4529 {
4532 
4533  const double x = xi[0];
4534  const double y = xi[1];
4535  const double z = xi[2];
4536  switch (i)
4537  {
4538  case 0:
4539  return Point<3>(-(1 - y) * (1 - z),
4540  -(1 - x) * (1 - z),
4541  -(1 - x) * (1 - y));
4542  case 1:
4543  return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4544  case 2:
4545  return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4546  case 3:
4547  return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4548  case 4:
4549  return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4550  case 5:
4551  return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4552  case 6:
4553  return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4554  case 7:
4555  return Point<3>(y * z, x * z, x * y);
4556  }
4557 
4558  return Point<3>(-1e9, -1e9, -1e9);
4559 }
4560 
4561 
4562 
4563 template <int dim>
4564 inline Tensor<1, dim>
4566  const unsigned int)
4567 {
4568  Assert(false, ExcNotImplemented());
4569  return Tensor<1, dim>();
4570 }
4571 
4572 
4573 unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4574 {
4575  return 0;
4576 }
4577 
4578 
4579 namespace internal
4580 {
4581  namespace GeometryInfoHelper
4582  {
4583  // wedge product of a single
4584  // vector in 2d: we just have to
4585  // rotate it by 90 degrees to the
4586  // right
4587  inline Tensor<1, 2>
4588  wedge_product(const Tensor<1, 2> (&derivative)[1])
4589  {
4590  Tensor<1, 2> result;
4591  result[0] = derivative[0][1];
4592  result[1] = -derivative[0][0];
4593 
4594  return result;
4595  }
4596 
4597 
4598  // wedge product of 2 vectors in
4599  // 3d is the cross product
4600  inline Tensor<1, 3>
4601  wedge_product(const Tensor<1, 3> (&derivative)[2])
4602  {
4603  return cross_product_3d(derivative[0], derivative[1]);
4604  }
4605 
4606 
4607  // wedge product of dim vectors
4608  // in dim-d: that's the
4609  // determinant of the matrix
4610  template <int dim>
4611  inline Tensor<0, dim>
4612  wedge_product(const Tensor<1, dim> (&derivative)[dim])
4613  {
4614  Tensor<2, dim> jacobian;
4615  for (unsigned int i = 0; i < dim; ++i)
4616  jacobian[i] = derivative[i];
4617 
4618  return determinant(jacobian);
4619  }
4620  } // namespace GeometryInfoHelper
4621 } // namespace internal
4622 
4623 
4624 template <int dim>
4625 template <int spacedim>
4626 inline void
4628 # ifndef DEAL_II_CONSTEXPR_BUG
4629  (const Point<spacedim> (&vertices)[vertices_per_cell],
4631 # else
4632  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms)
4633 # endif
4634 {
4635  // for each of the vertices,
4636  // compute the alternating form
4637  // of the mapped unit
4638  // vectors. consider for
4639  // example the case of a quad
4640  // in spacedim==3: to do so, we
4641  // need to see how the
4642  // infinitesimal vectors
4643  // (d\xi_1,0) and (0,d\xi_2)
4644  // are transformed into
4645  // spacedim-dimensional space
4646  // and then form their cross
4647  // product (i.e. the wedge product
4648  // of two vectors). to this end, note
4649  // that
4650  // \vec x = sum_i \vec v_i phi_i(\vec xi)
4651  // so the transformed vectors are
4652  // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
4653  // and
4654  // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
4655  // which boils down to the columns
4656  // of the 3x2 matrix \grad_\xi x(\xi)
4657  //
4658  // a similar reasoning would
4659  // hold for all dim,spacedim
4660  // pairs -- we only have to
4661  // compute the wedge product of
4662  // the columns of the
4663  // derivatives
4664  for (unsigned int i = 0; i < vertices_per_cell; ++i)
4665  {
4666  Tensor<1, spacedim> derivatives[dim];
4667 
4668  for (unsigned int j = 0; j < vertices_per_cell; ++j)
4669  {
4670  const Tensor<1, dim> grad_phi_j =
4672  for (unsigned int l = 0; l < dim; ++l)
4673  derivatives[l] += vertices[j] * grad_phi_j[l];
4674  }
4675 
4676  forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
4677  }
4678 }
4679 
4680 #endif // DOXYGEN
4681 DEAL_II_NAMESPACE_CLOSE
4682 
4683 #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)
static const unsigned int invalid_unsigned_int
Definition: types.h:187
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 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)
void serialize(Archive &ar, const unsigned int version)
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
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 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 ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
static const std::array< unsigned int, vertices_per_cell > ucd_to_deal
static constexpr unsigned int vertices_per_cell
RefinementCase operator~() const
static double distance_to_unit_cell(const Point< dim > &p)
UpdateFlags operator|(const UpdateFlags f1, const UpdateFlags f2)
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:518
static constexpr unsigned int lines_per_cell
#define Assert(cond, exc)
Definition: exceptions.h:1407
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
static ::ExceptionBase & ExcInvalidCoordinate(double arg1)
static std::size_t memory_consumption()
static unsigned int standard_to_real_face_vertex(const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static 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
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
Definition: tensor.h:422
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)
static unsigned int n_subfaces(const internal::SubfaceCase< dim > &subface_case)
DEAL_II_CONSTEXPR Number determinant(const SymmetricTensor< 2, dim, Number > &)
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 ::ExceptionBase & ExcNotImplemented()
std::uint8_t value
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:564
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
static RefinementCase cut_axis(const unsigned int i)
static ::ExceptionBase & ExcInternalError()