Reference documentation for deal.II version Git c4a4607bd1 2020-01-22 19:33:00 -0500
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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  AssertIndexRange(i, 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  AssertIndexRange(i, 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  AssertIndexRange(i, 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 {
2640 
2641  return Point<1>(static_cast<double>(vertex));
2642 }
2643 
2644 
2645 
2646 template <>
2647 inline Point<2>
2648 GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2649 {
2651 
2652  return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2653 }
2654 
2655 
2656 
2657 template <>
2658 inline Point<3>
2659 GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2660 {
2662 
2663  return {static_cast<double>(vertex % 2),
2664  static_cast<double>(vertex / 2 % 2),
2665  static_cast<double>(vertex / 4)};
2666 }
2667 
2668 
2669 
2670 template <int dim>
2671 inline Point<dim>
2672 GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2673 {
2674  Assert(false, ExcNotImplemented());
2675 
2676  return Point<dim>();
2677 }
2678 
2679 
2680 
2681 template <>
2682 inline unsigned int
2684 {
2685  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2686 
2687  return (p[0] <= 0.5 ? 0 : 1);
2688 }
2689 
2690 
2691 
2692 template <>
2693 inline unsigned int
2695 {
2696  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2697  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2698 
2699  return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
2700 }
2701 
2702 
2703 
2704 template <>
2705 inline unsigned int
2707 {
2708  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2709  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2710  Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
2711 
2712  return (p[0] <= 0.5 ?
2713  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
2714  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
2715 }
2716 
2717 
2718 template <int dim>
2719 inline unsigned int
2721 {
2722  Assert(false, ExcNotImplemented());
2723 
2724  return 0;
2725 }
2726 
2727 
2728 
2729 template <>
2730 inline Point<1>
2732  const unsigned int child_index,
2733  const RefinementCase<1> refine_case)
2734 
2735 {
2736  AssertIndexRange(child_index, 2);
2738  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
2739 
2740  return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
2741 }
2742 
2743 
2744 
2745 template <>
2746 inline Point<2>
2748  const unsigned int child_index,
2749  const RefinementCase<2> refine_case)
2750 
2751 {
2752  AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
2753 
2754  Point<2> point = p;
2755  switch (refine_case)
2756  {
2758  point[0] *= 2.0;
2759  if (child_index == 1)
2760  point[0] -= 1.0;
2761  break;
2763  point[1] *= 2.0;
2764  if (child_index == 1)
2765  point[1] -= 1.0;
2766  break;
2768  point *= 2.0;
2769  point -= unit_cell_vertex(child_index);
2770  break;
2771  default:
2772  Assert(false, ExcInternalError());
2773  }
2774 
2775  return point;
2776 }
2777 
2778 
2779 
2780 template <>
2781 inline Point<3>
2783  const unsigned int child_index,
2784  const RefinementCase<3> refine_case)
2785 
2786 {
2787  AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
2788 
2789  Point<3> point = p;
2790  // there might be a cleverer way to do
2791  // this, but since this function is called
2792  // in very few places for initialization
2793  // purposes only, I don't care at the
2794  // moment
2795  switch (refine_case)
2796  {
2798  point[0] *= 2.0;
2799  if (child_index == 1)
2800  point[0] -= 1.0;
2801  break;
2803  point[1] *= 2.0;
2804  if (child_index == 1)
2805  point[1] -= 1.0;
2806  break;
2808  point[2] *= 2.0;
2809  if (child_index == 1)
2810  point[2] -= 1.0;
2811  break;
2813  point[0] *= 2.0;
2814  point[1] *= 2.0;
2815  if (child_index % 2 == 1)
2816  point[0] -= 1.0;
2817  if (child_index / 2 == 1)
2818  point[1] -= 1.0;
2819  break;
2821  // careful, this is slightly
2822  // different from xy and yz due to
2823  // different internal numbering of
2824  // children!
2825  point[0] *= 2.0;
2826  point[2] *= 2.0;
2827  if (child_index / 2 == 1)
2828  point[0] -= 1.0;
2829  if (child_index % 2 == 1)
2830  point[2] -= 1.0;
2831  break;
2833  point[1] *= 2.0;
2834  point[2] *= 2.0;
2835  if (child_index % 2 == 1)
2836  point[1] -= 1.0;
2837  if (child_index / 2 == 1)
2838  point[2] -= 1.0;
2839  break;
2841  point *= 2.0;
2842  point -= unit_cell_vertex(child_index);
2843  break;
2844  default:
2845  Assert(false, ExcInternalError());
2846  }
2847 
2848  return point;
2849 }
2850 
2851 
2852 
2853 template <int dim>
2854 inline Point<dim>
2856  const Point<dim> & /*p*/,
2857  const unsigned int /*child_index*/,
2858  const RefinementCase<dim> /*refine_case*/)
2859 
2860 {
2861  Assert(false, ExcNotImplemented());
2862  return Point<dim>();
2863 }
2864 
2865 
2866 
2867 template <>
2868 inline Point<1>
2870  const unsigned int child_index,
2871  const RefinementCase<1> refine_case)
2872 
2873 {
2874  AssertIndexRange(child_index, 2);
2876  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
2877 
2878  return (p + unit_cell_vertex(child_index)) * 0.5;
2879 }
2880 
2881 
2882 
2883 template <>
2884 inline Point<3>
2886  const unsigned int child_index,
2887  const RefinementCase<3> refine_case)
2888 
2889 {
2890  AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
2891 
2892  Point<3> point = p;
2893  // there might be a cleverer way to do
2894  // this, but since this function is called
2895  // in very few places for initialization
2896  // purposes only, I don't care at the
2897  // moment
2898  switch (refine_case)
2899  {
2901  if (child_index == 1)
2902  point[0] += 1.0;
2903  point[0] *= 0.5;
2904  break;
2906  if (child_index == 1)
2907  point[1] += 1.0;
2908  point[1] *= 0.5;
2909  break;
2911  if (child_index == 1)
2912  point[2] += 1.0;
2913  point[2] *= 0.5;
2914  break;
2916  if (child_index % 2 == 1)
2917  point[0] += 1.0;
2918  if (child_index / 2 == 1)
2919  point[1] += 1.0;
2920  point[0] *= 0.5;
2921  point[1] *= 0.5;
2922  break;
2924  // careful, this is slightly
2925  // different from xy and yz due to
2926  // different internal numbering of
2927  // children!
2928  if (child_index / 2 == 1)
2929  point[0] += 1.0;
2930  if (child_index % 2 == 1)
2931  point[2] += 1.0;
2932  point[0] *= 0.5;
2933  point[2] *= 0.5;
2934  break;
2936  if (child_index % 2 == 1)
2937  point[1] += 1.0;
2938  if (child_index / 2 == 1)
2939  point[2] += 1.0;
2940  point[1] *= 0.5;
2941  point[2] *= 0.5;
2942  break;
2944  point += unit_cell_vertex(child_index);
2945  point *= 0.5;
2946  break;
2947  default:
2948  Assert(false, ExcInternalError());
2949  }
2950 
2951  return point;
2952 }
2953 
2954 
2955 
2956 template <>
2957 inline Point<2>
2959  const unsigned int child_index,
2960  const RefinementCase<2> refine_case)
2961 {
2962  AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
2963 
2964  Point<2> point = p;
2965  switch (refine_case)
2966  {
2968  if (child_index == 1)
2969  point[0] += 1.0;
2970  point[0] *= 0.5;
2971  break;
2973  if (child_index == 1)
2974  point[1] += 1.0;
2975  point[1] *= 0.5;
2976  break;
2978  point += unit_cell_vertex(child_index);
2979  point *= 0.5;
2980  break;
2981  default:
2982  Assert(false, ExcInternalError());
2983  }
2984 
2985  return point;
2986 }
2987 
2988 
2989 
2990 template <int dim>
2991 inline Point<dim>
2993  const Point<dim> & /*p*/,
2994  const unsigned int /*child_index*/,
2995  const RefinementCase<dim> /*refine_case*/)
2996 {
2997  Assert(false, ExcNotImplemented());
2998  return Point<dim>();
2999 }
3000 
3001 
3002 
3003 template <int dim>
3004 inline bool
3006 {
3007  Assert(false, ExcNotImplemented());
3008  return false;
3009 }
3010 
3011 template <>
3012 inline bool
3014 {
3015  return (p[0] >= 0.) && (p[0] <= 1.);
3016 }
3017 
3018 
3019 
3020 template <>
3021 inline bool
3023 {
3024  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
3025 }
3026 
3027 
3028 
3029 template <>
3030 inline bool
3032 {
3033  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
3034  (p[2] >= 0.) && (p[2] <= 1.);
3035 }
3036 
3037 
3038 
3039 template <int dim>
3040 inline bool
3042 {
3043  Assert(false, ExcNotImplemented());
3044  return false;
3045 }
3046 
3047 template <>
3048 inline bool
3049 GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
3050 {
3051  return (p[0] >= -eps) && (p[0] <= 1. + eps);
3052 }
3053 
3054 
3055 
3056 template <>
3057 inline bool
3058 GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
3059 {
3060  const double l = -eps, u = 1 + eps;
3061  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
3062 }
3063 
3064 
3065 
3066 template <>
3067 inline bool
3068 GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
3069 {
3070  const double l = -eps, u = 1.0 + eps;
3071  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
3072  (p[2] >= l) && (p[2] <= u);
3073 }
3074 
3075 
3076 
3077 template <>
3078 inline unsigned int
3079 GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
3080  const unsigned int vertex)
3081 {
3082  (void)line;
3084  AssertIndexRange(vertex, 2);
3085 
3086  return vertex;
3087 }
3088 
3089 
3090 template <>
3091 inline unsigned int
3092 GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
3093  const unsigned int vertex)
3094 {
3095  constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
3096  return cell_vertices[line][vertex];
3097 }
3098 
3099 
3100 
3101 template <>
3102 inline unsigned int
3103 GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
3104  const unsigned int vertex)
3105 {
3107  AssertIndexRange(vertex, 2);
3108 
3109  constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3110  {1, 3},
3111  {0, 1},
3112  {2, 3},
3113  {4, 6}, // top face
3114  {5, 7},
3115  {4, 5},
3116  {6, 7},
3117  {0,
3118  4}, // connects of bottom
3119  {1, 5}, // top face
3120  {2, 6},
3121  {3, 7}};
3122  return vertices[line][vertex];
3123 }
3124 
3125 
3126 template <>
3127 inline unsigned int
3128 GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3129 {
3130  Assert(false, ExcNotImplemented());
3132 }
3133 
3134 template <>
3135 inline unsigned int
3136 GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3137  const bool face_orientation,
3138  const bool face_flip,
3139  const bool face_rotation)
3140 {
3142 
3143  // set up a table to make sure that
3144  // we handle non-standard faces correctly
3145  //
3146  // so set up a table that for each vertex (of
3147  // a quad in standard position) describes
3148  // which vertex to take
3149  //
3150  // first index: four vertices 0...3
3151  //
3152  // second index: face_orientation; 0:
3153  // opposite normal, 1: standard
3154  //
3155  // third index: face_flip; 0: standard, 1:
3156  // face rotated by 180 degrees
3157  //
3158  // forth index: face_rotation: 0: standard,
3159  // 1: face rotated by 90 degrees
3160 
3161  constexpr unsigned int vertex_translation[4][2][2][2] = {
3162  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3163  // face_rotation=false and true
3164  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3165  // face_rotation=false and true
3166  {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3167  // face_rotation=false and true
3168  {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3169  // face_rotation=false and true
3170 
3171  {{{2, 3}, // vertex 1 ...
3172  {1, 0}},
3173  {{1, 0}, {2, 3}}},
3174 
3175  {{{1, 0}, // vertex 2 ...
3176  {2, 3}},
3177  {{2, 3}, {1, 0}}},
3178 
3179  {{{3, 1}, // vertex 3 ...
3180  {0, 2}},
3181  {{3, 1}, {0, 2}}}};
3182 
3183  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3184 }
3185 
3186 
3187 
3188 template <int dim>
3189 inline unsigned int
3190 GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3191  const bool,
3192  const bool,
3193  const bool)
3194 {
3195  Assert(dim > 1, ExcImpossibleInDim(dim));
3197  return vertex;
3198 }
3199 
3200 template <int dim>
3201 inline unsigned int
3203 {
3204  constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3205  0, 2, 2, 4, 2, 4, 4, 8};
3206 
3207  return n_children[ref_case];
3208 }
3209 
3210 
3211 
3212 template <int dim>
3213 inline unsigned int
3215 {
3216  Assert(false, ExcNotImplemented());
3217  return 0;
3218 }
3219 
3220 template <>
3221 inline unsigned int
3223 {
3224  Assert(false, ExcImpossibleInDim(1));
3225  return 0;
3226 }
3227 
3228 template <>
3229 inline unsigned int
3231 {
3232  return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3233 }
3234 
3235 
3236 
3237 template <>
3238 inline unsigned int
3240 {
3241  const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3242  0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3243  return nsubs[subface_case];
3244 }
3245 
3246 
3247 
3248 template <int dim>
3249 inline double
3251  const unsigned int)
3252 {
3253  Assert(false, ExcNotImplemented());
3254  return 0.;
3255 }
3256 
3257 template <>
3258 inline double
3260  const unsigned int)
3261 {
3262  return 1;
3263 }
3264 
3265 
3266 template <>
3267 inline double
3269  const unsigned int)
3270 {
3271  double ratio = 1;
3272  switch (subface_case)
3273  {
3275  // Here, an
3276  // Assert(false,ExcInternalError())
3277  // would be the right
3278  // choice, but
3279  // unfortunately the
3280  // current function is
3281  // also called for faces
3282  // without children (see
3283  // tests/fe/mapping.cc).
3284  // Assert(false, ExcMessage("Face has no subfaces."));
3285  // Furthermore, assign
3286  // following value as
3287  // otherwise the
3288  // bits/volume_x tests
3289  // break
3291  break;
3293  ratio = 0.5;
3294  break;
3295  default:
3296  // there should be no
3297  // cases left
3298  Assert(false, ExcInternalError());
3299  break;
3300  }
3301 
3302  return ratio;
3303 }
3304 
3305 
3306 template <>
3307 inline double
3309  const unsigned int subface_no)
3310 {
3311  double ratio = 1;
3312  switch (subface_case)
3313  {
3315  // Here, an
3316  // Assert(false,ExcInternalError())
3317  // would be the right
3318  // choice, but
3319  // unfortunately the
3320  // current function is
3321  // also called for faces
3322  // without children (see
3323  // tests/bits/mesh_3d_16.cc). Add
3324  // following switch to
3325  // avoid diffs in
3326  // tests/bits/mesh_3d_16
3328  break;
3331  ratio = 0.5;
3332  break;
3336  ratio = 0.25;
3337  break;
3340  if (subface_no < 2)
3341  ratio = 0.25;
3342  else
3343  ratio = 0.5;
3344  break;
3347  if (subface_no == 0)
3348  ratio = 0.5;
3349  else
3350  ratio = 0.25;
3351  break;
3352  default:
3353  // there should be no
3354  // cases left
3355  Assert(false, ExcInternalError());
3356  break;
3357  }
3358 
3359  return ratio;
3360 }
3361 
3362 
3363 
3364 template <int dim>
3366  const RefinementCase<dim> &,
3367  const unsigned int,
3368  const bool,
3369  const bool,
3370  const bool)
3371 {
3372  Assert(false, ExcNotImplemented());
3373  return RefinementCase<dim - 1>::no_refinement;
3374 }
3375 
3376 template <>
3378  const RefinementCase<1> &,
3379  const unsigned int,
3380  const bool,
3381  const bool,
3382  const bool)
3383 {
3384  Assert(false, ExcImpossibleInDim(1));
3385 
3387 }
3388 
3389 
3390 template <>
3391 inline RefinementCase<1>
3393  const RefinementCase<2> &cell_refinement_case,
3394  const unsigned int face_no,
3395  const bool,
3396  const bool,
3397  const bool)
3398 {
3399  const unsigned int dim = 2;
3400  AssertIndexRange(cell_refinement_case,
3403 
3404  const RefinementCase<dim - 1>
3407  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3408  RefinementCase<dim - 1>::no_refinement},
3409 
3410  {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3411 
3412  {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3413 
3414  {RefinementCase<dim - 1>::cut_x, // cut_xy
3415  RefinementCase<dim - 1>::cut_x}};
3416 
3417  return ref_cases[cell_refinement_case][face_no / 2];
3418 }
3419 
3420 
3421 template <>
3422 inline RefinementCase<2>
3424  const RefinementCase<3> &cell_refinement_case,
3425  const unsigned int face_no,
3426  const bool face_orientation,
3427  const bool /*face_flip*/,
3428  const bool face_rotation)
3429 {
3430  const unsigned int dim = 3;
3431  AssertIndexRange(cell_refinement_case,
3434 
3435  const RefinementCase<dim - 1>
3438  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3439  RefinementCase<dim - 1>::no_refinement,
3440  RefinementCase<dim - 1>::no_refinement},
3441 
3442  {RefinementCase<dim - 1>::no_refinement, // cut_x
3443  RefinementCase<dim - 1>::cut_y,
3444  RefinementCase<dim - 1>::cut_x},
3445 
3446  {RefinementCase<dim - 1>::cut_x, // cut_y
3447  RefinementCase<dim - 1>::no_refinement,
3448  RefinementCase<dim - 1>::cut_y},
3449 
3450  {RefinementCase<dim - 1>::cut_x, // cut_xy
3451  RefinementCase<dim - 1>::cut_y,
3452  RefinementCase<dim - 1>::cut_xy},
3453 
3454  {RefinementCase<dim - 1>::cut_y, // cut_z
3455  RefinementCase<dim - 1>::cut_x,
3456  RefinementCase<dim - 1>::no_refinement},
3457 
3458  {RefinementCase<dim - 1>::cut_y, // cut_xz
3459  RefinementCase<dim - 1>::cut_xy,
3460  RefinementCase<dim - 1>::cut_x},
3461 
3462  {RefinementCase<dim - 1>::cut_xy, // cut_yz
3463  RefinementCase<dim - 1>::cut_x,
3464  RefinementCase<dim - 1>::cut_y},
3465 
3466  {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3467  RefinementCase<dim - 1>::cut_xy,
3468  RefinementCase<dim - 1>::cut_xy},
3469  };
3470 
3471  const RefinementCase<dim - 1> ref_case =
3472  ref_cases[cell_refinement_case][face_no / 2];
3473 
3474  const RefinementCase<dim - 1> flip[4] = {
3475  RefinementCase<dim - 1>::no_refinement,
3476  RefinementCase<dim - 1>::cut_y,
3477  RefinementCase<dim - 1>::cut_x,
3478  RefinementCase<dim - 1>::cut_xy};
3479 
3480  // correct the ref_case for face_orientation
3481  // and face_rotation. for face_orientation,
3482  // 'true' is the default value whereas for
3483  // face_rotation, 'false' is standard. If
3484  // <tt>face_rotation==face_orientation</tt>,
3485  // then one of them is non-standard and we
3486  // have to swap cut_x and cut_y, otherwise no
3487  // change is necessary. face_flip has no
3488  // influence. however, in order to keep the
3489  // interface consistent with other functions,
3490  // we still include it as an argument to this
3491  // function
3492  return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3493 }
3494 
3495 
3496 
3497 template <int dim>
3498 inline RefinementCase<1>
3500  const unsigned int)
3501 {
3502  Assert(false, ExcNotImplemented());
3504 }
3505 
3506 template <>
3507 inline RefinementCase<1>
3509  const RefinementCase<1> &cell_refinement_case,
3510  const unsigned int line_no)
3511 {
3512  (void)line_no;
3513  const unsigned int dim = 1;
3514  (void)dim;
3515  AssertIndexRange(cell_refinement_case,
3518 
3519  return cell_refinement_case;
3520 }
3521 
3522 
3523 template <>
3524 inline RefinementCase<1>
3526  const RefinementCase<2> &cell_refinement_case,
3527  const unsigned int line_no)
3528 {
3529  // Assertions are in face_refinement_case()
3530  return face_refinement_case(cell_refinement_case, line_no);
3531 }
3532 
3533 
3534 template <>
3535 inline RefinementCase<1>
3537  const RefinementCase<3> &cell_refinement_case,
3538  const unsigned int line_no)
3539 {
3540  const unsigned int dim = 3;
3541  AssertIndexRange(cell_refinement_case,
3544 
3545  // array indicating, which simple refine
3546  // case cuts a line in direction x, y or
3547  // z. For example, cut_y and everything
3548  // containing cut_y (cut_xy, cut_yz,
3549  // cut_xyz) cuts lines, which are in y
3550  // direction.
3551  const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3554 
3555  // order the direction of lines
3556  // 0->x, 1->y, 2->z
3557  const unsigned int direction[lines_per_cell] = {
3558  1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3559 
3560  return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3563 }
3564 
3565 
3566 
3567 template <int dim>
3568 inline RefinementCase<dim>
3570  const RefinementCase<dim - 1> &,
3571  const unsigned int,
3572  const bool,
3573  const bool,
3574  const bool)
3575 {
3576  Assert(false, ExcNotImplemented());
3577 
3579 }
3580 
3581 template <>
3582 inline RefinementCase<1>
3584  const RefinementCase<0> &,
3585  const unsigned int,
3586  const bool,
3587  const bool,
3588  const bool)
3589 {
3590  const unsigned int dim = 1;
3591  Assert(false, ExcImpossibleInDim(dim));
3592 
3594 }
3595 
3596 
3597 template <>
3598 inline RefinementCase<2>
3601  const unsigned int face_no,
3602  const bool,
3603  const bool,
3604  const bool)
3605 {
3606  const unsigned int dim = 2;
3607  AssertIndexRange(face_refinement_case,
3610 
3611  if (face_refinement_case == RefinementCase<dim>::cut_x)
3612  return (face_no / 2) ? RefinementCase<dim>::cut_x :
3614  else
3616 }
3617 
3618 
3619 template <>
3620 inline RefinementCase<3>
3622  const RefinementCase<2> &face_refinement_case,
3623  const unsigned int face_no,
3624  const bool face_orientation,
3625  const bool /*face_flip*/,
3626  const bool face_rotation)
3627 {
3628  const unsigned int dim = 3;
3629  AssertIndexRange(face_refinement_case,
3632 
3637 
3638  // correct the face_refinement_case for
3639  // face_orientation and face_rotation. for
3640  // face_orientation, 'true' is the default
3641  // value whereas for face_rotation, 'false'
3642  // is standard. If
3643  // <tt>face_rotation==face_orientation</tt>,
3644  // then one of them is non-standard and we
3645  // have to swap cut_x and cut_y, otherwise no
3646  // change is necessary. face_flip has no
3647  // influence. however, in order to keep the
3648  // interface consistent with other functions,
3649  // we still include it as an argument to this
3650  // function
3651  const RefinementCase<dim - 1> std_face_ref =
3652  (face_orientation == face_rotation) ? flip[face_refinement_case] :
3653  face_refinement_case;
3654 
3655  const RefinementCase<dim> face_to_cell[3][4] = {
3656  {RefinementCase<dim>::no_refinement, // faces 0 and 1
3657  RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3658  RefinementCase<dim>::cut_z,
3660 
3661  {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3662  // "exchanged on faces 2 and 3")
3663  RefinementCase<dim>::cut_z,
3666 
3667  {RefinementCase<dim>::no_refinement, // faces 4 and 5
3671 
3672  return face_to_cell[face_no / 2][std_face_ref];
3673 }
3674 
3675 
3676 
3677 template <int dim>
3678 inline RefinementCase<dim>
3680  const unsigned int)
3681 {
3682  Assert(false, ExcNotImplemented());
3683 
3685 }
3686 
3687 template <>
3688 inline RefinementCase<1>
3690  const unsigned int line_no)
3691 {
3692  (void)line_no;
3693  AssertIndexRange(line_no, 1);
3694 
3695  return RefinementCase<1>::cut_x;
3696 }
3697 
3698 
3699 template <>
3700 inline RefinementCase<2>
3702  const unsigned int line_no)
3703 {
3704  const unsigned int dim = 2;
3705  (void)dim;
3707 
3708  return (line_no / 2) ? RefinementCase<2>::cut_x : RefinementCase<2>::cut_y;
3709 }
3710 
3711 
3712 template <>
3713 inline RefinementCase<3>
3715  const unsigned int line_no)
3716 {
3717  const unsigned int dim = 3;
3719 
3720  const RefinementCase<dim> ref_cases[6] = {
3721  RefinementCase<dim>::cut_y, // lines 0 and 1
3722  RefinementCase<dim>::cut_x, // lines 2 and 3
3723  RefinementCase<dim>::cut_y, // lines 4 and 5
3724  RefinementCase<dim>::cut_x, // lines 6 and 7
3725  RefinementCase<dim>::cut_z, // lines 8 and 9
3726  RefinementCase<dim>::cut_z}; // lines 10 and 11
3727 
3728  return ref_cases[line_no / 2];
3729 }
3730 
3731 
3732 
3733 template <>
3734 inline unsigned int
3735 GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
3736  const bool face_orientation,
3737  const bool face_flip,
3738  const bool face_rotation)
3739 {
3741 
3742  // set up a table to make sure that
3743  // we handle non-standard faces correctly
3744  //
3745  // so set up a table that for each vertex (of
3746  // a quad in standard position) describes
3747  // which vertex to take
3748  //
3749  // first index: four vertices 0...3
3750  //
3751  // second index: face_orientation; 0:
3752  // opposite normal, 1: standard
3753  //
3754  // third index: face_flip; 0: standard, 1:
3755  // face rotated by 180 degrees
3756  //
3757  // forth index: face_rotation: 0: standard,
3758  // 1: face rotated by 90 degrees
3759 
3760  const unsigned int vertex_translation[4][2][2][2] = {
3761  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3762  // face_rotation=false and true
3763  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3764  // face_rotation=false and true
3765  {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
3766  // face_rotation=false and true
3767  {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
3768  // face_rotation=false and true
3769 
3770  {{{2, 3}, // vertex 1 ...
3771  {1, 0}},
3772  {{1, 3}, {2, 0}}},
3773 
3774  {{{1, 0}, // vertex 2 ...
3775  {2, 3}},
3776  {{2, 0}, {1, 3}}},
3777 
3778  {{{3, 1}, // vertex 3 ...
3779  {0, 2}},
3780  {{3, 2}, {0, 1}}}};
3781 
3782  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3783 }
3784 
3785 
3786 
3787 template <int dim>
3788 inline unsigned int
3789 GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
3790  const bool,
3791  const bool,
3792  const bool)
3793 {
3794  Assert(dim > 1, ExcImpossibleInDim(dim));
3796  return vertex;
3797 }
3798 
3799 
3800 
3801 template <>
3802 inline unsigned int
3803 GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
3804  const bool face_orientation,
3805  const bool face_flip,
3806  const bool face_rotation)
3807 {
3809 
3810 
3811  // make sure we handle
3812  // non-standard faces correctly
3813  //
3814  // so set up a table that for each line (of a
3815  // quad) describes which line to take
3816  //
3817  // first index: four lines 0...3
3818  //
3819  // second index: face_orientation; 0:
3820  // opposite normal, 1: standard
3821  //
3822  // third index: face_flip; 0: standard, 1:
3823  // face rotated by 180 degrees
3824  //
3825  // forth index: face_rotation: 0: standard,
3826  // 1: face rotated by 90 degrees
3827 
3828  const unsigned int line_translation[4][2][2][2] = {
3829  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
3830  // face_rotation=false and true
3831  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
3832  // face_rotation=false and true
3833  {{0, 3}, // line 0, face_orientation=true, face_flip=false,
3834  // face_rotation=false and true
3835  {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
3836  // face_rotation=false and true
3837 
3838  {{{3, 1}, // line 1 ...
3839  {2, 0}},
3840  {{1, 2}, {0, 3}}},
3841 
3842  {{{0, 3}, // line 2 ...
3843  {1, 2}},
3844  {{2, 0}, {3, 1}}},
3845 
3846  {{{1, 2}, // line 3 ...
3847  {0, 3}},
3848  {{3, 1}, {2, 0}}}};
3849 
3850  return line_translation[line][face_orientation][face_flip][face_rotation];
3851 }
3852 
3853 
3854 
3855 template <int dim>
3856 inline unsigned int
3857 GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
3858  const bool,
3859  const bool,
3860  const bool)
3861 {
3862  Assert(false, ExcNotImplemented());
3863  return line;
3864 }
3865 
3866 
3867 
3868 template <>
3869 inline unsigned int
3870 GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
3871  const bool face_orientation,
3872  const bool face_flip,
3873  const bool face_rotation)
3874 {
3876 
3877 
3878  // make sure we handle
3879  // non-standard faces correctly
3880  //
3881  // so set up a table that for each line (of a
3882  // quad) describes which line to take
3883  //
3884  // first index: four lines 0...3
3885  //
3886  // second index: face_orientation; 0:
3887  // opposite normal, 1: standard
3888  //
3889  // third index: face_flip; 0: standard, 1:
3890  // face rotated by 180 degrees
3891  //
3892  // forth index: face_rotation: 0: standard,
3893  // 1: face rotated by 90 degrees
3894 
3895  const unsigned int line_translation[4][2][2][2] = {
3896  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
3897  // face_rotation=false and true
3898  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
3899  // face_rotation=false and true
3900  {{0, 2}, // line 0, face_orientation=true, face_flip=false,
3901  // face_rotation=false and true
3902  {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
3903  // face_rotation=false and true
3904 
3905  {{{3, 1}, // line 1 ...
3906  {2, 0}},
3907  {{1, 3}, {0, 2}}},
3908 
3909  {{{0, 3}, // line 2 ...
3910  {1, 2}},
3911  {{2, 1}, {3, 0}}},
3912 
3913  {{{1, 2}, // line 3 ...
3914  {0, 3}},
3915  {{3, 0}, {2, 1}}}};
3916 
3917  return line_translation[line][face_orientation][face_flip][face_rotation];
3918 }
3919 
3920 
3921 
3922 template <int dim>
3923 inline unsigned int
3924 GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
3925  const bool,
3926  const bool,
3927  const bool)
3928 {
3929  Assert(false, ExcNotImplemented());
3930  return line;
3931 }
3932 
3933 
3934 
3935 template <>
3936 inline unsigned int
3938  const unsigned int face,
3939  const unsigned int subface,
3940  const bool,
3941  const bool,
3942  const bool,
3943  const RefinementCase<0> &)
3944 {
3945  (void)subface;
3948 
3949  return face;
3950 }
3951 
3952 
3953 
3954 template <>
3955 inline unsigned int
3957  const unsigned int face,
3958  const unsigned int subface,
3959  const bool /*face_orientation*/,
3960  const bool face_flip,
3961  const bool /*face_rotation*/,
3962  const RefinementCase<1> &)
3963 {
3966 
3967  // always return the child adjacent to the specified
3968  // subface. if the face of a cell is not refined, don't
3969  // throw an assertion but deliver the child adjacent to
3970  // the face nevertheless, i.e. deliver the child of
3971  // this cell adjacent to the subface of a possibly
3972  // refined neighbor. this simplifies setting neighbor
3973  // information in execute_refinement.
3974  const unsigned int
3976  [max_children_per_face] = {
3977  {
3978  // Normal orientation (face_flip = false)
3979  {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
3980  {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
3981  {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
3982  },
3983  {
3984  // Flipped orientation (face_flip = true)
3985  {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
3986  {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
3987  {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
3988  }};
3989 
3990  return subcells[face_flip][ref_case - 1][face][subface];
3991 }
3992 
3993 
3994 
3995 template <>
3996 inline unsigned int
3998  const unsigned int face,
3999  const unsigned int subface,
4000  const bool face_orientation,
4001  const bool face_flip,
4002  const bool face_rotation,
4003  const RefinementCase<2> &face_ref_case)
4004 {
4005  const unsigned int dim = 3;
4006 
4008  ExcMessage("Cell has no children."));
4010  if (!(subface == 0 &&
4011  face_ref_case == RefinementCase<dim - 1>::no_refinement))
4012  {
4013  AssertIndexRange(subface,
4014  GeometryInfo<dim - 1>::n_children(face_ref_case));
4015  }
4016 
4017  // invalid number used for invalid cases,
4018  // e.g. when the children are more refined at
4019  // a given face than the face itself
4020  const unsigned int e = numbers::invalid_unsigned_int;
4021 
4022  // the whole process of finding a child cell
4023  // at a given subface considering the
4024  // possibly anisotropic refinement cases of
4025  // the cell and the face as well as
4026  // orientation, flip and rotation of the face
4027  // is quite complicated. thus, we break it
4028  // down into several steps.
4029 
4030  // first step: convert the given face refine
4031  // case to a face refine case concerning the
4032  // face in standard orientation (, flip and
4033  // rotation). This only affects cut_x and
4034  // cut_y
4035  const RefinementCase<dim - 1> flip[4] = {
4036  RefinementCase<dim - 1>::no_refinement,
4037  RefinementCase<dim - 1>::cut_y,
4038  RefinementCase<dim - 1>::cut_x,
4039  RefinementCase<dim - 1>::cut_xy};
4040  // for face_orientation, 'true' is the
4041  // default value whereas for face_rotation,
4042  // 'false' is standard. If
4043  // <tt>face_rotation==face_orientation</tt>,
4044  // then one of them is non-standard and we
4045  // have to swap cut_x and cut_y, otherwise no
4046  // change is necessary.
4047  const RefinementCase<dim - 1> std_face_ref =
4048  (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4049 
4050  // second step: convert the given subface
4051  // index to the one for a standard face
4052  // respecting face_orientation, face_flip and
4053  // face_rotation
4054 
4055  // first index: face_ref_case
4056  // second index: face_orientation
4057  // third index: face_flip
4058  // forth index: face_rotation
4059  // fifth index: subface index
4060  const unsigned int subface_exchange[4][2][2][2][4] = {
4061  // no_refinement (subface 0 stays 0,
4062  // all others are invalid)
4063  {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4064  {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4065  // cut_x (here, if the face is only
4066  // rotated OR only falsely oriented,
4067  // then subface 0 of the non-standard
4068  // face does NOT correspond to one of
4069  // the subfaces of a standard
4070  // face. Thus we indicate the subface
4071  // which is located at the lower left
4072  // corner (the origin of the face's
4073  // local coordinate system) with
4074  // '0'. The rest of this issue is
4075  // taken care of using the above
4076  // conversion to a 'standard face
4077  // refine case')
4078  {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4079  {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4080  // cut_y (the same applies as for
4081  // cut_x)
4082  {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4083  {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4084  // cut_xyz: this information is
4085  // identical to the information
4086  // returned by
4087  // GeometryInfo<3>::real_to_standard_face_vertex()
4088  {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4089  // face_rotation=false, subfaces 0,1,2,3
4090  {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4091  // face_rotation=true, subfaces 0,1,2,3
4092  {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4093  // face_rotation=false, subfaces 0,1,2,3
4094  {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4095  // face_rotation=true, subfaces 0,1,2,3
4096  {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4097  // face_rotation=false, subfaces 0,1,2,3
4098  {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4099  // face_rotation=true, subfaces 0,1,2,3
4100  {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4101  // face_rotation=false, subfaces 0,1,2,3
4102  {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4103  // face_rotation=true, subfaces 0,1,2,3
4104 
4105  const unsigned int std_subface =
4106  subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4107  [subface];
4108  Assert(std_subface != e, ExcInternalError());
4109 
4110  // third step: these are the children, which
4111  // can be found at the given subfaces of an
4112  // isotropically refined (standard) face
4113  //
4114  // first index: (refinement_case-1)
4115  // second index: face_index
4116  // third index: subface_index (isotropic refinement)
4117  const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4118  [max_children_per_face] = {
4119  // cut_x
4120  {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4121  {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4122  {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4123  {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4124  {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4125  {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4126  // cut_y
4127  {{0, 1, 0, 1},
4128  {0, 1, 0, 1},
4129  {0, 0, 0, 0},
4130  {1, 1, 1, 1},
4131  {0, 0, 1, 1},
4132  {0, 0, 1, 1}},
4133  // cut_xy
4134  {{0, 2, 0, 2},
4135  {1, 3, 1, 3},
4136  {0, 0, 1, 1},
4137  {2, 2, 3, 3},
4138  {0, 1, 2, 3},
4139  {0, 1, 2, 3}},
4140  // cut_z
4141  {{0, 0, 1, 1},
4142  {0, 0, 1, 1},
4143  {0, 1, 0, 1},
4144  {0, 1, 0, 1},
4145  {0, 0, 0, 0},
4146  {1, 1, 1, 1}},
4147  // cut_xz
4148  {{0, 0, 1, 1},
4149  {2, 2, 3, 3},
4150  {0, 1, 2, 3},
4151  {0, 1, 2, 3},
4152  {0, 2, 0, 2},
4153  {1, 3, 1, 3}},
4154  // cut_yz
4155  {{0, 1, 2, 3},
4156  {0, 1, 2, 3},
4157  {0, 2, 0, 2},
4158  {1, 3, 1, 3},
4159  {0, 0, 1, 1},
4160  {2, 2, 3, 3}},
4161  // cut_xyz
4162  {{0, 2, 4, 6},
4163  {1, 3, 5, 7},
4164  {0, 4, 1, 5},
4165  {2, 6, 3, 7},
4166  {0, 1, 2, 3},
4167  {4, 5, 6, 7}}};
4168 
4169  // forth step: check, whether the given face
4170  // refine case is valid for the given cell
4171  // refine case. this is the case, if the
4172  // given face refine case is at least as
4173  // refined as the face is for the given cell
4174  // refine case
4175 
4176  // note, that we are considering standard
4177  // face refinement cases here and thus must
4178  // not pass the given orientation, flip and
4179  // rotation flags
4180  if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4181  face_refinement_case(ref_case, face))
4182  {
4183  // all is fine. for anisotropic face
4184  // refine cases, select one of the
4185  // isotropic subfaces which neighbors the
4186  // same child
4187 
4188  // first index: (standard) face refine case
4189  // second index: subface index
4190  const unsigned int equivalent_iso_subface[4][4] = {
4191  {0, e, e, e}, // no_refinement
4192  {0, 3, e, e}, // cut_x
4193  {0, 3, e, e}, // cut_y
4194  {0, 1, 2, 3}}; // cut_xy
4195 
4196  const unsigned int equ_std_subface =
4197  equivalent_iso_subface[std_face_ref][std_subface];
4198  Assert(equ_std_subface != e, ExcInternalError());
4199 
4200  return iso_children[ref_case - 1][face][equ_std_subface];
4201  }
4202  else
4203  {
4204  // the face_ref_case was too coarse,
4205  // throw an error
4206  Assert(false,
4207  ExcMessage("The face RefineCase is too coarse "
4208  "for the given cell RefineCase."));
4209  }
4210  // we only get here in case of an error
4211  return e;
4212 }
4213 
4214 
4215 
4216 template <>
4217 inline unsigned int
4219  const unsigned int,
4220  const unsigned int,
4221  const bool,
4222  const bool,
4223  const bool,
4224  const RefinementCase<3> &)
4225 {
4226  Assert(false, ExcNotImplemented());
4228 }
4229 
4230 
4231 
4232 template <>
4233 inline unsigned int
4234 GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4235  const unsigned int line,
4236  const bool,
4237  const bool,
4238  const bool)
4239 {
4240  (void)line;
4243 
4244  // The face is a line itself.
4245  return face;
4246 }
4247 
4248 
4249 
4250 template <>
4251 inline unsigned int
4252 GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4253  const unsigned int line,
4254  const bool face_orientation,
4255  const bool face_flip,
4256  const bool face_rotation)
4257 {
4260 
4261  const unsigned lines[faces_per_cell][lines_per_face] = {
4262  {8, 10, 0, 4}, // left face
4263  {9, 11, 1, 5}, // right face
4264  {2, 6, 8, 9}, // front face
4265  {3, 7, 10, 11}, // back face
4266  {0, 1, 2, 3}, // bottom face
4267  {4, 5, 6, 7}}; // top face
4268  return lines[face][real_to_standard_face_line(
4269  line, face_orientation, face_flip, face_rotation)];
4270 }
4271 
4272 
4273 
4274 template <int dim>
4275 inline unsigned int
4276 GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4277  const unsigned int,
4278  const bool,
4279  const bool,
4280  const bool)
4281 {
4282  Assert(false, ExcNotImplemented());
4284 }
4285 
4286 
4287 
4288 template <int dim>
4289 inline unsigned int
4290 GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4291  const unsigned int vertex,
4292  const bool face_orientation,
4293  const bool face_flip,
4294  const bool face_rotation)
4295 {
4297  face,
4298  vertex,
4299  face_orientation,
4300  face_flip,
4301  face_rotation);
4302 }
4303 
4304 
4305 
4306 template <int dim>
4307 inline Point<dim>
4309 {
4310  Point<dim> p = q;
4311  for (unsigned int i = 0; i < dim; i++)
4312  if (p[i] < 0.)
4313  p[i] = 0.;
4314  else if (p[i] > 1.)
4315  p[i] = 1.;
4316 
4317  return p;
4318 }
4319 
4320 
4321 
4322 template <int dim>
4323 inline double
4325 {
4326  double result = 0.0;
4327 
4328  for (unsigned int i = 0; i < dim; i++)
4329  if ((-p[i]) > result)
4330  result = -p[i];
4331  else if ((p[i] - 1.) > result)
4332  result = (p[i] - 1.);
4333 
4334  return result;
4335 }
4336 
4337 
4338 
4339 template <int dim>
4340 inline double
4342  const unsigned int i)
4343 {
4345 
4346  switch (dim)
4347  {
4348  case 1:
4349  {
4350  const double x = xi[0];
4351  switch (i)
4352  {
4353  case 0:
4354  return 1 - x;
4355  case 1:
4356  return x;
4357  }
4358  break;
4359  }
4360 
4361  case 2:
4362  {
4363  const double x = xi[0];
4364  const double y = xi[1];
4365  switch (i)
4366  {
4367  case 0:
4368  return (1 - x) * (1 - y);
4369  case 1:
4370  return x * (1 - y);
4371  case 2:
4372  return (1 - x) * y;
4373  case 3:
4374  return x * y;
4375  }
4376  break;
4377  }
4378 
4379  case 3:
4380  {
4381  const double x = xi[0];
4382  const double y = xi[1];
4383  const double z = xi[2];
4384  switch (i)
4385  {
4386  case 0:
4387  return (1 - x) * (1 - y) * (1 - z);
4388  case 1:
4389  return x * (1 - y) * (1 - z);
4390  case 2:
4391  return (1 - x) * y * (1 - z);
4392  case 3:
4393  return x * y * (1 - z);
4394  case 4:
4395  return (1 - x) * (1 - y) * z;
4396  case 5:
4397  return x * (1 - y) * z;
4398  case 6:
4399  return (1 - x) * y * z;
4400  case 7:
4401  return x * y * z;
4402  }
4403  break;
4404  }
4405 
4406  default:
4407  Assert(false, ExcNotImplemented());
4408  }
4409  return -1e9;
4410 }
4411 
4412 
4413 
4414 template <>
4416  const Point<1> &,
4417  const unsigned int i)
4418 {
4420 
4421  switch (i)
4422  {
4423  case 0:
4424  return Point<1>(-1.);
4425  case 1:
4426  return Point<1>(1.);
4427  }
4428 
4429  return Point<1>(-1e9);
4430 }
4431 
4432 
4433 
4434 template <>
4436  const Point<2> & xi,
4437  const unsigned int i)
4438 {
4440 
4441  const double x = xi[0];
4442  const double y = xi[1];
4443  switch (i)
4444  {
4445  case 0:
4446  return Point<2>(-(1 - y), -(1 - x));
4447  case 1:
4448  return Point<2>(1 - y, -x);
4449  case 2:
4450  return Point<2>(-y, 1 - x);
4451  case 3:
4452  return Point<2>(y, x);
4453  }
4454  return Point<2>(-1e9, -1e9);
4455 }
4456 
4457 
4458 
4459 template <>
4461  const Point<3> & xi,
4462  const unsigned int i)
4463 {
4465 
4466  const double x = xi[0];
4467  const double y = xi[1];
4468  const double z = xi[2];
4469  switch (i)
4470  {
4471  case 0:
4472  return Point<3>(-(1 - y) * (1 - z),
4473  -(1 - x) * (1 - z),
4474  -(1 - x) * (1 - y));
4475  case 1:
4476  return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4477  case 2:
4478  return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4479  case 3:
4480  return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4481  case 4:
4482  return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4483  case 5:
4484  return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4485  case 6:
4486  return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4487  case 7:
4488  return Point<3>(y * z, x * z, x * y);
4489  }
4490 
4491  return Point<3>(-1e9, -1e9, -1e9);
4492 }
4493 
4494 
4495 
4496 template <int dim>
4497 inline Tensor<1, dim>
4499  const unsigned int)
4500 {
4501  Assert(false, ExcNotImplemented());
4502  return Tensor<1, dim>();
4503 }
4504 
4505 
4506 unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4507 {
4508  return 0;
4509 }
4510 
4511 
4512 namespace internal
4513 {
4514  namespace GeometryInfoHelper
4515  {
4516  // wedge product of a single
4517  // vector in 2d: we just have to
4518  // rotate it by 90 degrees to the
4519  // right
4520  inline Tensor<1, 2>
4521  wedge_product(const Tensor<1, 2> (&derivative)[1])
4522  {
4523  Tensor<1, 2> result;
4524  result[0] = derivative[0][1];
4525  result[1] = -derivative[0][0];
4526 
4527  return result;
4528  }
4529 
4530 
4531  // wedge product of 2 vectors in
4532  // 3d is the cross product
4533  inline Tensor<1, 3>
4534  wedge_product(const Tensor<1, 3> (&derivative)[2])
4535  {
4536  return cross_product_3d(derivative[0], derivative[1]);
4537  }
4538 
4539 
4540  // wedge product of dim vectors
4541  // in dim-d: that's the
4542  // determinant of the matrix
4543  template <int dim>
4544  inline Tensor<0, dim>
4545  wedge_product(const Tensor<1, dim> (&derivative)[dim])
4546  {
4547  Tensor<2, dim> jacobian;
4548  for (unsigned int i = 0; i < dim; ++i)
4549  jacobian[i] = derivative[i];
4550 
4551  return determinant(jacobian);
4552  }
4553  } // namespace GeometryInfoHelper
4554 } // namespace internal
4555 
4556 
4557 template <int dim>
4558 template <int spacedim>
4559 inline void
4561 # ifndef DEAL_II_CONSTEXPR_BUG
4562  (const Point<spacedim> (&vertices)[vertices_per_cell],
4564 # else
4565  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms)
4566 # endif
4567 {
4568  // for each of the vertices,
4569  // compute the alternating form
4570  // of the mapped unit
4571  // vectors. consider for
4572  // example the case of a quad
4573  // in spacedim==3: to do so, we
4574  // need to see how the
4575  // infinitesimal vectors
4576  // (d\xi_1,0) and (0,d\xi_2)
4577  // are transformed into
4578  // spacedim-dimensional space
4579  // and then form their cross
4580  // product (i.e. the wedge product
4581  // of two vectors). to this end, note
4582  // that
4583  // \vec x = sum_i \vec v_i phi_i(\vec xi)
4584  // so the transformed vectors are
4585  // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
4586  // and
4587  // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
4588  // which boils down to the columns
4589  // of the 3x2 matrix \grad_\xi x(\xi)
4590  //
4591  // a similar reasoning would
4592  // hold for all dim,spacedim
4593  // pairs -- we only have to
4594  // compute the wedge product of
4595  // the columns of the
4596  // derivatives
4597  for (unsigned int i = 0; i < vertices_per_cell; ++i)
4598  {
4599  Tensor<1, spacedim> derivatives[dim];
4600 
4601  for (unsigned int j = 0; j < vertices_per_cell; ++j)
4602  {
4603  const Tensor<1, dim> grad_phi_j =
4605  for (unsigned int l = 0; l < dim; ++l)
4606  derivatives[l] += vertices[j] * grad_phi_j[l];
4607  }
4608 
4609  forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
4610  }
4611 }
4612 
4613 #endif // DOXYGEN
4614 DEAL_II_NAMESPACE_CLOSE
4615 
4616 #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)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1641
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 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: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
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()