Reference documentation for deal.II version 9.1.0-pre
geometry_info.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2018 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> ucd_to_deal{{0, 1}};
45 
46  static constexpr std::array<unsigned int, 2> unit_normal_direction{
47  {0, 0}};
48 
49  static constexpr std::array<int, 2> unit_normal_orientation{{-1, 1}};
50 
51  static constexpr std::array<unsigned int, 2> opposite_face{{1, 0}};
52 
53  static constexpr std::array<unsigned int, 2> dx_to_deal{{0, 1}};
54 
55  static constexpr std::array<std::array<unsigned int, 1>, 2>
56  vertex_to_face{{{{0}}, {{1}}}};
57  };
58 
59  template <>
60  struct Initializers<2>
61  {
62  static constexpr std::array<unsigned int, 4> ucd_to_deal{{0, 1, 3, 2}};
63 
64  static constexpr std::array<unsigned int, 4> unit_normal_direction{
65  {0, 0, 1, 1}};
66 
67  static constexpr std::array<int, 4> unit_normal_orientation{
68  {-1, 1, -1, 1}};
69 
70  static constexpr std::array<unsigned int, 4> opposite_face{{1, 0, 3, 2}};
71 
72  static constexpr std::array<unsigned int, 4> dx_to_deal{{0, 2, 1, 3}};
73 
74  static constexpr std::array<std::array<unsigned int, 2>, 4>
75  vertex_to_face{{{{0, 2}}, {{1, 2}}, {{0, 3}}, {{1, 3}}}};
76  };
77 
78  template <>
79  struct Initializers<3>
80  {
81  static constexpr std::array<unsigned int, 8> ucd_to_deal{
82  {0, 1, 5, 4, 2, 3, 7, 6}};
83 
84  static constexpr std::array<unsigned int, 6> unit_normal_direction{
85  {0, 0, 1, 1, 2, 2}};
86 
87  static constexpr std::array<int, 6> unit_normal_orientation{
88  {-1, 1, -1, 1, -1, 1}};
89 
90  static constexpr std::array<unsigned int, 6> opposite_face{
91  {1, 0, 3, 2, 5, 4}};
92 
93  static constexpr std::array<unsigned int, 8> dx_to_deal{
94  {0, 4, 2, 6, 1, 5, 3, 7}};
95 
96  static constexpr std::array<std::array<unsigned int, 3>, 8>
97  vertex_to_face{{{{0, 2, 4}},
98  {{1, 2, 4}},
99  {{0, 3, 4}},
100  {{1, 3, 4}},
101  {{0, 2, 5}},
102  {{1, 2, 5}},
103  {{0, 3, 5}},
104  {{1, 3, 5}}}};
105  };
106 
107  template <>
108  struct Initializers<4>
109  {
110  static constexpr std::array<unsigned int, 16> ucd_to_deal{
111  {numbers::invalid_unsigned_int,
112  numbers::invalid_unsigned_int,
113  numbers::invalid_unsigned_int,
114  numbers::invalid_unsigned_int,
115  numbers::invalid_unsigned_int,
116  numbers::invalid_unsigned_int,
117  numbers::invalid_unsigned_int,
118  numbers::invalid_unsigned_int,
119  numbers::invalid_unsigned_int,
120  numbers::invalid_unsigned_int,
121  numbers::invalid_unsigned_int,
122  numbers::invalid_unsigned_int,
123  numbers::invalid_unsigned_int,
124  numbers::invalid_unsigned_int,
125  numbers::invalid_unsigned_int,
126  numbers::invalid_unsigned_int}};
127 
128  static constexpr std::array<unsigned int, 8> unit_normal_direction{
129  {0, 0, 1, 1, 2, 2, 3, 3}};
130 
131  static constexpr std::array<int, 8> unit_normal_orientation{
132  {-1, 1, -1, 1, -1, 1, -1, 1}};
133 
134  static constexpr std::array<unsigned int, 8> opposite_face{
135  {1, 0, 3, 2, 5, 4, 7, 6}};
136 
137  static constexpr std::array<unsigned int, 16> dx_to_deal{
138  {numbers::invalid_unsigned_int,
139  numbers::invalid_unsigned_int,
140  numbers::invalid_unsigned_int,
141  numbers::invalid_unsigned_int,
142  numbers::invalid_unsigned_int,
143  numbers::invalid_unsigned_int,
144  numbers::invalid_unsigned_int,
145  numbers::invalid_unsigned_int,
146  numbers::invalid_unsigned_int,
147  numbers::invalid_unsigned_int,
148  numbers::invalid_unsigned_int,
149  numbers::invalid_unsigned_int,
150  numbers::invalid_unsigned_int,
151  numbers::invalid_unsigned_int,
152  numbers::invalid_unsigned_int,
153  numbers::invalid_unsigned_int}};
154 
155  static constexpr std::array<std::array<unsigned int, 4>, 16>
156  vertex_to_face{{{{numbers::invalid_unsigned_int,
157  numbers::invalid_unsigned_int,
158  numbers::invalid_unsigned_int,
159  numbers::invalid_unsigned_int}},
160  {{numbers::invalid_unsigned_int,
161  numbers::invalid_unsigned_int,
162  numbers::invalid_unsigned_int,
163  numbers::invalid_unsigned_int}},
164  {{numbers::invalid_unsigned_int,
165  numbers::invalid_unsigned_int,
166  numbers::invalid_unsigned_int,
167  numbers::invalid_unsigned_int}},
168  {{numbers::invalid_unsigned_int,
169  numbers::invalid_unsigned_int,
170  numbers::invalid_unsigned_int,
171  numbers::invalid_unsigned_int}},
172  {{numbers::invalid_unsigned_int,
173  numbers::invalid_unsigned_int,
174  numbers::invalid_unsigned_int,
175  numbers::invalid_unsigned_int}},
176  {{numbers::invalid_unsigned_int,
177  numbers::invalid_unsigned_int,
178  numbers::invalid_unsigned_int,
179  numbers::invalid_unsigned_int}},
180  {{numbers::invalid_unsigned_int,
181  numbers::invalid_unsigned_int,
182  numbers::invalid_unsigned_int,
183  numbers::invalid_unsigned_int}},
184  {{numbers::invalid_unsigned_int,
185  numbers::invalid_unsigned_int,
186  numbers::invalid_unsigned_int,
187  numbers::invalid_unsigned_int}},
188  {{numbers::invalid_unsigned_int,
189  numbers::invalid_unsigned_int,
190  numbers::invalid_unsigned_int,
191  numbers::invalid_unsigned_int}},
192  {{numbers::invalid_unsigned_int,
193  numbers::invalid_unsigned_int,
194  numbers::invalid_unsigned_int,
195  numbers::invalid_unsigned_int}},
196  {{numbers::invalid_unsigned_int,
197  numbers::invalid_unsigned_int,
198  numbers::invalid_unsigned_int,
199  numbers::invalid_unsigned_int}},
200  {{numbers::invalid_unsigned_int,
201  numbers::invalid_unsigned_int,
202  numbers::invalid_unsigned_int,
203  numbers::invalid_unsigned_int}},
204  {{numbers::invalid_unsigned_int,
205  numbers::invalid_unsigned_int,
206  numbers::invalid_unsigned_int,
207  numbers::invalid_unsigned_int}},
208  {{numbers::invalid_unsigned_int,
209  numbers::invalid_unsigned_int,
210  numbers::invalid_unsigned_int,
211  numbers::invalid_unsigned_int}},
212  {{numbers::invalid_unsigned_int,
213  numbers::invalid_unsigned_int,
214  numbers::invalid_unsigned_int,
215  numbers::invalid_unsigned_int}},
216  {{numbers::invalid_unsigned_int,
217  numbers::invalid_unsigned_int,
218  numbers::invalid_unsigned_int,
219  numbers::invalid_unsigned_int}}}};
220  };
221  } // namespace GeometryInfoHelper
222 } // namespace internal
223 #endif // DOXYGEN
224 
225 
243 class GeometryPrimitive
244 {
245 public:
252  enum Object
253  {
257  vertex = 0,
261  line = 1,
265  quad = 2,
269  hex = 3
270  };
271 
276  GeometryPrimitive(const Object object);
277 
283  GeometryPrimitive(const unsigned int object_dimension);
284 
289  operator unsigned int() const;
290 
291 private:
295  Object object;
296 };
297 
298 
299 
313 template <int dim>
314 struct RefinementPossibilities
315 {
350  enum Possibilities
351  {
355  no_refinement = 0,
356 
360  isotropic_refinement = static_cast<std::uint8_t>(-1)
361  };
362 };
363 
364 
365 
376 template <>
377 struct RefinementPossibilities<1>
378 {
413  enum Possibilities
414  {
418  no_refinement = 0,
422  cut_x = 1,
426  isotropic_refinement = cut_x
427  };
428 };
429 
430 
431 
443 template <>
444 struct RefinementPossibilities<2>
445 {
480  enum Possibilities
481  {
485  no_refinement = 0,
489  cut_x = 1,
493  cut_y = 2,
497  cut_xy = cut_x | cut_y,
498 
502  isotropic_refinement = cut_xy
503  };
504 };
505 
506 
507 
519 template <>
520 struct RefinementPossibilities<3>
521 {
556  enum Possibilities
557  {
561  no_refinement = 0,
565  cut_x = 1,
569  cut_y = 2,
573  cut_xy = cut_x | cut_y,
577  cut_z = 4,
581  cut_xz = cut_x | cut_z,
585  cut_yz = cut_y | cut_z,
589  cut_xyz = cut_x | cut_y | cut_z,
590 
594  isotropic_refinement = cut_xyz
595  };
596 };
597 
598 
599 
611 template <int dim>
612 class RefinementCase : public RefinementPossibilities<dim>
613 {
614 public:
618  RefinementCase();
619 
624  RefinementCase(
625  const typename RefinementPossibilities<dim>::Possibilities refinement_case);
626 
632  explicit RefinementCase(const std::uint8_t refinement_case);
633 
645  operator std::uint8_t() const;
646 
651  RefinementCase
652  operator|(const RefinementCase &r) const;
653 
658  RefinementCase operator&(const RefinementCase &r) const;
659 
667  RefinementCase
668  operator~() const;
669 
670 
676  static RefinementCase
677  cut_axis(const unsigned int i);
678 
682  static std::size_t
683  memory_consumption();
684 
689  template <class Archive>
690  void
691  serialize(Archive &ar, const unsigned int version);
692 
696  DeclException1(
697  ExcInvalidRefinementCase,
698  int,
699  << "The refinement flags given (" << arg1
700  << ") contain set bits that do not "
701  << "make sense for the space dimension of the object to which they are applied.");
702 
703 private:
708  std::uint8_t value : (dim > 0 ? dim : 1);
709 };
710 
711 
712 namespace internal
713 {
734  template <int dim>
735  struct SubfacePossibilities
736  {
740  enum Possibilities
741  {
745  case_none = 0,
746 
750  case_isotropic = static_cast<std::uint8_t>(-1)
751  };
752  };
753 
754 
764  template <>
765  struct SubfacePossibilities<0>
766  {
772  enum Possibilities
773  {
777  case_none = 0,
778 
782  case_isotropic = case_none
783  };
784  };
785 
786 
787 
798  template <>
799  struct SubfacePossibilities<1>
800  {
806  enum Possibilities
807  {
811  case_none = 0,
812 
816  case_isotropic = case_none
817  };
818  };
819 
820 
821 
833  template <>
834  struct SubfacePossibilities<2>
835  {
842  enum Possibilities
843  {
847  case_none = 0,
851  case_x = 1,
855  case_isotropic = case_x
856  };
857  };
858 
859 
860 
953  template <>
954  struct SubfacePossibilities<3>
955  {
962  enum Possibilities
963  {
964  case_none = 0,
965  case_x = 1,
966  case_x1y = 2,
967  case_x2y = 3,
968  case_x1y2y = 4,
969  case_y = 5,
970  case_y1x = 6,
971  case_y2x = 7,
972  case_y1x2x = 8,
973  case_xy = 9,
974 
975  case_isotropic = case_xy
976  };
977  };
978 
979 
980 
988  template <int dim>
989  class SubfaceCase : public SubfacePossibilities<dim>
990  {
991  public:
997  SubfaceCase(const typename SubfacePossibilities<dim>::Possibilities
998  subface_possibility);
999 
1011  operator std::uint8_t() const;
1012 
1016  static constexpr std::size_t
1017  memory_consumption();
1018 
1022  DeclException1(
1023  ExcInvalidSubfaceCase,
1024  int,
1025  << "The subface case given (" << arg1 << ") does not make sense "
1026  << "for the space dimension of the object to which they are applied.");
1027 
1028  private:
1033  std::uint8_t value : (dim == 3 ? 4 : 1);
1034  };
1035 
1036 } // namespace internal
1037 
1038 
1039 
1040 template <int dim>
1041 struct GeometryInfo;
1042 
1043 
1044 
1062 template <>
1063 struct GeometryInfo<0>
1064 {
1072  static constexpr unsigned int max_children_per_cell = 1;
1073 
1077  static constexpr unsigned int faces_per_cell = 0;
1078 
1086  static constexpr unsigned int max_children_per_face = 0;
1087 
1093  static unsigned int
1094  n_children(const RefinementCase<0> &refinement_case);
1095 
1099  static constexpr unsigned int vertices_per_cell = 1;
1100 
1107  static constexpr unsigned int vertices_per_face = 0;
1108 
1112  static constexpr unsigned int lines_per_face = 0;
1113 
1117  static constexpr unsigned int quads_per_face = 0;
1118 
1122  static constexpr unsigned int lines_per_cell = 0;
1123 
1127  static constexpr unsigned int quads_per_cell = 0;
1128 
1132  static constexpr unsigned int hexes_per_cell = 0;
1133 
1151  static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1152 
1166  static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1167 };
1168 
1169 
1170 
1695 template <int dim>
1696 struct GeometryInfo
1697 {
1705  static constexpr unsigned int max_children_per_cell = 1 << dim;
1706 
1710  static constexpr unsigned int faces_per_cell = 2 * dim;
1711 
1719  static constexpr unsigned int max_children_per_face =
1720  GeometryInfo<dim - 1>::max_children_per_cell;
1721 
1725  static constexpr unsigned int vertices_per_cell = 1 << dim;
1726 
1730  static constexpr unsigned int vertices_per_face =
1731  GeometryInfo<dim - 1>::vertices_per_cell;
1732 
1736  static constexpr unsigned int lines_per_face =
1737  GeometryInfo<dim - 1>::lines_per_cell;
1738 
1742  static constexpr unsigned int quads_per_face =
1743  GeometryInfo<dim - 1>::quads_per_cell;
1744 
1754  static constexpr unsigned int lines_per_cell =
1755  (2 * GeometryInfo<dim - 1>::lines_per_cell +
1756  GeometryInfo<dim - 1>::vertices_per_cell);
1757 
1765  static constexpr unsigned int quads_per_cell =
1766  (2 * GeometryInfo<dim - 1>::quads_per_cell +
1767  GeometryInfo<dim - 1>::lines_per_cell);
1768 
1772  static constexpr unsigned int hexes_per_cell =
1773  (2 * GeometryInfo<dim - 1>::hexes_per_cell +
1774  GeometryInfo<dim - 1>::quads_per_cell);
1775 
1793  static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
1794  internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal;
1795 
1809  static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
1810  internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal;
1811 
1822  static constexpr std::array<std::array<unsigned int, dim>, vertices_per_cell>
1823  vertex_to_face =
1824  internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face;
1825 
1830  static unsigned int
1831  n_children(const RefinementCase<dim> &refinement_case);
1832 
1837  static unsigned int
1838  n_subfaces(const internal::SubfaceCase<dim> &subface_case);
1839 
1849  static double
1850  subface_ratio(const internal::SubfaceCase<dim> &subface_case,
1851  const unsigned int subface_no);
1852 
1858  static RefinementCase<dim - 1>
1859  face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
1860  const unsigned int face_no,
1861  const bool face_orientation = true,
1862  const bool face_flip = false,
1863  const bool face_rotation = false);
1864 
1870  static RefinementCase<dim>
1871  min_cell_refinement_case_for_face_refinement(
1872  const RefinementCase<dim - 1> &face_refinement_case,
1873  const unsigned int face_no,
1874  const bool face_orientation = true,
1875  const bool face_flip = false,
1876  const bool face_rotation = false);
1877 
1882  static RefinementCase<1>
1883  line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
1884  const unsigned int line_no);
1885 
1890  static RefinementCase<dim>
1891  min_cell_refinement_case_for_line_refinement(const unsigned int line_no);
1892 
1939  static unsigned int
1940  child_cell_on_face(const RefinementCase<dim> & ref_case,
1941  const unsigned int face,
1942  const unsigned int subface,
1943  const bool face_orientation = true,
1944  const bool face_flip = false,
1945  const bool face_rotation = false,
1946  const RefinementCase<dim - 1> &face_refinement_case =
1947  RefinementCase<dim - 1>::isotropic_refinement);
1948 
1962  static unsigned int
1963  line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
1964 
1985  static unsigned int
1986  face_to_cell_vertices(const unsigned int face,
1987  const unsigned int vertex,
1988  const bool face_orientation = true,
1989  const bool face_flip = false,
1990  const bool face_rotation = false);
1991 
2003  static unsigned int
2004  face_to_cell_lines(const unsigned int face,
2005  const unsigned int line,
2006  const bool face_orientation = true,
2007  const bool face_flip = false,
2008  const bool face_rotation = false);
2009 
2019  static unsigned int
2020  standard_to_real_face_vertex(const unsigned int vertex,
2021  const bool face_orientation = true,
2022  const bool face_flip = false,
2023  const bool face_rotation = false);
2024 
2034  static unsigned int
2035  real_to_standard_face_vertex(const unsigned int vertex,
2036  const bool face_orientation = true,
2037  const bool face_flip = false,
2038  const bool face_rotation = false);
2039 
2049  static unsigned int
2050  standard_to_real_face_line(const unsigned int line,
2051  const bool face_orientation = true,
2052  const bool face_flip = false,
2053  const bool face_rotation = false);
2054 
2064  static unsigned int
2065  real_to_standard_face_line(const unsigned int line,
2066  const bool face_orientation = true,
2067  const bool face_flip = false,
2068  const bool face_rotation = false);
2069 
2075  static Point<dim>
2076  unit_cell_vertex(const unsigned int vertex);
2077 
2087  static unsigned int
2088  child_cell_from_point(const Point<dim> &p);
2089 
2097  static Point<dim>
2098  cell_to_child_coordinates(const Point<dim> & p,
2099  const unsigned int child_index,
2100  const RefinementCase<dim> refine_case =
2101  RefinementCase<dim>::isotropic_refinement);
2102 
2108  static Point<dim>
2109  child_to_cell_coordinates(const Point<dim> & p,
2110  const unsigned int child_index,
2111  const RefinementCase<dim> refine_case =
2112  RefinementCase<dim>::isotropic_refinement);
2113 
2118  static bool
2119  is_inside_unit_cell(const Point<dim> &p);
2120 
2132  static bool
2133  is_inside_unit_cell(const Point<dim> &p, const double eps);
2134 
2139  static Point<dim>
2140  project_to_unit_cell(const Point<dim> &p);
2141 
2147  static double
2148  distance_to_unit_cell(const Point<dim> &p);
2149 
2154  static double
2155  d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2156 
2161  static Tensor<1, dim>
2162  d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2163 
2215  template <int spacedim>
2216  static void
2217  alternating_form_at_vertices
2218 #ifndef DEAL_II_CONSTEXPR_BUG
2219  (const Point<spacedim> (&vertices)[vertices_per_cell],
2220  Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2221 #else
2222  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2223 #endif
2224 
2234  static constexpr std::array<unsigned int, faces_per_cell>
2235  unit_normal_direction =
2236  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction;
2237 
2254  static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2255  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation;
2256 
2262  static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2263  internal::GeometryInfoHelper::Initializers<dim>::opposite_face;
2264 
2265 
2269  DeclException1(ExcInvalidCoordinate,
2270  double,
2271  << "The coordinates must satisfy 0 <= x_i <= 1, "
2272  << "but here we have x_i=" << arg1);
2273 
2277  DeclException3(ExcInvalidSubface,
2278  int,
2279  int,
2280  int,
2281  << "RefinementCase<dim> " << arg1 << ": face " << arg2
2282  << " has no subface " << arg3);
2283 };
2284 
2285 
2286 
2287 #ifndef DOXYGEN
2288 
2289 
2290 /* -------------- declaration of explicit specializations ------------- */
2291 
2292 
2293 template <>
2294 Tensor<1, 1>
2295 GeometryInfo<1>::d_linear_shape_function_gradient(const Point<1> & xi,
2296  const unsigned int i);
2297 template <>
2298 Tensor<1, 2>
2299 GeometryInfo<2>::d_linear_shape_function_gradient(const Point<2> & xi,
2300  const unsigned int i);
2301 template <>
2302 Tensor<1, 3>
2303 GeometryInfo<3>::d_linear_shape_function_gradient(const Point<3> & xi,
2304  const unsigned int i);
2305 
2306 
2307 
2308 /* -------------- inline functions ------------- */
2309 
2310 
2311 inline GeometryPrimitive::GeometryPrimitive(const Object object)
2312  : object(object)
2313 {}
2314 
2315 
2316 
2317 inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2318  : object(static_cast<Object>(object_dimension))
2319 {}
2320 
2321 
2322 inline GeometryPrimitive::operator unsigned int() const
2323 {
2324  return static_cast<unsigned int>(object);
2325 }
2326 
2327 
2328 
2329 namespace internal
2330 {
2331  template <int dim>
2332  inline SubfaceCase<dim>::SubfaceCase(
2333  const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2334  : value(subface_possibility)
2335  {}
2336 
2337 
2338  template <int dim>
2339  inline SubfaceCase<dim>::operator std::uint8_t() const
2340  {
2341  return value;
2342  }
2343 
2344 
2345 } // namespace internal
2346 
2347 
2348 template <int dim>
2349 inline RefinementCase<dim>
2350 RefinementCase<dim>::cut_axis(const unsigned int)
2351 {
2352  Assert(false, ExcInternalError());
2353  return static_cast<std::uint8_t>(-1);
2354 }
2355 
2356 
2357 template <>
2358 inline RefinementCase<1>
2359 RefinementCase<1>::cut_axis(const unsigned int i)
2360 {
2361  Assert(i < 1, ExcIndexRange(i, 0, 1));
2362 
2363  const RefinementCase options[1] = {RefinementPossibilities<1>::cut_x};
2364  return options[i];
2365 }
2366 
2367 
2368 
2369 template <>
2370 inline RefinementCase<2>
2371 RefinementCase<2>::cut_axis(const unsigned int i)
2372 {
2373  Assert(i < 2, ExcIndexRange(i, 0, 2));
2374 
2375  const RefinementCase options[2] = {RefinementPossibilities<2>::cut_x,
2376  RefinementPossibilities<2>::cut_y};
2377  return options[i];
2378 }
2379 
2380 
2381 
2382 template <>
2383 inline RefinementCase<3>
2384 RefinementCase<3>::cut_axis(const unsigned int i)
2385 {
2386  Assert(i < 3, ExcIndexRange(i, 0, 3));
2387 
2388  const RefinementCase options[3] = {RefinementPossibilities<3>::cut_x,
2389  RefinementPossibilities<3>::cut_y,
2390  RefinementPossibilities<3>::cut_z};
2391  return options[i];
2392 }
2393 
2394 
2395 
2396 template <int dim>
2397 inline RefinementCase<dim>::RefinementCase()
2398  : value(RefinementPossibilities<dim>::no_refinement)
2399 {}
2400 
2401 
2402 
2403 template <int dim>
2404 inline RefinementCase<dim>::RefinementCase(
2405  const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2406  : value(refinement_case)
2407 {
2408  // check that only those bits of
2409  // the given argument are set that
2410  // make sense for a given space
2411  // dimension
2412  Assert((refinement_case &
2413  RefinementPossibilities<dim>::isotropic_refinement) ==
2414  refinement_case,
2415  ExcInvalidRefinementCase(refinement_case));
2416 }
2417 
2418 
2419 
2420 template <int dim>
2421 inline RefinementCase<dim>::RefinementCase(const std::uint8_t refinement_case)
2422  : value(refinement_case)
2423 {
2424  // check that only those bits of
2425  // the given argument are set that
2426  // make sense for a given space
2427  // dimension
2428  Assert((refinement_case &
2429  RefinementPossibilities<dim>::isotropic_refinement) ==
2430  refinement_case,
2431  ExcInvalidRefinementCase(refinement_case));
2432 }
2433 
2434 
2435 
2436 template <int dim>
2437 inline RefinementCase<dim>::operator std::uint8_t() const
2438 {
2439  return value;
2440 }
2441 
2442 
2443 
2444 template <int dim>
2445 inline RefinementCase<dim>
2446 RefinementCase<dim>::operator|(const RefinementCase<dim> &r) const
2447 {
2448  return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2449 }
2450 
2451 
2452 
2453 template <int dim>
2454 inline RefinementCase<dim> RefinementCase<dim>::
2455  operator&(const RefinementCase<dim> &r) const
2456 {
2457  return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2458 }
2459 
2460 
2461 
2462 template <int dim>
2463 inline RefinementCase<dim>
2464 RefinementCase<dim>::operator~() const
2465 {
2466  return RefinementCase<dim>(static_cast<std::uint8_t>(
2467  (~value) & RefinementPossibilities<dim>::isotropic_refinement));
2468 }
2469 
2470 
2471 
2472 template <int dim>
2473 inline std::size_t
2474 RefinementCase<dim>::memory_consumption()
2475 {
2476  return sizeof(RefinementCase<dim>);
2477 }
2478 
2479 
2480 
2481 template <int dim>
2482 template <class Archive>
2483 inline void
2484 RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2485 {
2486  // serialization can't deal with bitfields, so copy from/to a full sized
2487  // std::uint8_t
2488  std::uint8_t uchar_value = value;
2489  ar & uchar_value;
2490  value = uchar_value;
2491 }
2492 
2493 
2494 
2495 template <>
2496 inline Point<1>
2497 GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2498 {
2499  Assert(vertex < vertices_per_cell,
2500  ExcIndexRange(vertex, 0, vertices_per_cell));
2501 
2502  return Point<1>(static_cast<double>(vertex));
2503 }
2504 
2505 
2506 
2507 template <>
2508 inline Point<2>
2509 GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2510 {
2511  Assert(vertex < vertices_per_cell,
2512  ExcIndexRange(vertex, 0, vertices_per_cell));
2513 
2514  return Point<2>(vertex % 2, vertex / 2);
2515 }
2516 
2517 
2518 
2519 template <>
2520 inline Point<3>
2521 GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2522 {
2523  Assert(vertex < vertices_per_cell,
2524  ExcIndexRange(vertex, 0, vertices_per_cell));
2525 
2526  return Point<3>(vertex % 2, vertex / 2 % 2, vertex / 4);
2527 }
2528 
2529 
2530 
2531 template <int dim>
2532 inline Point<dim>
2533 GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2534 {
2535  Assert(false, ExcNotImplemented());
2536 
2537  return Point<dim>();
2538 }
2539 
2540 
2541 
2542 template <>
2543 inline unsigned int
2544 GeometryInfo<1>::child_cell_from_point(const Point<1> &p)
2545 {
2546  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2547 
2548  return (p[0] <= 0.5 ? 0 : 1);
2549 }
2550 
2551 
2552 
2553 template <>
2554 inline unsigned int
2555 GeometryInfo<2>::child_cell_from_point(const Point<2> &p)
2556 {
2557  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2558  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2559 
2560  return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
2561 }
2562 
2563 
2564 
2565 template <>
2566 inline unsigned int
2567 GeometryInfo<3>::child_cell_from_point(const Point<3> &p)
2568 {
2569  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2570  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2571  Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
2572 
2573  return (p[0] <= 0.5 ?
2574  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
2575  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
2576 }
2577 
2578 
2579 template <int dim>
2580 inline unsigned int
2581 GeometryInfo<dim>::child_cell_from_point(const Point<dim> &)
2582 {
2583  Assert(false, ExcNotImplemented());
2584 
2585  return 0;
2586 }
2587 
2588 
2589 
2590 template <>
2591 inline Point<1>
2592 GeometryInfo<1>::cell_to_child_coordinates(const Point<1> & p,
2593  const unsigned int child_index,
2594  const RefinementCase<1> refine_case)
2595 
2596 {
2597  Assert(child_index < 2, ExcIndexRange(child_index, 0, 2));
2598  Assert(refine_case == RefinementCase<1>::cut_x, ExcInternalError());
2599  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
2600 
2601  return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
2602 }
2603 
2604 
2605 
2606 template <>
2607 inline Point<2>
2608 GeometryInfo<2>::cell_to_child_coordinates(const Point<2> & p,
2609  const unsigned int child_index,
2610  const RefinementCase<2> refine_case)
2611 
2612 {
2613  Assert(child_index < GeometryInfo<2>::n_children(refine_case),
2614  ExcIndexRange(child_index,
2615  0,
2616  GeometryInfo<2>::n_children(refine_case)));
2617 
2618  Point<2> point = p;
2619  switch (refine_case)
2620  {
2621  case RefinementCase<2>::cut_x:
2622  point[0] *= 2.0;
2623  if (child_index == 1)
2624  point[0] -= 1.0;
2625  break;
2626  case RefinementCase<2>::cut_y:
2627  point[1] *= 2.0;
2628  if (child_index == 1)
2629  point[1] -= 1.0;
2630  break;
2631  case RefinementCase<2>::cut_xy:
2632  point *= 2.0;
2633  point -= unit_cell_vertex(child_index);
2634  break;
2635  default:
2636  Assert(false, ExcInternalError());
2637  }
2638 
2639  return point;
2640 }
2641 
2642 
2643 
2644 template <>
2645 inline Point<3>
2646 GeometryInfo<3>::cell_to_child_coordinates(const Point<3> & p,
2647  const unsigned int child_index,
2648  const RefinementCase<3> refine_case)
2649 
2650 {
2651  Assert(child_index < GeometryInfo<3>::n_children(refine_case),
2652  ExcIndexRange(child_index,
2653  0,
2654  GeometryInfo<3>::n_children(refine_case)));
2655 
2656  Point<3> point = p;
2657  // there might be a cleverer way to do
2658  // this, but since this function is called
2659  // in very few places for initialization
2660  // purposes only, I don't care at the
2661  // moment
2662  switch (refine_case)
2663  {
2664  case RefinementCase<3>::cut_x:
2665  point[0] *= 2.0;
2666  if (child_index == 1)
2667  point[0] -= 1.0;
2668  break;
2669  case RefinementCase<3>::cut_y:
2670  point[1] *= 2.0;
2671  if (child_index == 1)
2672  point[1] -= 1.0;
2673  break;
2674  case RefinementCase<3>::cut_z:
2675  point[2] *= 2.0;
2676  if (child_index == 1)
2677  point[2] -= 1.0;
2678  break;
2679  case RefinementCase<3>::cut_xy:
2680  point[0] *= 2.0;
2681  point[1] *= 2.0;
2682  if (child_index % 2 == 1)
2683  point[0] -= 1.0;
2684  if (child_index / 2 == 1)
2685  point[1] -= 1.0;
2686  break;
2687  case RefinementCase<3>::cut_xz:
2688  // careful, this is slightly
2689  // different from xy and yz due to
2690  // different internal numbering of
2691  // children!
2692  point[0] *= 2.0;
2693  point[2] *= 2.0;
2694  if (child_index / 2 == 1)
2695  point[0] -= 1.0;
2696  if (child_index % 2 == 1)
2697  point[2] -= 1.0;
2698  break;
2699  case RefinementCase<3>::cut_yz:
2700  point[1] *= 2.0;
2701  point[2] *= 2.0;
2702  if (child_index % 2 == 1)
2703  point[1] -= 1.0;
2704  if (child_index / 2 == 1)
2705  point[2] -= 1.0;
2706  break;
2707  case RefinementCase<3>::cut_xyz:
2708  point *= 2.0;
2709  point -= unit_cell_vertex(child_index);
2710  break;
2711  default:
2712  Assert(false, ExcInternalError());
2713  }
2714 
2715  return point;
2716 }
2717 
2718 
2719 
2720 template <int dim>
2721 inline Point<dim>
2722 GeometryInfo<dim>::cell_to_child_coordinates(
2723  const Point<dim> & /*p*/,
2724  const unsigned int /*child_index*/,
2725  const RefinementCase<dim> /*refine_case*/)
2726 
2727 {
2728  Assert(false, ExcNotImplemented());
2729  return Point<dim>();
2730 }
2731 
2732 
2733 
2734 template <>
2735 inline Point<1>
2736 GeometryInfo<1>::child_to_cell_coordinates(const Point<1> & p,
2737  const unsigned int child_index,
2738  const RefinementCase<1> refine_case)
2739 
2740 {
2741  Assert(child_index < 2, ExcIndexRange(child_index, 0, 2));
2742  Assert(refine_case == RefinementCase<1>::cut_x, ExcInternalError());
2743  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
2744 
2745  return (p + unit_cell_vertex(child_index)) * 0.5;
2746 }
2747 
2748 
2749 
2750 template <>
2751 inline Point<3>
2752 GeometryInfo<3>::child_to_cell_coordinates(const Point<3> & p,
2753  const unsigned int child_index,
2754  const RefinementCase<3> refine_case)
2755 
2756 {
2757  Assert(child_index < GeometryInfo<3>::n_children(refine_case),
2758  ExcIndexRange(child_index,
2759  0,
2760  GeometryInfo<3>::n_children(refine_case)));
2761 
2762  Point<3> point = p;
2763  // there might be a cleverer way to do
2764  // this, but since this function is called
2765  // in very few places for initialization
2766  // purposes only, I don't care at the
2767  // moment
2768  switch (refine_case)
2769  {
2770  case RefinementCase<3>::cut_x:
2771  if (child_index == 1)
2772  point[0] += 1.0;
2773  point[0] *= 0.5;
2774  break;
2775  case RefinementCase<3>::cut_y:
2776  if (child_index == 1)
2777  point[1] += 1.0;
2778  point[1] *= 0.5;
2779  break;
2780  case RefinementCase<3>::cut_z:
2781  if (child_index == 1)
2782  point[2] += 1.0;
2783  point[2] *= 0.5;
2784  break;
2785  case RefinementCase<3>::cut_xy:
2786  if (child_index % 2 == 1)
2787  point[0] += 1.0;
2788  if (child_index / 2 == 1)
2789  point[1] += 1.0;
2790  point[0] *= 0.5;
2791  point[1] *= 0.5;
2792  break;
2793  case RefinementCase<3>::cut_xz:
2794  // careful, this is slightly
2795  // different from xy and yz due to
2796  // different internal numbering of
2797  // children!
2798  if (child_index / 2 == 1)
2799  point[0] += 1.0;
2800  if (child_index % 2 == 1)
2801  point[2] += 1.0;
2802  point[0] *= 0.5;
2803  point[2] *= 0.5;
2804  break;
2805  case RefinementCase<3>::cut_yz:
2806  if (child_index % 2 == 1)
2807  point[1] += 1.0;
2808  if (child_index / 2 == 1)
2809  point[2] += 1.0;
2810  point[1] *= 0.5;
2811  point[2] *= 0.5;
2812  break;
2813  case RefinementCase<3>::cut_xyz:
2814  point += unit_cell_vertex(child_index);
2815  point *= 0.5;
2816  break;
2817  default:
2818  Assert(false, ExcInternalError());
2819  }
2820 
2821  return point;
2822 }
2823 
2824 
2825 
2826 template <>
2827 inline Point<2>
2828 GeometryInfo<2>::child_to_cell_coordinates(const Point<2> & p,
2829  const unsigned int child_index,
2830  const RefinementCase<2> refine_case)
2831 {
2832  Assert(child_index < GeometryInfo<2>::n_children(refine_case),
2833  ExcIndexRange(child_index,
2834  0,
2835  GeometryInfo<2>::n_children(refine_case)));
2836 
2837  Point<2> point = p;
2838  switch (refine_case)
2839  {
2840  case RefinementCase<2>::cut_x:
2841  if (child_index == 1)
2842  point[0] += 1.0;
2843  point[0] *= 0.5;
2844  break;
2845  case RefinementCase<2>::cut_y:
2846  if (child_index == 1)
2847  point[1] += 1.0;
2848  point[1] *= 0.5;
2849  break;
2850  case RefinementCase<2>::cut_xy:
2851  point += unit_cell_vertex(child_index);
2852  point *= 0.5;
2853  break;
2854  default:
2855  Assert(false, ExcInternalError());
2856  }
2857 
2858  return point;
2859 }
2860 
2861 
2862 
2863 template <int dim>
2864 inline Point<dim>
2865 GeometryInfo<dim>::child_to_cell_coordinates(
2866  const Point<dim> & /*p*/,
2867  const unsigned int /*child_index*/,
2868  const RefinementCase<dim> /*refine_case*/)
2869 {
2870  Assert(false, ExcNotImplemented());
2871  return Point<dim>();
2872 }
2873 
2874 
2875 
2876 template <int dim>
2877 inline bool
2878 GeometryInfo<dim>::is_inside_unit_cell(const Point<dim> &)
2879 {
2880  Assert(false, ExcNotImplemented());
2881  return false;
2882 }
2883 
2884 template <>
2885 inline bool
2886 GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p)
2887 {
2888  return (p[0] >= 0.) && (p[0] <= 1.);
2889 }
2890 
2891 
2892 
2893 template <>
2894 inline bool
2895 GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p)
2896 {
2897  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
2898 }
2899 
2900 
2901 
2902 template <>
2903 inline bool
2904 GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p)
2905 {
2906  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
2907  (p[2] >= 0.) && (p[2] <= 1.);
2908 }
2909 
2910 
2911 
2912 template <int dim>
2913 inline bool
2914 GeometryInfo<dim>::is_inside_unit_cell(const Point<dim> &, const double)
2915 {
2916  Assert(false, ExcNotImplemented());
2917  return false;
2918 }
2919 
2920 template <>
2921 inline bool
2922 GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
2923 {
2924  return (p[0] >= -eps) && (p[0] <= 1. + eps);
2925 }
2926 
2927 
2928 
2929 template <>
2930 inline bool
2931 GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
2932 {
2933  const double l = -eps, u = 1 + eps;
2934  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
2935 }
2936 
2937 
2938 
2939 template <>
2940 inline bool
2941 GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
2942 {
2943  const double l = -eps, u = 1.0 + eps;
2944  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
2945  (p[2] >= l) && (p[2] <= u);
2946 }
2947 
2948 
2949 
2950 template <>
2951 inline unsigned int
2952 GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
2953  const unsigned int vertex)
2954 {
2955  (void)line;
2956  Assert(line < lines_per_cell, ExcIndexRange(line, 0, lines_per_cell));
2957  Assert(vertex < 2, ExcIndexRange(vertex, 0, 2));
2958 
2959  return vertex;
2960 }
2961 
2962 
2963 template <>
2964 inline unsigned int
2965 GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
2966  const unsigned int vertex)
2967 {
2968  constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
2969  return cell_vertices[line][vertex];
2970 }
2971 
2972 
2973 
2974 template <>
2975 inline unsigned int
2976 GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
2977  const unsigned int vertex)
2978 {
2979  Assert(line < lines_per_cell, ExcIndexRange(line, 0, lines_per_cell));
2980  Assert(vertex < 2, ExcIndexRange(vertex, 0, 2));
2981 
2982  constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
2983  {1, 3},
2984  {0, 1},
2985  {2, 3},
2986  {4, 6}, // top face
2987  {5, 7},
2988  {4, 5},
2989  {6, 7},
2990  {0,
2991  4}, // connects of bottom
2992  {1, 5}, // top face
2993  {2, 6},
2994  {3, 7}};
2995  return vertices[line][vertex];
2996 }
2997 
2998 
2999 template <>
3000 inline unsigned int
3001 GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3002 {
3003  Assert(false, ExcNotImplemented());
3004  return numbers::invalid_unsigned_int;
3005 }
3006 
3007 template <>
3008 inline unsigned int
3009 GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3010  const bool face_orientation,
3011  const bool face_flip,
3012  const bool face_rotation)
3013 {
3014  Assert(vertex < GeometryInfo<3>::vertices_per_face,
3015  ExcIndexRange(vertex, 0, GeometryInfo<3>::vertices_per_face));
3016 
3017  // set up a table to make sure that
3018  // we handle non-standard faces correctly
3019  //
3020  // so set up a table that for each vertex (of
3021  // a quad in standard position) describes
3022  // which vertex to take
3023  //
3024  // first index: four vertices 0...3
3025  //
3026  // second index: face_orientation; 0:
3027  // opposite normal, 1: standard
3028  //
3029  // third index: face_flip; 0: standard, 1:
3030  // face rotated by 180 degrees
3031  //
3032  // forth index: face_rotation: 0: standard,
3033  // 1: face rotated by 90 degrees
3034 
3035  constexpr unsigned int vertex_translation[4][2][2][2] = {
3036  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3037  // face_rotation=false and true
3038  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3039  // face_rotation=false and true
3040  {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3041  // face_rotation=false and true
3042  {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3043  // face_rotation=false and true
3044 
3045  {{{2, 3}, // vertex 1 ...
3046  {1, 0}},
3047  {{1, 0}, {2, 3}}},
3048 
3049  {{{1, 0}, // vertex 2 ...
3050  {2, 3}},
3051  {{2, 3}, {1, 0}}},
3052 
3053  {{{3, 1}, // vertex 3 ...
3054  {0, 2}},
3055  {{3, 1}, {0, 2}}}};
3056 
3057  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3058 }
3059 
3060 
3061 
3062 template <int dim>
3063 inline unsigned int
3064 GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3065  const bool,
3066  const bool,
3067  const bool)
3068 {
3069  Assert(dim > 1, ExcImpossibleInDim(dim));
3070  Assert(vertex < GeometryInfo<dim>::vertices_per_face,
3071  ExcIndexRange(vertex, 0, GeometryInfo<dim>::vertices_per_face));
3072  return vertex;
3073 }
3074 
3075 template <int dim>
3076 inline unsigned int
3077 GeometryInfo<dim>::n_children(const RefinementCase<dim> &ref_case)
3078 {
3079  constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3080  0, 2, 2, 4, 2, 4, 4, 8};
3081 
3082  return n_children[ref_case];
3083 }
3084 
3085 
3086 
3087 template <int dim>
3088 inline unsigned int
3089 GeometryInfo<dim>::n_subfaces(const internal::SubfaceCase<dim> &)
3090 {
3091  Assert(false, ExcNotImplemented());
3092  return 0;
3093 }
3094 
3095 template <>
3096 inline unsigned int
3097 GeometryInfo<1>::n_subfaces(const internal::SubfaceCase<1> &)
3098 {
3099  Assert(false, ExcImpossibleInDim(1));
3100  return 0;
3101 }
3102 
3103 template <>
3104 inline unsigned int
3105 GeometryInfo<2>::n_subfaces(const internal::SubfaceCase<2> &subface_case)
3106 {
3107  return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3108 }
3109 
3110 
3111 
3112 template <>
3113 inline unsigned int
3114 GeometryInfo<3>::n_subfaces(const internal::SubfaceCase<3> &subface_case)
3115 {
3116  const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3117  0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3118  return nsubs[subface_case];
3119 }
3120 
3121 
3122 
3123 template <int dim>
3124 inline double
3125 GeometryInfo<dim>::subface_ratio(const internal::SubfaceCase<dim> &,
3126  const unsigned int)
3127 {
3128  Assert(false, ExcNotImplemented());
3129  return 0.;
3130 }
3131 
3132 template <>
3133 inline double
3134 GeometryInfo<1>::subface_ratio(const internal::SubfaceCase<1> &,
3135  const unsigned int)
3136 {
3137  return 1;
3138 }
3139 
3140 
3141 template <>
3142 inline double
3143 GeometryInfo<2>::subface_ratio(const internal::SubfaceCase<2> &subface_case,
3144  const unsigned int)
3145 {
3146  double ratio = 1;
3147  switch (subface_case)
3148  {
3149  case internal::SubfaceCase<2>::case_none:
3150  // Here, an
3151  // Assert(false,ExcInternalError())
3152  // would be the right
3153  // choice, but
3154  // unfortunately the
3155  // current function is
3156  // also called for faces
3157  // without children (see
3158  // tests/fe/mapping.cc).
3159  // Assert(false, ExcMessage("Face has no subfaces."));
3160  // Furthermore, assign
3161  // following value as
3162  // otherwise the
3163  // bits/volume_x tests
3164  // break
3165  ratio = 1. / GeometryInfo<2>::max_children_per_face;
3166  break;
3167  case internal::SubfaceCase<2>::case_x:
3168  ratio = 0.5;
3169  break;
3170  default:
3171  // there should be no
3172  // cases left
3173  Assert(false, ExcInternalError());
3174  break;
3175  }
3176 
3177  return ratio;
3178 }
3179 
3180 
3181 template <>
3182 inline double
3183 GeometryInfo<3>::subface_ratio(const internal::SubfaceCase<3> &subface_case,
3184  const unsigned int subface_no)
3185 {
3186  double ratio = 1;
3187  switch (subface_case)
3188  {
3189  case internal::SubfaceCase<3>::case_none:
3190  // Here, an
3191  // Assert(false,ExcInternalError())
3192  // would be the right
3193  // choice, but
3194  // unfortunately the
3195  // current function is
3196  // also called for faces
3197  // without children (see
3198  // tests/bits/mesh_3d_16.cc). Add
3199  // following switch to
3200  // avoid diffs in
3201  // tests/bits/mesh_3d_16
3202  ratio = 1. / GeometryInfo<3>::max_children_per_face;
3203  break;
3204  case internal::SubfaceCase<3>::case_x:
3205  case internal::SubfaceCase<3>::case_y:
3206  ratio = 0.5;
3207  break;
3208  case internal::SubfaceCase<3>::case_xy:
3209  case internal::SubfaceCase<3>::case_x1y2y:
3210  case internal::SubfaceCase<3>::case_y1x2x:
3211  ratio = 0.25;
3212  break;
3213  case internal::SubfaceCase<3>::case_x1y:
3214  case internal::SubfaceCase<3>::case_y1x:
3215  if (subface_no < 2)
3216  ratio = 0.25;
3217  else
3218  ratio = 0.5;
3219  break;
3220  case internal::SubfaceCase<3>::case_x2y:
3221  case internal::SubfaceCase<3>::case_y2x:
3222  if (subface_no == 0)
3223  ratio = 0.5;
3224  else
3225  ratio = 0.25;
3226  break;
3227  default:
3228  // there should be no
3229  // cases left
3230  Assert(false, ExcInternalError());
3231  break;
3232  }
3233 
3234  return ratio;
3235 }
3236 
3237 
3238 
3239 template <int dim>
3240 RefinementCase<dim - 1> inline GeometryInfo<dim>::face_refinement_case(
3241  const RefinementCase<dim> &,
3242  const unsigned int,
3243  const bool,
3244  const bool,
3245  const bool)
3246 {
3247  Assert(false, ExcNotImplemented());
3248  return RefinementCase<dim - 1>::no_refinement;
3249 }
3250 
3251 template <>
3252 RefinementCase<0> inline GeometryInfo<1>::face_refinement_case(
3253  const RefinementCase<1> &,
3254  const unsigned int,
3255  const bool,
3256  const bool,
3257  const bool)
3258 {
3259  Assert(false, ExcImpossibleInDim(1));
3260 
3261  return RefinementCase<0>::no_refinement;
3262 }
3263 
3264 
3265 template <>
3266 inline RefinementCase<1>
3267 GeometryInfo<2>::face_refinement_case(
3268  const RefinementCase<2> &cell_refinement_case,
3269  const unsigned int face_no,
3270  const bool,
3271  const bool,
3272  const bool)
3273 {
3274  const unsigned int dim = 2;
3275  Assert(cell_refinement_case < RefinementCase<dim>::isotropic_refinement + 1,
3276  ExcIndexRange(cell_refinement_case,
3277  0,
3278  RefinementCase<dim>::isotropic_refinement + 1));
3279  Assert(face_no < GeometryInfo<dim>::faces_per_cell,
3280  ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
3281 
3282  const RefinementCase<dim - 1>
3283  ref_cases[RefinementCase<dim>::isotropic_refinement +
3284  1][GeometryInfo<dim>::faces_per_cell / 2] = {
3285  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3286  RefinementCase<dim - 1>::no_refinement},
3287 
3288  {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3289 
3290  {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3291 
3292  {RefinementCase<dim - 1>::cut_x, // cut_xy
3293  RefinementCase<dim - 1>::cut_x}};
3294 
3295  return ref_cases[cell_refinement_case][face_no / 2];
3296 }
3297 
3298 
3299 template <>
3300 inline RefinementCase<2>
3301 GeometryInfo<3>::face_refinement_case(
3302  const RefinementCase<3> &cell_refinement_case,
3303  const unsigned int face_no,
3304  const bool face_orientation,
3305  const bool /*face_flip*/,
3306  const bool face_rotation)
3307 {
3308  const unsigned int dim = 3;
3309  Assert(cell_refinement_case < RefinementCase<dim>::isotropic_refinement + 1,
3310  ExcIndexRange(cell_refinement_case,
3311  0,
3312  RefinementCase<dim>::isotropic_refinement + 1));
3313  Assert(face_no < GeometryInfo<dim>::faces_per_cell,
3314  ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
3315 
3316  const RefinementCase<dim - 1>
3317  ref_cases[RefinementCase<dim>::isotropic_refinement + 1]
3318  [GeometryInfo<dim>::faces_per_cell / 2] = {
3319  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3320  RefinementCase<dim - 1>::no_refinement,
3321  RefinementCase<dim - 1>::no_refinement},
3322 
3323  {RefinementCase<dim - 1>::no_refinement, // cut_x
3324  RefinementCase<dim - 1>::cut_y,
3325  RefinementCase<dim - 1>::cut_x},
3326 
3327  {RefinementCase<dim - 1>::cut_x, // cut_y
3328  RefinementCase<dim - 1>::no_refinement,
3329  RefinementCase<dim - 1>::cut_y},
3330 
3331  {RefinementCase<dim - 1>::cut_x, // cut_xy
3332  RefinementCase<dim - 1>::cut_y,
3333  RefinementCase<dim - 1>::cut_xy},
3334 
3335  {RefinementCase<dim - 1>::cut_y, // cut_z
3336  RefinementCase<dim - 1>::cut_x,
3337  RefinementCase<dim - 1>::no_refinement},
3338 
3339  {RefinementCase<dim - 1>::cut_y, // cut_xz
3340  RefinementCase<dim - 1>::cut_xy,
3341  RefinementCase<dim - 1>::cut_x},
3342 
3343  {RefinementCase<dim - 1>::cut_xy, // cut_yz
3344  RefinementCase<dim - 1>::cut_x,
3345  RefinementCase<dim - 1>::cut_y},
3346 
3347  {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3348  RefinementCase<dim - 1>::cut_xy,
3349  RefinementCase<dim - 1>::cut_xy},
3350  };
3351 
3352  const RefinementCase<dim - 1> ref_case =
3353  ref_cases[cell_refinement_case][face_no / 2];
3354 
3355  const RefinementCase<dim - 1> flip[4] = {
3356  RefinementCase<dim - 1>::no_refinement,
3357  RefinementCase<dim - 1>::cut_y,
3358  RefinementCase<dim - 1>::cut_x,
3359  RefinementCase<dim - 1>::cut_xy};
3360 
3361  // correct the ref_case for face_orientation
3362  // and face_rotation. for face_orientation,
3363  // 'true' is the default value whereas for
3364  // face_rotation, 'false' is standard. If
3365  // <tt>face_rotation==face_orientation</tt>,
3366  // then one of them is non-standard and we
3367  // have to swap cut_x and cut_y, otherwise no
3368  // change is necessary. face_flip has no
3369  // influence. however, in order to keep the
3370  // interface consistent with other functions,
3371  // we still include it as an argument to this
3372  // function
3373  return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3374 }
3375 
3376 
3377 
3378 template <int dim>
3379 inline RefinementCase<1>
3380 GeometryInfo<dim>::line_refinement_case(const RefinementCase<dim> &,
3381  const unsigned int)
3382 {
3383  Assert(false, ExcNotImplemented());
3384  return RefinementCase<1>::no_refinement;
3385 }
3386 
3387 template <>
3388 inline RefinementCase<1>
3389 GeometryInfo<1>::line_refinement_case(
3390  const RefinementCase<1> &cell_refinement_case,
3391  const unsigned int line_no)
3392 {
3393  (void)line_no;
3394  const unsigned int dim = 1;
3395  (void)dim;
3396  Assert(cell_refinement_case < RefinementCase<dim>::isotropic_refinement + 1,
3397  ExcIndexRange(cell_refinement_case,
3398  0,
3399  RefinementCase<dim>::isotropic_refinement + 1));
3400  Assert(line_no < GeometryInfo<dim>::lines_per_cell,
3401  ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
3402 
3403  return cell_refinement_case;
3404 }
3405 
3406 
3407 template <>
3408 inline RefinementCase<1>
3409 GeometryInfo<2>::line_refinement_case(
3410  const RefinementCase<2> &cell_refinement_case,
3411  const unsigned int line_no)
3412 {
3413  // Assertions are in face_refinement_case()
3414  return face_refinement_case(cell_refinement_case, line_no);
3415 }
3416 
3417 
3418 template <>
3419 inline RefinementCase<1>
3420 GeometryInfo<3>::line_refinement_case(
3421  const RefinementCase<3> &cell_refinement_case,
3422  const unsigned int line_no)
3423 {
3424  const unsigned int dim = 3;
3425  Assert(cell_refinement_case < RefinementCase<dim>::isotropic_refinement + 1,
3426  ExcIndexRange(cell_refinement_case,
3427  0,
3428  RefinementCase<dim>::isotropic_refinement + 1));
3429  Assert(line_no < GeometryInfo<dim>::lines_per_cell,
3430  ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
3431 
3432  // array indicating, which simple refine
3433  // case cuts a line in direction x, y or
3434  // z. For example, cut_y and everything
3435  // containing cut_y (cut_xy, cut_yz,
3436  // cut_xyz) cuts lines, which are in y
3437  // direction.
3438  const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3439  RefinementCase<dim>::cut_y,
3440  RefinementCase<dim>::cut_z};
3441 
3442  // order the direction of lines
3443  // 0->x, 1->y, 2->z
3444  const unsigned int direction[lines_per_cell] = {
3445  1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3446 
3447  return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3448  RefinementCase<1>::cut_x :
3449  RefinementCase<1>::no_refinement);
3450 }
3451 
3452 
3453 
3454 template <int dim>
3455 inline RefinementCase<dim>
3456 GeometryInfo<dim>::min_cell_refinement_case_for_face_refinement(
3457  const RefinementCase<dim - 1> &,
3458  const unsigned int,
3459  const bool,
3460  const bool,
3461  const bool)
3462 {
3463  Assert(false, ExcNotImplemented());
3464 
3465  return RefinementCase<dim>::no_refinement;
3466 }
3467 
3468 template <>
3469 inline RefinementCase<1>
3470 GeometryInfo<1>::min_cell_refinement_case_for_face_refinement(
3471  const RefinementCase<0> &,
3472  const unsigned int,
3473  const bool,
3474  const bool,
3475  const bool)
3476 {
3477  const unsigned int dim = 1;
3478  Assert(false, ExcImpossibleInDim(dim));
3479 
3480  return RefinementCase<dim>::no_refinement;
3481 }
3482 
3483 
3484 template <>
3485 inline RefinementCase<2>
3486 GeometryInfo<2>::min_cell_refinement_case_for_face_refinement(
3487  const RefinementCase<1> &face_refinement_case,
3488  const unsigned int face_no,
3489  const bool,
3490  const bool,
3491  const bool)
3492 {
3493  const unsigned int dim = 2;
3494  Assert(face_refinement_case <
3495  RefinementCase<dim - 1>::isotropic_refinement + 1,
3496  ExcIndexRange(face_refinement_case,
3497  0,
3498  RefinementCase<dim - 1>::isotropic_refinement + 1));
3499  Assert(face_no < GeometryInfo<dim>::faces_per_cell,
3500  ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
3501 
3502  if (face_refinement_case == RefinementCase<dim>::cut_x)
3503  return (face_no / 2) ? RefinementCase<dim>::cut_x :
3504  RefinementCase<dim>::cut_y;
3505  else
3506  return RefinementCase<dim>::no_refinement;
3507 }
3508 
3509 
3510 template <>
3511 inline RefinementCase<3>
3512 GeometryInfo<3>::min_cell_refinement_case_for_face_refinement(
3513  const RefinementCase<2> &face_refinement_case,
3514  const unsigned int face_no,
3515  const bool face_orientation,
3516  const bool /*face_flip*/,
3517  const bool face_rotation)
3518 {
3519  const unsigned int dim = 3;
3520  Assert(face_refinement_case <
3521  RefinementCase<dim - 1>::isotropic_refinement + 1,
3522  ExcIndexRange(face_refinement_case,
3523  0,
3524  RefinementCase<dim - 1>::isotropic_refinement + 1));
3525  Assert(face_no < GeometryInfo<dim>::faces_per_cell,
3526  ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
3527 
3528  const RefinementCase<2> flip[4] = {RefinementCase<2>::no_refinement,
3529  RefinementCase<2>::cut_y,
3530  RefinementCase<2>::cut_x,
3531  RefinementCase<2>::cut_xy};
3532 
3533  // correct the face_refinement_case for
3534  // face_orientation and face_rotation. for
3535  // face_orientation, 'true' is the default
3536  // value whereas for face_rotation, 'false'
3537  // is standard. If
3538  // <tt>face_rotation==face_orientation</tt>,
3539  // then one of them is non-standard and we
3540  // have to swap cut_x and cut_y, otherwise no
3541  // change is necessary. face_flip has no
3542  // influence. however, in order to keep the
3543  // interface consistent with other functions,
3544  // we still include it as an argument to this
3545  // function
3546  const RefinementCase<dim - 1> std_face_ref =
3547  (face_orientation == face_rotation) ? flip[face_refinement_case] :
3548  face_refinement_case;
3549 
3550  const RefinementCase<dim> face_to_cell[3][4] = {
3551  {RefinementCase<dim>::no_refinement, // faces 0 and 1
3552  RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3553  RefinementCase<dim>::cut_z,
3554  RefinementCase<dim>::cut_yz},
3555 
3556  {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3557  // "exchanged on faces 2 and 3")
3558  RefinementCase<dim>::cut_z,
3559  RefinementCase<dim>::cut_x,
3560  RefinementCase<dim>::cut_xz},
3561 
3562  {RefinementCase<dim>::no_refinement, // faces 4 and 5
3563  RefinementCase<dim>::cut_x,
3564  RefinementCase<dim>::cut_y,
3565  RefinementCase<dim>::cut_xy}};
3566 
3567  return face_to_cell[face_no / 2][std_face_ref];
3568 }
3569 
3570 
3571 
3572 template <int dim>
3573 inline RefinementCase<dim>
3574 GeometryInfo<dim>::min_cell_refinement_case_for_line_refinement(
3575  const unsigned int)
3576 {
3577  Assert(false, ExcNotImplemented());
3578 
3579  return RefinementCase<dim>::no_refinement;
3580 }
3581 
3582 template <>
3583 inline RefinementCase<1>
3584 GeometryInfo<1>::min_cell_refinement_case_for_line_refinement(
3585  const unsigned int line_no)
3586 {
3587  (void)line_no;
3588  Assert(line_no == 0, ExcIndexRange(line_no, 0, 1));
3589 
3590  return RefinementCase<1>::cut_x;
3591 }
3592 
3593 
3594 template <>
3595 inline RefinementCase<2>
3596 GeometryInfo<2>::min_cell_refinement_case_for_line_refinement(
3597  const unsigned int line_no)
3598 {
3599  const unsigned int dim = 2;
3600  (void)dim;
3601  Assert(line_no < GeometryInfo<dim>::lines_per_cell,
3602  ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
3603 
3604  return (line_no / 2) ? RefinementCase<2>::cut_x : RefinementCase<2>::cut_y;
3605 }
3606 
3607 
3608 template <>
3609 inline RefinementCase<3>
3610 GeometryInfo<3>::min_cell_refinement_case_for_line_refinement(
3611  const unsigned int line_no)
3612 {
3613  const unsigned int dim = 3;
3614  Assert(line_no < GeometryInfo<dim>::lines_per_cell,
3615  ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
3616 
3617  const RefinementCase<dim> ref_cases[6] = {
3618  RefinementCase<dim>::cut_y, // lines 0 and 1
3619  RefinementCase<dim>::cut_x, // lines 2 and 3
3620  RefinementCase<dim>::cut_y, // lines 4 and 5
3621  RefinementCase<dim>::cut_x, // lines 6 and 7
3622  RefinementCase<dim>::cut_z, // lines 8 and 9
3623  RefinementCase<dim>::cut_z}; // lines 10 and 11
3624 
3625  return ref_cases[line_no / 2];
3626 }
3627 
3628 
3629 
3630 template <>
3631 inline unsigned int
3632 GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
3633  const bool face_orientation,
3634  const bool face_flip,
3635  const bool face_rotation)
3636 {
3637  Assert(vertex < GeometryInfo<3>::vertices_per_face,
3638  ExcIndexRange(vertex, 0, GeometryInfo<3>::vertices_per_face));
3639 
3640  // set up a table to make sure that
3641  // we handle non-standard faces correctly
3642  //
3643  // so set up a table that for each vertex (of
3644  // a quad in standard position) describes
3645  // which vertex to take
3646  //
3647  // first index: four vertices 0...3
3648  //
3649  // second index: face_orientation; 0:
3650  // opposite normal, 1: standard
3651  //
3652  // third index: face_flip; 0: standard, 1:
3653  // face rotated by 180 degrees
3654  //
3655  // forth index: face_rotation: 0: standard,
3656  // 1: face rotated by 90 degrees
3657 
3658  const unsigned int vertex_translation[4][2][2][2] = {
3659  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3660  // face_rotation=false and true
3661  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3662  // face_rotation=false and true
3663  {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
3664  // face_rotation=false and true
3665  {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
3666  // face_rotation=false and true
3667 
3668  {{{2, 3}, // vertex 1 ...
3669  {1, 0}},
3670  {{1, 3}, {2, 0}}},
3671 
3672  {{{1, 0}, // vertex 2 ...
3673  {2, 3}},
3674  {{2, 0}, {1, 3}}},
3675 
3676  {{{3, 1}, // vertex 3 ...
3677  {0, 2}},
3678  {{3, 2}, {0, 1}}}};
3679 
3680  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3681 }
3682 
3683 
3684 
3685 template <int dim>
3686 inline unsigned int
3687 GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
3688  const bool,
3689  const bool,
3690  const bool)
3691 {
3692  Assert(dim > 1, ExcImpossibleInDim(dim));
3693  Assert(vertex < GeometryInfo<dim>::vertices_per_face,
3694  ExcIndexRange(vertex, 0, GeometryInfo<dim>::vertices_per_face));
3695  return vertex;
3696 }
3697 
3698 
3699 
3700 template <>
3701 inline unsigned int
3702 GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
3703  const bool face_orientation,
3704  const bool face_flip,
3705  const bool face_rotation)
3706 {
3707  Assert(line < GeometryInfo<3>::lines_per_face,
3708  ExcIndexRange(line, 0, GeometryInfo<3>::lines_per_face));
3709 
3710 
3711  // make sure we handle
3712  // non-standard faces correctly
3713  //
3714  // so set up a table that for each line (of a
3715  // quad) describes which line to take
3716  //
3717  // first index: four lines 0...3
3718  //
3719  // second index: face_orientation; 0:
3720  // opposite normal, 1: standard
3721  //
3722  // third index: face_flip; 0: standard, 1:
3723  // face rotated by 180 degrees
3724  //
3725  // forth index: face_rotation: 0: standard,
3726  // 1: face rotated by 90 degrees
3727 
3728  const unsigned int line_translation[4][2][2][2] = {
3729  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
3730  // face_rotation=false and true
3731  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
3732  // face_rotation=false and true
3733  {{0, 3}, // line 0, face_orientation=true, face_flip=false,
3734  // face_rotation=false and true
3735  {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
3736  // face_rotation=false and true
3737 
3738  {{{3, 1}, // line 1 ...
3739  {2, 0}},
3740  {{1, 2}, {0, 3}}},
3741 
3742  {{{0, 3}, // line 2 ...
3743  {1, 2}},
3744  {{2, 0}, {3, 1}}},
3745 
3746  {{{1, 2}, // line 3 ...
3747  {0, 3}},
3748  {{3, 1}, {2, 0}}}};
3749 
3750  return line_translation[line][face_orientation][face_flip][face_rotation];
3751 }
3752 
3753 
3754 
3755 template <int dim>
3756 inline unsigned int
3757 GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
3758  const bool,
3759  const bool,
3760  const bool)
3761 {
3762  Assert(false, ExcNotImplemented());
3763  return line;
3764 }
3765 
3766 
3767 
3768 template <>
3769 inline unsigned int
3770 GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
3771  const bool face_orientation,
3772  const bool face_flip,
3773  const bool face_rotation)
3774 {
3775  Assert(line < GeometryInfo<3>::lines_per_face,
3776  ExcIndexRange(line, 0, GeometryInfo<3>::lines_per_face));
3777 
3778 
3779  // make sure we handle
3780  // non-standard faces correctly
3781  //
3782  // so set up a table that for each line (of a
3783  // quad) describes which line to take
3784  //
3785  // first index: four lines 0...3
3786  //
3787  // second index: face_orientation; 0:
3788  // opposite normal, 1: standard
3789  //
3790  // third index: face_flip; 0: standard, 1:
3791  // face rotated by 180 degrees
3792  //
3793  // forth index: face_rotation: 0: standard,
3794  // 1: face rotated by 90 degrees
3795 
3796  const unsigned int line_translation[4][2][2][2] = {
3797  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
3798  // face_rotation=false and true
3799  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
3800  // face_rotation=false and true
3801  {{0, 2}, // line 0, face_orientation=true, face_flip=false,
3802  // face_rotation=false and true
3803  {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
3804  // face_rotation=false and true
3805 
3806  {{{3, 1}, // line 1 ...
3807  {2, 0}},
3808  {{1, 3}, {0, 2}}},
3809 
3810  {{{0, 3}, // line 2 ...
3811  {1, 2}},
3812  {{2, 1}, {3, 0}}},
3813 
3814  {{{1, 2}, // line 3 ...
3815  {0, 3}},
3816  {{3, 0}, {2, 1}}}};
3817 
3818  return line_translation[line][face_orientation][face_flip][face_rotation];
3819 }
3820 
3821 
3822 
3823 template <int dim>
3824 inline unsigned int
3825 GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
3826  const bool,
3827  const bool,
3828  const bool)
3829 {
3830  Assert(false, ExcNotImplemented());
3831  return line;
3832 }
3833 
3834 
3835 
3836 template <>
3837 inline unsigned int
3838 GeometryInfo<1>::child_cell_on_face(const RefinementCase<1> &,
3839  const unsigned int face,
3840  const unsigned int subface,
3841  const bool,
3842  const bool,
3843  const bool,
3844  const RefinementCase<0> &)
3845 {
3846  (void)subface;
3847  Assert(face < faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
3848  Assert(subface < max_children_per_face,
3849  ExcIndexRange(subface, 0, max_children_per_face));
3850 
3851  return face;
3852 }
3853 
3854 
3855 
3856 template <>
3857 inline unsigned int
3858 GeometryInfo<2>::child_cell_on_face(const RefinementCase<2> &ref_case,
3859  const unsigned int face,
3860  const unsigned int subface,
3861  const bool /*face_orientation*/,
3862  const bool face_flip,
3863  const bool /*face_rotation*/,
3864  const RefinementCase<1> &)
3865 {
3866  Assert(face < faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
3867  Assert(subface < max_children_per_face,
3868  ExcIndexRange(subface, 0, max_children_per_face));
3869 
3870  // always return the child adjacent to the specified
3871  // subface. if the face of a cell is not refined, don't
3872  // throw an assertion but deliver the child adjacent to
3873  // the face nevertheless, i.e. deliver the child of
3874  // this cell adjacent to the subface of a possibly
3875  // refined neighbor. this simplifies setting neighbor
3876  // information in execute_refinement.
3877  const unsigned int
3878  subcells[2][RefinementCase<2>::isotropic_refinement][faces_per_cell]
3879  [max_children_per_face] = {
3880  {
3881  // Normal orientation (face_flip = false)
3882  {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
3883  {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
3884  {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
3885  },
3886  {
3887  // Flipped orientation (face_flip = true)
3888  {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
3889  {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
3890  {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
3891  }};
3892 
3893  return subcells[face_flip][ref_case - 1][face][subface];
3894 }
3895 
3896 
3897 
3898 template <>
3899 inline unsigned int
3900 GeometryInfo<3>::child_cell_on_face(const RefinementCase<3> &ref_case,
3901  const unsigned int face,
3902  const unsigned int subface,
3903  const bool face_orientation,
3904  const bool face_flip,
3905  const bool face_rotation,
3906  const RefinementCase<2> &face_ref_case)
3907 {
3908  const unsigned int dim = 3;
3909 
3910  Assert(ref_case > RefinementCase<dim - 1>::no_refinement,
3911  ExcMessage("Cell has no children."));
3912  Assert(face < faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
3913  Assert(subface < GeometryInfo<dim - 1>::n_children(face_ref_case) ||
3914  (subface == 0 &&
3915  face_ref_case == RefinementCase<dim - 1>::no_refinement),
3916  ExcIndexRange(subface, 0, GeometryInfo<2>::n_children(face_ref_case)));
3917 
3918  // invalid number used for invalid cases,
3919  // e.g. when the children are more refined at
3920  // a given face than the face itself
3921  const unsigned int e = numbers::invalid_unsigned_int;
3922 
3923  // the whole process of finding a child cell
3924  // at a given subface considering the
3925  // possibly anisotropic refinement cases of
3926  // the cell and the face as well as
3927  // orientation, flip and rotation of the face
3928  // is quite complicated. thus, we break it
3929  // down into several steps.
3930 
3931  // first step: convert the given face refine
3932  // case to a face refine case concerning the
3933  // face in standard orientation (, flip and
3934  // rotation). This only affects cut_x and
3935  // cut_y
3936  const RefinementCase<dim - 1> flip[4] = {
3937  RefinementCase<dim - 1>::no_refinement,
3938  RefinementCase<dim - 1>::cut_y,
3939  RefinementCase<dim - 1>::cut_x,
3940  RefinementCase<dim - 1>::cut_xy};
3941  // for face_orientation, 'true' is the
3942  // default value whereas for face_rotation,
3943  // 'false' is standard. If
3944  // <tt>face_rotation==face_orientation</tt>,
3945  // then one of them is non-standard and we
3946  // have to swap cut_x and cut_y, otherwise no
3947  // change is necessary.
3948  const RefinementCase<dim - 1> std_face_ref =
3949  (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
3950 
3951  // second step: convert the given subface
3952  // index to the one for a standard face
3953  // respecting face_orientation, face_flip and
3954  // face_rotation
3955 
3956  // first index: face_ref_case
3957  // second index: face_orientation
3958  // third index: face_flip
3959  // forth index: face_rotation
3960  // fifth index: subface index
3961  const unsigned int subface_exchange[4][2][2][2][4] = {
3962  // no_refinement (subface 0 stays 0,
3963  // all others are invalid)
3964  {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
3965  {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
3966  // cut_x (here, if the face is only
3967  // rotated OR only falsely oriented,
3968  // then subface 0 of the non-standard
3969  // face does NOT correspond to one of
3970  // the subfaces of a standard
3971  // face. Thus we indicate the subface
3972  // which is located at the lower left
3973  // corner (the origin of the face's
3974  // local coordinate system) with
3975  // '0'. The rest of this issue is
3976  // taken care of using the above
3977  // conversion to a 'standard face
3978  // refine case')
3979  {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
3980  {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
3981  // cut_y (the same applies as for
3982  // cut_x)
3983  {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
3984  {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
3985  // cut_xyz: this information is
3986  // identical to the information
3987  // returned by
3988  // GeometryInfo<3>::real_to_standard_face_vertex()
3989  {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
3990  // face_rotation=false, subfaces 0,1,2,3
3991  {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
3992  // face_rotation=true, subfaces 0,1,2,3
3993  {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
3994  // face_rotation=false, subfaces 0,1,2,3
3995  {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
3996  // face_rotation=true, subfaces 0,1,2,3
3997  {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
3998  // face_rotation=false, subfaces 0,1,2,3
3999  {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4000  // face_rotation=true, subfaces 0,1,2,3
4001  {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4002  // face_rotation=false, subfaces 0,1,2,3
4003  {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4004  // face_rotation=true, subfaces 0,1,2,3
4005 
4006  const unsigned int std_subface =
4007  subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4008  [subface];
4009  Assert(std_subface != e, ExcInternalError());
4010 
4011  // third step: these are the children, which
4012  // can be found at the given subfaces of an
4013  // isotropically refined (standard) face
4014  //
4015  // first index: (refinement_case-1)
4016  // second index: face_index
4017  // third index: subface_index (isotropic refinement)
4018  const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4019  [max_children_per_face] = {
4020  // cut_x
4021  {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4022  {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4023  {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4024  {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4025  {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4026  {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4027  // cut_y
4028  {{0, 1, 0, 1},
4029  {0, 1, 0, 1},
4030  {0, 0, 0, 0},
4031  {1, 1, 1, 1},
4032  {0, 0, 1, 1},
4033  {0, 0, 1, 1}},
4034  // cut_xy
4035  {{0, 2, 0, 2},
4036  {1, 3, 1, 3},
4037  {0, 0, 1, 1},
4038  {2, 2, 3, 3},
4039  {0, 1, 2, 3},
4040  {0, 1, 2, 3}},
4041  // cut_z
4042  {{0, 0, 1, 1},
4043  {0, 0, 1, 1},
4044  {0, 1, 0, 1},
4045  {0, 1, 0, 1},
4046  {0, 0, 0, 0},
4047  {1, 1, 1, 1}},
4048  // cut_xz
4049  {{0, 0, 1, 1},
4050  {2, 2, 3, 3},
4051  {0, 1, 2, 3},
4052  {0, 1, 2, 3},
4053  {0, 2, 0, 2},
4054  {1, 3, 1, 3}},
4055  // cut_yz
4056  {{0, 1, 2, 3},
4057  {0, 1, 2, 3},
4058  {0, 2, 0, 2},
4059  {1, 3, 1, 3},
4060  {0, 0, 1, 1},
4061  {2, 2, 3, 3}},
4062  // cut_xyz
4063  {{0, 2, 4, 6},
4064  {1, 3, 5, 7},
4065  {0, 4, 1, 5},
4066  {2, 6, 3, 7},
4067  {0, 1, 2, 3},
4068  {4, 5, 6, 7}}};
4069 
4070  // forth step: check, whether the given face
4071  // refine case is valid for the given cell
4072  // refine case. this is the case, if the
4073  // given face refine case is at least as
4074  // refined as the face is for the given cell
4075  // refine case
4076 
4077  // note, that we are considering standard
4078  // face refinement cases here and thus must
4079  // not pass the given orientation, flip and
4080  // rotation flags
4081  if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4082  face_refinement_case(ref_case, face))
4083  {
4084  // all is fine. for anisotropic face
4085  // refine cases, select one of the
4086  // isotropic subfaces which neighbors the
4087  // same child
4088 
4089  // first index: (standard) face refine case
4090  // second index: subface index
4091  const unsigned int equivalent_iso_subface[4][4] = {
4092  {0, e, e, e}, // no_refinement
4093  {0, 3, e, e}, // cut_x
4094  {0, 3, e, e}, // cut_y
4095  {0, 1, 2, 3}}; // cut_xy
4096 
4097  const unsigned int equ_std_subface =
4098  equivalent_iso_subface[std_face_ref][std_subface];
4099  Assert(equ_std_subface != e, ExcInternalError());
4100 
4101  return iso_children[ref_case - 1][face][equ_std_subface];
4102  }
4103  else
4104  {
4105  // the face_ref_case was too coarse,
4106  // throw an error
4107  Assert(false,
4108  ExcMessage("The face RefineCase is too coarse "
4109  "for the given cell RefineCase."));
4110  }
4111  // we only get here in case of an error
4112  return e;
4113 }
4114 
4115 
4116 
4117 template <>
4118 inline unsigned int
4119 GeometryInfo<4>::child_cell_on_face(const RefinementCase<4> &,
4120  const unsigned int,
4121  const unsigned int,
4122  const bool,
4123  const bool,
4124  const bool,
4125  const RefinementCase<3> &)
4126 {
4127  Assert(false, ExcNotImplemented());
4128  return numbers::invalid_unsigned_int;
4129 }
4130 
4131 
4132 
4133 template <>
4134 inline unsigned int
4135 GeometryInfo<1>::face_to_cell_lines(const unsigned int face,
4136  const unsigned int line,
4137  const bool,
4138  const bool,
4139  const bool)
4140 {
4141  (void)face;
4142  (void)line;
4143  Assert(face + 1 < faces_per_cell + 1, ExcIndexRange(face, 0, faces_per_cell));
4144  Assert(line + 1 < lines_per_face + 1, ExcIndexRange(line, 0, lines_per_face));
4145 
4146  // There is only a single line, so
4147  // it must be this.
4148  return 0;
4149 }
4150 
4151 
4152 
4153 template <>
4154 inline unsigned int
4155 GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4156  const unsigned int line,
4157  const bool,
4158  const bool,
4159  const bool)
4160 {
4161  (void)line;
4162  Assert(face < faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
4163  Assert(line < lines_per_face, ExcIndexRange(line, 0, lines_per_face));
4164 
4165  // The face is a line itself.
4166  return face;
4167 }
4168 
4169 
4170 
4171 template <>
4172 inline unsigned int
4173 GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4174  const unsigned int line,
4175  const bool face_orientation,
4176  const bool face_flip,
4177  const bool face_rotation)
4178 {
4179  Assert(face < faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
4180  Assert(line < lines_per_face, ExcIndexRange(line, 0, lines_per_face));
4181 
4182  const unsigned lines[faces_per_cell][lines_per_face] = {
4183  {8, 10, 0, 4}, // left face
4184  {9, 11, 1, 5}, // right face
4185  {2, 6, 8, 9}, // front face
4186  {3, 7, 10, 11}, // back face
4187  {0, 1, 2, 3}, // bottom face
4188  {4, 5, 6, 7}}; // top face
4189  return lines[face][real_to_standard_face_line(
4190  line, face_orientation, face_flip, face_rotation)];
4191 }
4192 
4193 
4194 
4195 template <int dim>
4196 inline unsigned int
4197 GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4198  const unsigned int,
4199  const bool,
4200  const bool,
4201  const bool)
4202 {
4203  Assert(false, ExcNotImplemented());
4204  return numbers::invalid_unsigned_int;
4205 }
4206 
4207 
4208 
4209 template <int dim>
4210 inline unsigned int
4211 GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4212  const unsigned int vertex,
4213  const bool face_orientation,
4214  const bool face_flip,
4215  const bool face_rotation)
4216 {
4217  return child_cell_on_face(RefinementCase<dim>::isotropic_refinement,
4218  face,
4219  vertex,
4220  face_orientation,
4221  face_flip,
4222  face_rotation);
4223 }
4224 
4225 
4226 
4227 template <int dim>
4228 inline Point<dim>
4229 GeometryInfo<dim>::project_to_unit_cell(const Point<dim> &q)
4230 {
4231  Point<dim> p = q;
4232  for (unsigned int i = 0; i < dim; i++)
4233  if (p[i] < 0.)
4234  p[i] = 0.;
4235  else if (p[i] > 1.)
4236  p[i] = 1.;
4237 
4238  return p;
4239 }
4240 
4241 
4242 
4243 template <int dim>
4244 inline double
4245 GeometryInfo<dim>::distance_to_unit_cell(const Point<dim> &p)
4246 {
4247  double result = 0.0;
4248 
4249  for (unsigned int i = 0; i < dim; i++)
4250  if ((-p[i]) > result)
4251  result = -p[i];
4252  else if ((p[i] - 1.) > result)
4253  result = (p[i] - 1.);
4254 
4255  return result;
4256 }
4257 
4258 
4259 
4260 template <int dim>
4261 inline double
4262 GeometryInfo<dim>::d_linear_shape_function(const Point<dim> & xi,
4263  const unsigned int i)
4264 {
4265  Assert(i < GeometryInfo<dim>::vertices_per_cell,
4266  ExcIndexRange(i, 0, GeometryInfo<dim>::vertices_per_cell));
4267 
4268  switch (dim)
4269  {
4270  case 1:
4271  {
4272  const double x = xi[0];
4273  switch (i)
4274  {
4275  case 0:
4276  return 1 - x;
4277  case 1:
4278  return x;
4279  }
4280  break;
4281  }
4282 
4283  case 2:
4284  {
4285  const double x = xi[0];
4286  const double y = xi[1];
4287  switch (i)
4288  {
4289  case 0:
4290  return (1 - x) * (1 - y);
4291  case 1:
4292  return x * (1 - y);
4293  case 2:
4294  return (1 - x) * y;
4295  case 3:
4296  return x * y;
4297  }
4298  break;
4299  }
4300 
4301  case 3:
4302  {
4303  const double x = xi[0];
4304  const double y = xi[1];
4305  const double z = xi[2];
4306  switch (i)
4307  {
4308  case 0:
4309  return (1 - x) * (1 - y) * (1 - z);
4310  case 1:
4311  return x * (1 - y) * (1 - z);
4312  case 2:
4313  return (1 - x) * y * (1 - z);
4314  case 3:
4315  return x * y * (1 - z);
4316  case 4:
4317  return (1 - x) * (1 - y) * z;
4318  case 5:
4319  return x * (1 - y) * z;
4320  case 6:
4321  return (1 - x) * y * z;
4322  case 7:
4323  return x * y * z;
4324  }
4325  break;
4326  }
4327 
4328  default:
4329  Assert(false, ExcNotImplemented());
4330  }
4331  return -1e9;
4332 }
4333 
4334 
4335 
4336 template <>
4337 Tensor<1, 1> inline GeometryInfo<1>::d_linear_shape_function_gradient(
4338  const Point<1> &,
4339  const unsigned int i)
4340 {
4341  Assert(i < GeometryInfo<1>::vertices_per_cell,
4342  ExcIndexRange(i, 0, GeometryInfo<1>::vertices_per_cell));
4343 
4344  switch (i)
4345  {
4346  case 0:
4347  return Point<1>(-1.);
4348  case 1:
4349  return Point<1>(1.);
4350  }
4351 
4352  return Point<1>(-1e9);
4353 }
4354 
4355 
4356 
4357 template <>
4358 Tensor<1, 2> inline GeometryInfo<2>::d_linear_shape_function_gradient(
4359  const Point<2> & xi,
4360  const unsigned int i)
4361 {
4362  Assert(i < GeometryInfo<2>::vertices_per_cell,
4363  ExcIndexRange(i, 0, GeometryInfo<2>::vertices_per_cell));
4364 
4365  const double x = xi[0];
4366  const double y = xi[1];
4367  switch (i)
4368  {
4369  case 0:
4370  return Point<2>(-(1 - y), -(1 - x));
4371  case 1:
4372  return Point<2>(1 - y, -x);
4373  case 2:
4374  return Point<2>(-y, 1 - x);
4375  case 3:
4376  return Point<2>(y, x);
4377  }
4378  return Point<2>(-1e9, -1e9);
4379 }
4380 
4381 
4382 
4383 template <>
4384 Tensor<1, 3> inline GeometryInfo<3>::d_linear_shape_function_gradient(
4385  const Point<3> & xi,
4386  const unsigned int i)
4387 {
4388  Assert(i < GeometryInfo<3>::vertices_per_cell,
4389  ExcIndexRange(i, 0, GeometryInfo<3>::vertices_per_cell));
4390 
4391  const double x = xi[0];
4392  const double y = xi[1];
4393  const double z = xi[2];
4394  switch (i)
4395  {
4396  case 0:
4397  return Point<3>(-(1 - y) * (1 - z),
4398  -(1 - x) * (1 - z),
4399  -(1 - x) * (1 - y));
4400  case 1:
4401  return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4402  case 2:
4403  return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4404  case 3:
4405  return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4406  case 4:
4407  return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4408  case 5:
4409  return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4410  case 6:
4411  return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4412  case 7:
4413  return Point<3>(y * z, x * z, x * y);
4414  }
4415 
4416  return Point<3>(-1e9, -1e9, -1e9);
4417 }
4418 
4419 
4420 
4421 template <int dim>
4422 inline Tensor<1, dim>
4423 GeometryInfo<dim>::d_linear_shape_function_gradient(const Point<dim> &,
4424  const unsigned int)
4425 {
4426  Assert(false, ExcNotImplemented());
4427  return Tensor<1, dim>();
4428 }
4429 
4430 
4431 unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4432 {
4433  return 0;
4434 }
4435 
4436 
4437 namespace internal
4438 {
4439  namespace GeometryInfoHelper
4440  {
4441  // wedge product of a single
4442  // vector in 2d: we just have to
4443  // rotate it by 90 degrees to the
4444  // right
4445  inline Tensor<1, 2>
4446  wedge_product(const Tensor<1, 2> (&derivative)[1])
4447  {
4448  Tensor<1, 2> result;
4449  result[0] = derivative[0][1];
4450  result[1] = -derivative[0][0];
4451 
4452  return result;
4453  }
4454 
4455 
4456  // wedge product of 2 vectors in
4457  // 3d is the cross product
4458  inline Tensor<1, 3>
4459  wedge_product(const Tensor<1, 3> (&derivative)[2])
4460  {
4461  return cross_product_3d(derivative[0], derivative[1]);
4462  }
4463 
4464 
4465  // wedge product of dim vectors
4466  // in dim-d: that's the
4467  // determinant of the matrix
4468  template <int dim>
4469  inline Tensor<0, dim>
4470  wedge_product(const Tensor<1, dim> (&derivative)[dim])
4471  {
4472  Tensor<2, dim> jacobian;
4473  for (unsigned int i = 0; i < dim; ++i)
4474  jacobian[i] = derivative[i];
4475 
4476  return determinant(jacobian);
4477  }
4478  } // namespace GeometryInfoHelper
4479 } // namespace internal
4480 
4481 
4482 template <int dim>
4483 template <int spacedim>
4484 inline void
4485 GeometryInfo<dim>::alternating_form_at_vertices
4486 # ifndef DEAL_II_CONSTEXPR_BUG
4487  (const Point<spacedim> (&vertices)[vertices_per_cell],
4488  Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell])
4489 # else
4490  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms)
4491 # endif
4492 {
4493  // for each of the vertices,
4494  // compute the alternating form
4495  // of the mapped unit
4496  // vectors. consider for
4497  // example the case of a quad
4498  // in spacedim==3: to do so, we
4499  // need to see how the
4500  // infinitesimal vectors
4501  // (d\xi_1,0) and (0,d\xi_2)
4502  // are transformed into
4503  // spacedim-dimensional space
4504  // and then form their cross
4505  // product (i.e. the wedge product
4506  // of two vectors). to this end, note
4507  // that
4508  // \vec x = sum_i \vec v_i phi_i(\vec xi)
4509  // so the transformed vectors are
4510  // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
4511  // and
4512  // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
4513  // which boils down to the columns
4514  // of the 3x2 matrix \grad_\xi x(\xi)
4515  //
4516  // a similar reasoning would
4517  // hold for all dim,spacedim
4518  // pairs -- we only have to
4519  // compute the wedge product of
4520  // the columns of the
4521  // derivatives
4522  for (unsigned int i = 0; i < vertices_per_cell; ++i)
4523  {
4524  Tensor<1, spacedim> derivatives[dim];
4525 
4526  for (unsigned int j = 0; j < vertices_per_cell; ++j)
4527  {
4528  const Tensor<1, dim> grad_phi_j =
4529  d_linear_shape_function_gradient(unit_cell_vertex(i), j);
4530  for (unsigned int l = 0; l < dim; ++l)
4531  derivatives[l] += vertices[j] * grad_phi_j[l];
4532  }
4533 
4534  forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
4535  }
4536 }
4537 
4538 #endif // DOXYGEN
4539 DEAL_II_NAMESPACE_CLOSE
4540 
4541 #endif
GeometryInfo::child_to_cell_coordinates
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)
GeometryInfo::real_to_standard_face_vertex
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)
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:173
GeometryInfo::standard_to_real_face_line
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)
determinant
Number determinant(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:2676
GeometryInfo::child_cell_from_point
static unsigned int child_cell_from_point(const Point< dim > &p)
GeometryInfo::face_to_cell_vertices
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)
GeometryInfo::line_refinement_case
static RefinementCase< 1 > line_refinement_case(const RefinementCase< dim > &cell_refinement_case, const unsigned int line_no)
GeometryInfo::face_to_cell_lines
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)
GeometryInfo::max_children_per_face
static constexpr unsigned int max_children_per_face
Definition: geometry_info.h:1719
RefinementCase::operator|
RefinementCase operator|(const RefinementCase &r) const
GeometryInfo::is_inside_unit_cell
static bool is_inside_unit_cell(const Point< dim > &p)
GeometryPrimitive::GeometryPrimitive
GeometryPrimitive(const Object object)
RefinementCase::serialize
void serialize(Archive &ar, const unsigned int version)
GeometryInfo::line_to_cell_vertices
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
GeometryInfo::real_to_standard_face_line
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)
GeometryInfo::cell_to_child_coordinates
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)
GeometryInfo::unit_cell_vertex
static Point< dim > unit_cell_vertex(const unsigned int vertex)
GeometryInfo::d_linear_shape_function_gradient
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
StandardExceptions::ExcIndexRange
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
GeometryInfo< 0 >::ucd_to_deal
static const std::array< unsigned int, vertices_per_cell > ucd_to_deal
Definition: geometry_info.h:1151
GeometryInfo::vertices_per_cell
static constexpr unsigned int vertices_per_cell
Definition: geometry_info.h:1725
operator&
UpdateFlags operator&(const UpdateFlags f1, const UpdateFlags f2)
Definition: fe_update_flags.h:339
RefinementCase::operator~
RefinementCase operator~() const
GeometryInfo::distance_to_unit_cell
static double distance_to_unit_cell(const Point< dim > &p)
operator|
UpdateFlags operator|(const UpdateFlags f1, const UpdateFlags f2)
Definition: fe_update_flags.h:308
StandardExceptions::ExcMessage
static::ExceptionBase & ExcMessage(std::string arg1)
StandardExceptions::ExcImpossibleInDim
static::ExceptionBase & ExcImpossibleInDim(int arg1)
DeclException1
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:408
GeometryInfo::lines_per_cell
static constexpr unsigned int lines_per_cell
Definition: geometry_info.h:1754
GeometryInfo::child_cell_on_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)
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1227
GeometryInfo::face_refinement_case
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)
GeometryInfo::faces_per_cell
static constexpr unsigned int faces_per_cell
Definition: geometry_info.h:1710
GeometryInfo::ExcInvalidCoordinate
static::ExceptionBase & ExcInvalidCoordinate(double arg1)
RefinementCase::memory_consumption
static std::size_t memory_consumption()
GeometryInfo::standard_to_real_face_vertex
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)
GeometryInfo::min_cell_refinement_case_for_line_refinement
static RefinementCase< dim > min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
RefinementCase::operator&
RefinementCase operator&(const RefinementCase &r) const
GeometryInfo::n_children
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
GeometryInfo::project_to_unit_cell
static Point< dim > project_to_unit_cell(const Point< dim > &p)
GeometryInfo< 0 >::dx_to_deal
static const std::array< unsigned int, vertices_per_cell > dx_to_deal
Definition: geometry_info.h:1166
GeometryInfo::d_linear_shape_function
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
GeometryInfo::lines_per_face
static constexpr unsigned int lines_per_face
Definition: geometry_info.h:1736
GeometryInfo::n_subfaces
static unsigned int n_subfaces(const internal::SubfaceCase< dim > &subface_case)
StandardExceptions::ExcNotImplemented
static::ExceptionBase & ExcNotImplemented()
RefinementCase::value
std::uint8_t value
Definition: geometry_info.h:708
DeclException3
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:432
GeometryInfo::subface_ratio
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
GeometryInfo::min_cell_refinement_case_for_face_refinement
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)
RefinementCase::cut_axis
static RefinementCase cut_axis(const unsigned int i)
StandardExceptions::ExcInternalError
static::ExceptionBase & ExcInternalError()
Generated by   doxygen 1.8.11