Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
geometry_info.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_geometry_info_h
17 #define dealii_geometry_info_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/point.h>
24 
25 #include <cstdint>
26 
27 
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 #ifndef DOXYGEN
32 namespace internal
33 {
34  namespace GeometryInfoHelper
35  {
36  // A struct that holds the values for all the arrays we want to initialize
37  // in GeometryInfo
38  template <int dim>
39  struct Initializers;
40 
41  template <>
42  struct Initializers<1>
43  {
44  static constexpr std::array<unsigned int, 2> 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 
368  isotropic_refinement = 0xFF
369  };
370 };
371 
372 
373 
384 template <>
385 struct RefinementPossibilities<1>
386 {
421  enum Possibilities
422  {
426  no_refinement = 0,
430  cut_x = 1,
434  isotropic_refinement = cut_x
435  };
436 };
437 
438 
439 
451 template <>
452 struct RefinementPossibilities<2>
453 {
488  enum Possibilities
489  {
493  no_refinement = 0,
497  cut_x = 1,
501  cut_y = 2,
505  cut_xy = cut_x | cut_y,
506 
510  isotropic_refinement = cut_xy
511  };
512 };
513 
514 
515 
527 template <>
528 struct RefinementPossibilities<3>
529 {
564  enum Possibilities
565  {
569  no_refinement = 0,
573  cut_x = 1,
577  cut_y = 2,
581  cut_xy = cut_x | cut_y,
585  cut_z = 4,
589  cut_xz = cut_x | cut_z,
593  cut_yz = cut_y | cut_z,
597  cut_xyz = cut_x | cut_y | cut_z,
598 
602  isotropic_refinement = cut_xyz
603  };
604 };
605 
606 
607 
619 template <int dim>
620 class RefinementCase : public RefinementPossibilities<dim>
621 {
622 public:
626  RefinementCase();
627 
632  RefinementCase(
633  const typename RefinementPossibilities<dim>::Possibilities refinement_case);
634 
640  explicit RefinementCase(const std::uint8_t refinement_case);
641 
653  operator std::uint8_t() const;
654 
659  RefinementCase
660  operator|(const RefinementCase &r) const;
661 
666  RefinementCase operator&(const RefinementCase &r) const;
667 
675  RefinementCase
676  operator~() const;
677 
678 
684  static RefinementCase
685  cut_axis(const unsigned int i);
686 
690  static std::size_t
691  memory_consumption();
692 
697  template <class Archive>
698  void
699  serialize(Archive &ar, const unsigned int version);
700 
704  DeclException1(
705  ExcInvalidRefinementCase,
706  int,
707  << "The refinement flags given (" << arg1
708  << ") contain set bits that do not "
709  << "make sense for the space dimension of the object to which they are applied.");
710 
711 private:
716  std::uint8_t value : (dim > 0 ? dim : 1);
717 };
718 
719 
720 namespace internal
721 {
742  template <int dim>
743  struct SubfacePossibilities
744  {
748  enum Possibilities
749  {
753  case_none = 0,
754 
758  case_isotropic = static_cast<std::uint8_t>(-1)
759  };
760  };
761 
762 
772  template <>
773  struct SubfacePossibilities<0>
774  {
780  enum Possibilities
781  {
785  case_none = 0,
786 
790  case_isotropic = case_none
791  };
792  };
793 
794 
795 
806  template <>
807  struct SubfacePossibilities<1>
808  {
814  enum Possibilities
815  {
819  case_none = 0,
820 
824  case_isotropic = case_none
825  };
826  };
827 
828 
829 
841  template <>
842  struct SubfacePossibilities<2>
843  {
850  enum Possibilities
851  {
855  case_none = 0,
859  case_x = 1,
863  case_isotropic = case_x
864  };
865  };
866 
867 
868 
961  template <>
962  struct SubfacePossibilities<3>
963  {
970  enum Possibilities
971  {
972  case_none = 0,
973  case_x = 1,
974  case_x1y = 2,
975  case_x2y = 3,
976  case_x1y2y = 4,
977  case_y = 5,
978  case_y1x = 6,
979  case_y2x = 7,
980  case_y1x2x = 8,
981  case_xy = 9,
982 
983  case_isotropic = case_xy
984  };
985  };
986 
987 
988 
996  template <int dim>
997  class SubfaceCase : public SubfacePossibilities<dim>
998  {
999  public:
1005  SubfaceCase(const typename SubfacePossibilities<dim>::Possibilities
1006  subface_possibility);
1007 
1019  operator std::uint8_t() const;
1020 
1024  static constexpr std::size_t
1025  memory_consumption();
1026 
1030  DeclException1(
1031  ExcInvalidSubfaceCase,
1032  int,
1033  << "The subface case given (" << arg1 << ") does not make sense "
1034  << "for the space dimension of the object to which they are applied.");
1035 
1036  private:
1041  std::uint8_t value : (dim == 3 ? 4 : 1);
1042  };
1043 
1044 } // namespace internal
1045 
1046 
1047 
1048 template <int dim>
1049 struct GeometryInfo;
1050 
1051 
1052 
1070 template <>
1071 struct GeometryInfo<0>
1072 {
1080  static constexpr unsigned int max_children_per_cell = 1;
1081 
1085  static constexpr unsigned int faces_per_cell = 0;
1086 
1094  static constexpr unsigned int max_children_per_face = 0;
1095 
1101  static unsigned int
1102  n_children(const RefinementCase<0> &refinement_case);
1103 
1107  static constexpr unsigned int vertices_per_cell = 1;
1108 
1115  static constexpr unsigned int vertices_per_face = 0;
1116 
1120  static constexpr unsigned int lines_per_face = 0;
1121 
1125  static constexpr unsigned int quads_per_face = 0;
1126 
1130  static constexpr unsigned int lines_per_cell = 0;
1131 
1135  static constexpr unsigned int quads_per_cell = 0;
1136 
1140  static constexpr unsigned int hexes_per_cell = 0;
1141 
1159  static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1160 
1174  static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1175 };
1176 
1177 
1178 
1703 template <int dim>
1704 struct GeometryInfo
1705 {
1713  static constexpr unsigned int max_children_per_cell = 1 << dim;
1714 
1718  static constexpr unsigned int faces_per_cell = 2 * dim;
1719 
1727  static constexpr unsigned int max_children_per_face =
1728  GeometryInfo<dim - 1>::max_children_per_cell;
1729 
1733  static constexpr unsigned int vertices_per_cell = 1 << dim;
1734 
1738  static constexpr unsigned int vertices_per_face =
1739  GeometryInfo<dim - 1>::vertices_per_cell;
1740 
1744  static constexpr unsigned int lines_per_face =
1745  GeometryInfo<dim - 1>::lines_per_cell;
1746 
1750  static constexpr unsigned int quads_per_face =
1751  GeometryInfo<dim - 1>::quads_per_cell;
1752 
1762  static constexpr unsigned int lines_per_cell =
1763  (2 * GeometryInfo<dim - 1>::lines_per_cell +
1764  GeometryInfo<dim - 1>::vertices_per_cell);
1765 
1773  static constexpr unsigned int quads_per_cell =
1774  (2 * GeometryInfo<dim - 1>::quads_per_cell +
1775  GeometryInfo<dim - 1>::lines_per_cell);
1776 
1780  static constexpr unsigned int hexes_per_cell =
1781  (2 * GeometryInfo<dim - 1>::hexes_per_cell +
1782  GeometryInfo<dim - 1>::quads_per_cell);
1783 
1801  static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
1802  internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal;
1803 
1817  static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
1818  internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal;
1819 
1830  static constexpr std::array<std::array<unsigned int, dim>, vertices_per_cell>
1831  vertex_to_face =
1832  internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face;
1833 
1838  static unsigned int
1839  n_children(const RefinementCase<dim> &refinement_case);
1840 
1845  static unsigned int
1846  n_subfaces(const internal::SubfaceCase<dim> &subface_case);
1847 
1857  static double
1858  subface_ratio(const internal::SubfaceCase<dim> &subface_case,
1859  const unsigned int subface_no);
1860 
1866  static RefinementCase<dim - 1>
1867  face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
1868  const unsigned int face_no,
1869  const bool face_orientation = true,
1870  const bool face_flip = false,
1871  const bool face_rotation = false);
1872 
1878  static RefinementCase<dim>
1879  min_cell_refinement_case_for_face_refinement(
1880  const RefinementCase<dim - 1> &face_refinement_case,
1881  const unsigned int face_no,
1882  const bool face_orientation = true,
1883  const bool face_flip = false,
1884  const bool face_rotation = false);
1885 
1890  static RefinementCase<1>
1891  line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
1892  const unsigned int line_no);
1893 
1898  static RefinementCase<dim>
1899  min_cell_refinement_case_for_line_refinement(const unsigned int line_no);
1900 
1947  static unsigned int
1948  child_cell_on_face(const RefinementCase<dim> & ref_case,
1949  const unsigned int face,
1950  const unsigned int subface,
1951  const bool face_orientation = true,
1952  const bool face_flip = false,
1953  const bool face_rotation = false,
1954  const RefinementCase<dim - 1> &face_refinement_case =
1955  RefinementCase<dim - 1>::isotropic_refinement);
1956 
1970  static unsigned int
1971  line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
1972 
1993  static unsigned int
1994  face_to_cell_vertices(const unsigned int face,
1995  const unsigned int vertex,
1996  const bool face_orientation = true,
1997  const bool face_flip = false,
1998  const bool face_rotation = false);
1999 
2011  static unsigned int
2012  face_to_cell_lines(const unsigned int face,
2013  const unsigned int line,
2014  const bool face_orientation = true,
2015  const bool face_flip = false,
2016  const bool face_rotation = false);
2017 
2027  static unsigned int
2028  standard_to_real_face_vertex(const unsigned int vertex,
2029  const bool face_orientation = true,
2030  const bool face_flip = false,
2031  const bool face_rotation = false);
2032 
2042  static unsigned int
2043  real_to_standard_face_vertex(const unsigned int vertex,
2044  const bool face_orientation = true,
2045  const bool face_flip = false,
2046  const bool face_rotation = false);
2047 
2057  static unsigned int
2058  standard_to_real_face_line(const unsigned int line,
2059  const bool face_orientation = true,
2060  const bool face_flip = false,
2061  const bool face_rotation = false);
2062 
2072  static unsigned int
2073  real_to_standard_face_line(const unsigned int line,
2074  const bool face_orientation = true,
2075  const bool face_flip = false,
2076  const bool face_rotation = false);
2077 
2083  static Point<dim>
2084  unit_cell_vertex(const unsigned int vertex);
2085 
2095  static unsigned int
2096  child_cell_from_point(const Point<dim> &p);
2097 
2105  static Point<dim>
2106  cell_to_child_coordinates(const Point<dim> & p,
2107  const unsigned int child_index,
2108  const RefinementCase<dim> refine_case =
2109  RefinementCase<dim>::isotropic_refinement);
2110 
2116  static Point<dim>
2117  child_to_cell_coordinates(const Point<dim> & p,
2118  const unsigned int child_index,
2119  const RefinementCase<dim> refine_case =
2120  RefinementCase<dim>::isotropic_refinement);
2121 
2126  static bool
2127  is_inside_unit_cell(const Point<dim> &p);
2128 
2140  static bool
2141  is_inside_unit_cell(const Point<dim> &p, const double eps);
2142 
2147  static Point<dim>
2148  project_to_unit_cell(const Point<dim> &p);
2149 
2155  static double
2156  distance_to_unit_cell(const Point<dim> &p);
2157 
2162  static double
2163  d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2164 
2169  static Tensor<1, dim>
2170  d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2171 
2223  template <int spacedim>
2224  static void
2225  alternating_form_at_vertices
2226 #ifndef DEAL_II_CONSTEXPR_BUG
2227  (const Point<spacedim> (&vertices)[vertices_per_cell],
2228  Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2229 #else
2230  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2231 #endif
2232 
2242  static constexpr std::array<unsigned int, faces_per_cell>
2243  unit_normal_direction =
2244  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction;
2245 
2262  static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2263  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation;
2264 
2270  static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2271  internal::GeometryInfoHelper::Initializers<dim>::opposite_face;
2272 
2273 
2277  DeclException1(ExcInvalidCoordinate,
2278  double,
2279  << "The coordinates must satisfy 0 <= x_i <= 1, "
2280  << "but here we have x_i=" << arg1);
2281 
2285  DeclException3(ExcInvalidSubface,
2286  int,
2287  int,
2288  int,
2289  << "RefinementCase<dim> " << arg1 << ": face " << arg2
2290  << " has no subface " << arg3);
2291 };
2292 
2293 
2294 
2295 #ifndef DOXYGEN
2296 
2297 
2298 /* -------------- declaration of explicit specializations ------------- */
2299 
2300 
2301 template <>
2302 Tensor<1, 1>
2303 GeometryInfo<1>::d_linear_shape_function_gradient(const Point<1> & xi,
2304  const unsigned int i);
2305 template <>
2306 Tensor<1, 2>
2307 GeometryInfo<2>::d_linear_shape_function_gradient(const Point<2> & xi,
2308  const unsigned int i);
2309 template <>
2310 Tensor<1, 3>
2311 GeometryInfo<3>::d_linear_shape_function_gradient(const Point<3> & xi,
2312  const unsigned int i);
2313 
2314 
2315 
2316 /* -------------- inline functions ------------- */
2317 
2318 
2319 inline GeometryPrimitive::GeometryPrimitive(const Object object)
2320  : object(object)
2321 {}
2322 
2323 
2324 
2325 inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2326  : object(static_cast<Object>(object_dimension))
2327 {}
2328 
2329 
2330 inline GeometryPrimitive::operator unsigned int() const
2331 {
2332  return static_cast<unsigned int>(object);
2333 }
2334 
2335 
2336 
2337 namespace internal
2338 {
2339  template <int dim>
2340  inline SubfaceCase<dim>::SubfaceCase(
2341  const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2342  : value(subface_possibility)
2343  {}
2344 
2345 
2346  template <int dim>
2347  inline SubfaceCase<dim>::operator std::uint8_t() const
2348  {
2349  return value;
2350  }
2351 
2352 
2353 } // namespace internal
2354 
2355 
2356 template <int dim>
2357 inline RefinementCase<dim>
2358 RefinementCase<dim>::cut_axis(const unsigned int)
2359 {
2360  Assert(false, ExcInternalError());
2361  return static_cast<std::uint8_t>(-1);
2362 }
2363 
2364 
2365 template <>
2366 inline RefinementCase<1>
2367 RefinementCase<1>::cut_axis(const unsigned int i)
2368 {
2369  Assert(i < 1, ExcIndexRange(i, 0, 1));
2370 
2371  const RefinementCase options[1] = {RefinementPossibilities<1>::cut_x};
2372  return options[i];
2373 }
2374 
2375 
2376 
2377 template <>
2378 inline RefinementCase<2>
2379 RefinementCase<2>::cut_axis(const unsigned int i)
2380 {
2381  Assert(i < 2, ExcIndexRange(i, 0, 2));
2382 
2383  const RefinementCase options[2] = {RefinementPossibilities<2>::cut_x,
2384  RefinementPossibilities<2>::cut_y};
2385  return options[i];
2386 }
2387 
2388 
2389 
2390 template <>
2391 inline RefinementCase<3>
2392 RefinementCase<3>::cut_axis(const unsigned int i)
2393 {
2394  Assert(i < 3, ExcIndexRange(i, 0, 3));
2395 
2396  const RefinementCase options[3] = {RefinementPossibilities<3>::cut_x,
2397  RefinementPossibilities<3>::cut_y,
2398  RefinementPossibilities<3>::cut_z};
2399  return options[i];
2400 }
2401 
2402 
2403 
2404 template <int dim>
2405 inline RefinementCase<dim>::RefinementCase()
2406  : value(RefinementPossibilities<dim>::no_refinement)
2407 {}
2408 
2409 
2410 
2411 template <int dim>
2412 inline RefinementCase<dim>::RefinementCase(
2413  const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2414  : value(refinement_case)
2415 {
2416  // check that only those bits of
2417  // the given argument are set that
2418  // make sense for a given space
2419  // dimension
2420  Assert((refinement_case &
2421  RefinementPossibilities<dim>::isotropic_refinement) ==
2422  refinement_case,
2423  ExcInvalidRefinementCase(refinement_case));
2424 }
2425 
2426 
2427 
2428 template <int dim>
2429 inline RefinementCase<dim>::RefinementCase(const std::uint8_t refinement_case)
2430  : value(refinement_case)
2431 {
2432  // check that only those bits of
2433  // the given argument are set that
2434  // make sense for a given space
2435  // dimension
2436  Assert((refinement_case &
2437  RefinementPossibilities<dim>::isotropic_refinement) ==
2438  refinement_case,
2439  ExcInvalidRefinementCase(refinement_case));
2440 }
2441 
2442 
2443 
2444 template <int dim>
2445 inline RefinementCase<dim>::operator std::uint8_t() const
2446 {
2447  return value;
2448 }
2449 
2450 
2451 
2452 template <int dim>
2453 inline RefinementCase<dim>
2454 RefinementCase<dim>::operator|(const RefinementCase<dim> &r) const
2455 {
2456  return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2457 }
2458 
2459 
2460 
2461 template <int dim>
2462 inline RefinementCase<dim> RefinementCase<dim>::
2463  operator&(const RefinementCase<dim> &r) const
2464 {
2465  return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2466 }
2467 
2468 
2469 
2470 template <int dim>
2471 inline RefinementCase<dim>
2472 RefinementCase<dim>::operator~() const
2473 {
2474  return RefinementCase<dim>(static_cast<std::uint8_t>(
2475  (~value) & RefinementPossibilities<dim>::isotropic_refinement));
2476 }
2477 
2478 
2479 
2480 template <int dim>
2481 inline std::size_t
2482 RefinementCase<dim>::memory_consumption()
2483 {
2484  return sizeof(RefinementCase<dim>);
2485 }
2486 
2487 
2488 
2489 template <int dim>
2490 template <class Archive>
2491 inline void
2492 RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2493 {
2494  // serialization can't deal with bitfields, so copy from/to a full sized
2495  // std::uint8_t
2496  std::uint8_t uchar_value = value;
2497  ar & uchar_value;
2498  value = uchar_value;
2499 }
2500 
2501 
2502 
2503 template <>
2504 inline Point<1>
2505 GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2506 {
2507  Assert(vertex < vertices_per_cell,
2508  ExcIndexRange(vertex, 0, vertices_per_cell));
2509 
2510  return Point<1>(static_cast<double>(vertex));
2511 }
2512 
2513 
2514 
2515 template <>
2516 inline Point<2>
2517 GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2518 {
2519  Assert(vertex < vertices_per_cell,
2520  ExcIndexRange(vertex, 0, vertices_per_cell));
2521 
2522  return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2523 }
2524 
2525 
2526 
2527 template <>
2528 inline Point<3>
2529 GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2530 {
2531  Assert(vertex < vertices_per_cell,
2532  ExcIndexRange(vertex, 0, vertices_per_cell));
2533 
2534  return {static_cast<double>(vertex % 2),
2535  static_cast<double>(vertex / 2 % 2),
2536  static_cast<double>(vertex / 4)};
2537 }
2538 
2539 
2540 
2541 template <int dim>
2542 inline Point<dim>
2543 GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2544 {
2545  Assert(false, ExcNotImplemented());
2546 
2547  return Point<dim>();
2548 }
2549 
2550 
2551 
2552 template <>
2553 inline unsigned int
2554 GeometryInfo<1>::child_cell_from_point(const Point<1> &p)
2555 {
2556  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2557 
2558  return (p[0] <= 0.5 ? 0 : 1);
2559 }
2560 
2561 
2562 
2563 template <>
2564 inline unsigned int
2565 GeometryInfo<2>::child_cell_from_point(const Point<2> &p)
2566 {
2567  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2568  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2569 
2570  return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
2571 }
2572 
2573 
2574 
2575 template <>
2576 inline unsigned int
2577 GeometryInfo<3>::child_cell_from_point(const Point<3> &p)
2578 {
2579  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2580  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2581  Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
2582 
2583  return (p[0] <= 0.5 ?
2584  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
2585  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
2586 }
2587 
2588 
2589 template <int dim>
2590 inline unsigned int
2591 GeometryInfo<dim>::child_cell_from_point(const Point<dim> &)
2592 {
2593  Assert(false, ExcNotImplemented());
2594 
2595  return 0;
2596 }
2597 
2598 
2599 
2600 template <>
2601 inline Point<1>
2602 GeometryInfo<1>::cell_to_child_coordinates(const Point<1> & p,
2603  const unsigned int child_index,
2604  const RefinementCase<1> refine_case)
2605 
2606 {
2607  Assert(child_index < 2, ExcIndexRange(child_index, 0, 2));
2608  Assert(refine_case == RefinementCase<1>::cut_x, ExcInternalError());
2609  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
2610 
2611  return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
2612 }
2613 
2614 
2615 
2616 template <>
2617 inline Point<2>
2618 GeometryInfo<2>::cell_to_child_coordinates(const Point<2> & p,
2619  const unsigned int child_index,
2620  const RefinementCase<2> refine_case)
2621 
2622 {
2623  Assert(child_index < GeometryInfo<2>::n_children(refine_case),
2624  ExcIndexRange(child_index,
2625  0,
2626  GeometryInfo<2>::n_children(refine_case)));
2627 
2628  Point<2> point = p;
2629  switch (refine_case)
2630  {
2631  case RefinementCase<2>::cut_x:
2632  point[0] *= 2.0;
2633  if (child_index == 1)
2634  point[0] -= 1.0;
2635  break;
2636  case RefinementCase<2>::cut_y:
2637  point[1] *= 2.0;
2638  if (child_index == 1)
2639  point[1] -= 1.0;
2640  break;
2641  case RefinementCase<2>::cut_xy:
2642  point *= 2.0;
2643  point -= unit_cell_vertex(child_index);
2644  break;
2645  default:
2646  Assert(false, ExcInternalError());
2647  }
2648 
2649  return point;
2650 }
2651 
2652 
2653 
2654 template <>
2655 inline Point<3>
2656 GeometryInfo<3>::cell_to_child_coordinates(const Point<3> & p,
2657  const unsigned int child_index,
2658  const RefinementCase<3> refine_case)
2659 
2660 {
2661  Assert(child_index < GeometryInfo<3>::n_children(refine_case),
2662  ExcIndexRange(child_index,
2663  0,
2664  GeometryInfo<3>::n_children(refine_case)));
2665 
2666  Point<3> point = p;
2667  // there might be a cleverer way to do
2668  // this, but since this function is called
2669  // in very few places for initialization
2670  // purposes only, I don't care at the
2671  // moment
2672  switch (refine_case)
2673  {
2674  case RefinementCase<3>::cut_x:
2675  point[0] *= 2.0;
2676  if (child_index == 1)
2677  point[0] -= 1.0;
2678  break;
2679  case RefinementCase<3>::cut_y:
2680  point[1] *= 2.0;
2681  if (child_index == 1)
2682  point[1] -= 1.0;
2683  break;
2684  case RefinementCase<3>::cut_z:
2685  point[2] *= 2.0;
2686  if (child_index == 1)
2687  point[2] -= 1.0;
2688  break;
2689  case RefinementCase<3>::cut_xy:
2690  point[0] *= 2.0;
2691  point[1] *= 2.0;
2692  if (child_index % 2 == 1)
2693  point[0] -= 1.0;
2694  if (child_index / 2 == 1)
2695  point[1] -= 1.0;
2696  break;
2697  case RefinementCase<3>::cut_xz:
2698  // careful, this is slightly
2699  // different from xy and yz due to
2700  // different internal numbering of
2701  // children!
2702  point[0] *= 2.0;
2703  point[2] *= 2.0;
2704  if (child_index / 2 == 1)
2705  point[0] -= 1.0;
2706  if (child_index % 2 == 1)
2707  point[2] -= 1.0;
2708  break;
2709  case RefinementCase<3>::cut_yz:
2710  point[1] *= 2.0;
2711  point[2] *= 2.0;
2712  if (child_index % 2 == 1)
2713  point[1] -= 1.0;
2714  if (child_index / 2 == 1)
2715  point[2] -= 1.0;
2716  break;
2717  case RefinementCase<3>::cut_xyz:
2718  point *= 2.0;
2719  point -= unit_cell_vertex(child_index);
2720  break;
2721  default:
2722  Assert(false, ExcInternalError());
2723  }
2724 
2725  return point;
2726 }
2727 
2728 
2729 
2730 template <int dim>
2731 inline Point<dim>
2732 GeometryInfo<dim>::cell_to_child_coordinates(
2733  const Point<dim> & /*p*/,
2734  const unsigned int /*child_index*/,
2735  const RefinementCase<dim> /*refine_case*/)
2736 
2737 {
2738  Assert(false, ExcNotImplemented());
2739  return Point<dim>();
2740 }
2741 
2742 
2743 
2744 template <>
2745 inline Point<1>
2746 GeometryInfo<1>::child_to_cell_coordinates(const Point<1> & p,
2747  const unsigned int child_index,
2748  const RefinementCase<1> refine_case)
2749 
2750 {
2751  Assert(child_index < 2, ExcIndexRange(child_index, 0, 2));
2752  Assert(refine_case == RefinementCase<1>::cut_x, ExcInternalError());
2753  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
2754 
2755  return (p + unit_cell_vertex(child_index)) * 0.5;
2756 }
2757 
2758 
2759 
2760 template <>
2761 inline Point<3>
2762 GeometryInfo<3>::child_to_cell_coordinates(const Point<3> & p,
2763  const unsigned int child_index,
2764  const RefinementCase<3> refine_case)
2765 
2766 {
2767  Assert(child_index < GeometryInfo<3>::n_children(refine_case),
2768  ExcIndexRange(child_index,
2769  0,
2770  GeometryInfo<3>::n_children(refine_case)));
2771 
2772  Point<3> point = p;
2773  // there might be a cleverer way to do
2774  // this, but since this function is called
2775  // in very few places for initialization
2776  // purposes only, I don't care at the
2777  // moment
2778  switch (refine_case)
2779  {
2780  case RefinementCase<3>::cut_x:
2781  if (child_index == 1)
2782  point[0] += 1.0;
2783  point[0] *= 0.5;
2784  break;
2785  case RefinementCase<3>::cut_y:
2786  if (child_index == 1)
2787  point[1] += 1.0;
2788  point[1] *= 0.5;
2789  break;
2790  case RefinementCase<3>::cut_z:
2791  if (child_index == 1)
2792  point[2] += 1.0;
2793  point[2] *= 0.5;
2794  break;
2795  case RefinementCase<3>::cut_xy:
2796  if (child_index % 2 == 1)
2797  point[0] += 1.0;
2798  if (child_index / 2 == 1)
2799  point[1] += 1.0;
2800  point[0] *= 0.5;
2801  point[1] *= 0.5;
2802  break;
2803  case RefinementCase<3>::cut_xz:
2804  // careful, this is slightly
2805  // different from xy and yz due to
2806  // different internal numbering of
2807  // children!
2808  if (child_index / 2 == 1)
2809  point[0] += 1.0;
2810  if (child_index % 2 == 1)
2811  point[2] += 1.0;
2812  point[0] *= 0.5;
2813  point[2] *= 0.5;
2814  break;
2815  case RefinementCase<3>::cut_yz:
2816  if (child_index % 2 == 1)
2817  point[1] += 1.0;
2818  if (child_index / 2 == 1)
2819  point[2] += 1.0;
2820  point[1] *= 0.5;
2821  point[2] *= 0.5;
2822  break;
2823  case RefinementCase<3>::cut_xyz:
2824  point += unit_cell_vertex(child_index);
2825  point *= 0.5;
2826  break;
2827  default:
2828  Assert(false, ExcInternalError());
2829  }
2830 
2831  return point;
2832 }
2833 
2834 
2835 
2836 template <>
2837 inline Point<2>
2838 GeometryInfo<2>::child_to_cell_coordinates(const Point<2> & p,
2839  const unsigned int child_index,
2840  const RefinementCase<2> refine_case)
2841 {
2842  Assert(child_index < GeometryInfo<2>::n_children(refine_case),
2843  ExcIndexRange(child_index,
2844  0,
2845  GeometryInfo<2>::n_children(refine_case)));
2846 
2847  Point<2> point = p;
2848  switch (refine_case)
2849  {
2850  case RefinementCase<2>::cut_x:
2851  if (child_index == 1)
2852  point[0] += 1.0;
2853  point[0] *= 0.5;
2854  break;
2855  case RefinementCase<2>::cut_y:
2856  if (child_index == 1)
2857  point[1] += 1.0;
2858  point[1] *= 0.5;
2859  break;
2860  case RefinementCase<2>::cut_xy:
2861  point += unit_cell_vertex(child_index);
2862  point *= 0.5;
2863  break;
2864  default:
2865  Assert(false, ExcInternalError());
2866  }
2867 
2868  return point;
2869 }
2870 
2871 
2872 
2873 template <int dim>
2874 inline Point<dim>
2875 GeometryInfo<dim>::child_to_cell_coordinates(
2876  const Point<dim> & /*p*/,
2877  const unsigned int /*child_index*/,
2878  const RefinementCase<dim> /*refine_case*/)
2879 {
2880  Assert(false, ExcNotImplemented());
2881  return Point<dim>();
2882 }
2883 
2884 
2885 
2886 template <int dim>
2887 inline bool
2888 GeometryInfo<dim>::is_inside_unit_cell(const Point<dim> &)
2889 {
2890  Assert(false, ExcNotImplemented());
2891  return false;
2892 }
2893 
2894 template <>
2895 inline bool
2896 GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p)
2897 {
2898  return (p[0] >= 0.) && (p[0] <= 1.);
2899 }
2900 
2901 
2902 
2903 template <>
2904 inline bool
2905 GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p)
2906 {
2907  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
2908 }
2909 
2910 
2911 
2912 template <>
2913 inline bool
2914 GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p)
2915 {
2916  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
2917  (p[2] >= 0.) && (p[2] <= 1.);
2918 }
2919 
2920 
2921 
2922 template <int dim>
2923 inline bool
2924 GeometryInfo<dim>::is_inside_unit_cell(const Point<dim> &, const double)
2925 {
2926  Assert(false, ExcNotImplemented());
2927  return false;
2928 }
2929 
2930 template <>
2931 inline bool
2932 GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
2933 {
2934  return (p[0] >= -eps) && (p[0] <= 1. + eps);
2935 }
2936 
2937 
2938 
2939 template <>
2940 inline bool
2941 GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
2942 {
2943  const double l = -eps, u = 1 + eps;
2944  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
2945 }
2946 
2947 
2948 
2949 template <>
2950 inline bool
2951 GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
2952 {
2953  const double l = -eps, u = 1.0 + eps;
2954  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
2955  (p[2] >= l) && (p[2] <= u);
2956 }
2957 
2958 
2959 
2960 template <>
2961 inline unsigned int
2962 GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
2963  const unsigned int vertex)
2964 {
2965  (void)line;
2966  Assert(line < lines_per_cell, ExcIndexRange(line, 0, lines_per_cell));
2967  Assert(vertex < 2, ExcIndexRange(vertex, 0, 2));
2968 
2969  return vertex;
2970 }
2971 
2972 
2973 template <>
2974 inline unsigned int
2975 GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
2976  const unsigned int vertex)
2977 {
2978  constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
2979  return cell_vertices[line][vertex];
2980 }
2981 
2982 
2983 
2984 template <>
2985 inline unsigned int
2986 GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
2987  const unsigned int vertex)
2988 {
2989  Assert(line < lines_per_cell, ExcIndexRange(line, 0, lines_per_cell));
2990  Assert(vertex < 2, ExcIndexRange(vertex, 0, 2));
2991 
2992  constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
2993  {1, 3},
2994  {0, 1},
2995  {2, 3},
2996  {4, 6}, // top face
2997  {5, 7},
2998  {4, 5},
2999  {6, 7},
3000  {0,
3001  4}, // connects of bottom
3002  {1, 5}, // top face
3003  {2, 6},
3004  {3, 7}};
3005  return vertices[line][vertex];
3006 }
3007 
3008 
3009 template <>
3010 inline unsigned int
3011 GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3012 {
3013  Assert(false, ExcNotImplemented());
3014  return numbers::invalid_unsigned_int;
3015 }
3016 
3017 template <>
3018 inline unsigned int
3019 GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3020  const bool face_orientation,
3021  const bool face_flip,
3022  const bool face_rotation)
3023 {
3024  Assert(vertex < GeometryInfo<3>::vertices_per_face,
3025  ExcIndexRange(vertex, 0, GeometryInfo<3>::vertices_per_face));
3026 
3027  // set up a table to make sure that
3028  // we handle non-standard faces correctly
3029  //
3030  // so set up a table that for each vertex (of
3031  // a quad in standard position) describes
3032  // which vertex to take
3033  //
3034  // first index: four vertices 0...3
3035  //
3036  // second index: face_orientation; 0:
3037  // opposite normal, 1: standard
3038  //
3039  // third index: face_flip; 0: standard, 1:
3040  // face rotated by 180 degrees
3041  //
3042  // forth index: face_rotation: 0: standard,
3043  // 1: face rotated by 90 degrees
3044 
3045  constexpr unsigned int vertex_translation[4][2][2][2] = {
3046  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3047  // face_rotation=false and true
3048  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3049  // face_rotation=false and true
3050  {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3051  // face_rotation=false and true
3052  {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3053  // face_rotation=false and true
3054 
3055  {{{2, 3}, // vertex 1 ...
3056  {1, 0}},
3057  {{1, 0}, {2, 3}}},
3058 
3059  {{{1, 0}, // vertex 2 ...
3060  {2, 3}},
3061  {{2, 3}, {1, 0}}},
3062 
3063  {{{3, 1}, // vertex 3 ...
3064  {0, 2}},
3065  {{3, 1}, {0, 2}}}};
3066 
3067  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3068 }
3069 
3070 
3071 
3072 template <int dim>
3073 inline unsigned int
3074 GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3075  const bool,
3076  const bool,
3077  const bool)
3078 {
3079  Assert(dim > 1, ExcImpossibleInDim(dim));
3080  Assert(vertex < GeometryInfo<dim>::vertices_per_face,
3081  ExcIndexRange(vertex, 0, GeometryInfo<dim>::vertices_per_face));
3082  return vertex;
3083 }
3084 
3085 template <int dim>
3086 inline unsigned int
3087 GeometryInfo<dim>::n_children(const RefinementCase<dim> &ref_case)
3088 {
3089  constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3090  0, 2, 2, 4, 2, 4, 4, 8};
3091 
3092  return n_children[ref_case];
3093 }
3094 
3095 
3096 
3097 template <int dim>
3098 inline unsigned int
3099 GeometryInfo<dim>::n_subfaces(const internal::SubfaceCase<dim> &)
3100 {
3101  Assert(false, ExcNotImplemented());
3102  return 0;
3103 }
3104 
3105 template <>
3106 inline unsigned int
3107 GeometryInfo<1>::n_subfaces(const internal::SubfaceCase<1> &)
3108 {
3109  Assert(false, ExcImpossibleInDim(1));
3110  return 0;
3111 }
3112 
3113 template <>
3114 inline unsigned int
3115 GeometryInfo<2>::n_subfaces(const internal::SubfaceCase<2> &subface_case)
3116 {
3117  return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3118 }
3119 
3120 
3121 
3122 template <>
3123 inline unsigned int
3124 GeometryInfo<3>::n_subfaces(const internal::SubfaceCase<3> &subface_case)
3125 {
3126  const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3127  0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3128  return nsubs[subface_case];
3129 }
3130 
3131 
3132 
3133 template <int dim>
3134 inline double
3135 GeometryInfo<dim>::subface_ratio(const internal::SubfaceCase<dim> &,
3136  const unsigned int)
3137 {
3138  Assert(false, ExcNotImplemented());
3139  return 0.;
3140 }
3141 
3142 template <>
3143 inline double
3144 GeometryInfo<1>::subface_ratio(const internal::SubfaceCase<1> &,
3145  const unsigned int)
3146 {
3147  return 1;
3148 }
3149 
3150 
3151 template <>
3152 inline double
3153 GeometryInfo<2>::subface_ratio(const internal::SubfaceCase<2> &subface_case,
3154  const unsigned int)
3155 {
3156  double ratio = 1;
3157  switch (subface_case)
3158  {
3159  case internal::SubfaceCase<2>::case_none:
3160  // Here, an
3161  // Assert(false,ExcInternalError())
3162  // would be the right
3163  // choice, but
3164  // unfortunately the
3165  // current function is
3166  // also called for faces
3167  // without children (see
3168  // tests/fe/mapping.cc).
3169  // Assert(false, ExcMessage("Face has no subfaces."));
3170  // Furthermore, assign
3171  // following value as
3172  // otherwise the
3173  // bits/volume_x tests
3174  // break
3175  ratio = 1. / GeometryInfo<2>::max_children_per_face;
3176  break;
3177  case internal::SubfaceCase<2>::case_x:
3178  ratio = 0.5;
3179  break;
3180  default:
3181  // there should be no
3182  // cases left
3183  Assert(false, ExcInternalError());
3184  break;
3185  }
3186 
3187  return ratio;
3188 }
3189 
3190 
3191 template <>
3192 inline double
3193 GeometryInfo<3>::subface_ratio(const internal::SubfaceCase<3> &subface_case,
3194  const unsigned int subface_no)
3195 {
3196  double ratio = 1;
3197  switch (subface_case)
3198  {
3199  case internal::SubfaceCase<3>::case_none:
3200  // Here, an
3201  // Assert(false,ExcInternalError())
3202  // would be the right
3203  // choice, but
3204  // unfortunately the
3205  // current function is
3206  // also called for faces
3207  // without children (see
3208  // tests/bits/mesh_3d_16.cc). Add
3209  // following switch to
3210  // avoid diffs in
3211  // tests/bits/mesh_3d_16
3212  ratio = 1. / GeometryInfo<3>::max_children_per_face;
3213  break;
3214  case internal::SubfaceCase<3>::case_x:
3215  case internal::SubfaceCase<3>::case_y:
3216  ratio = 0.5;
3217  break;
3218  case internal::SubfaceCase<3>::case_xy:
3219  case internal::SubfaceCase<3>::case_x1y2y:
3220  case internal::SubfaceCase<3>::case_y1x2x:
3221  ratio = 0.25;
3222  break;
3223  case internal::SubfaceCase<3>::case_x1y:
3224  case internal::SubfaceCase<3>::case_y1x:
3225  if (subface_no < 2)
3226  ratio = 0.25;
3227  else
3228  ratio = 0.5;
3229  break;
3230  case internal::SubfaceCase<3>::case_x2y:
3231  case internal::SubfaceCase<3>::case_y2x:
3232  if (subface_no == 0)
3233  ratio = 0.5;
3234  else
3235  ratio = 0.25;
3236  break;
3237  default:
3238  // there should be no
3239  // cases left
3240  Assert(false, ExcInternalError());
3241  break;
3242  }
3243 
3244  return ratio;
3245 }
3246 
3247 
3248 
3249 template <int dim>
3250 RefinementCase<dim - 1> inline GeometryInfo<dim>::face_refinement_case(
3251  const RefinementCase<dim> &,
3252  const unsigned int,
3253  const bool,
3254  const bool,
3255  const bool)
3256 {
3257  Assert(false, ExcNotImplemented());
3258  return RefinementCase<dim - 1>::no_refinement;
3259 }
3260 
3261 template <>
3262 RefinementCase<0> inline GeometryInfo<1>::face_refinement_case(
3263  const RefinementCase<1> &,
3264  const unsigned int,
3265  const bool,
3266  const bool,
3267  const bool)
3268 {
3269  Assert(false, ExcImpossibleInDim(1));
3270 
3271  return RefinementCase<0>::no_refinement;
3272 }
3273 
3274 
3275 template <>
3276 inline RefinementCase<1>
3277 GeometryInfo<2>::face_refinement_case(
3278  const RefinementCase<2> &cell_refinement_case,
3279  const unsigned int face_no,
3280  const bool,
3281  const bool,
3282  const bool)
3283 {
3284  const unsigned int dim = 2;
3285  Assert(cell_refinement_case < RefinementCase<dim>::isotropic_refinement + 1,
3286  ExcIndexRange(cell_refinement_case,
3287  0,
3288  RefinementCase<dim>::isotropic_refinement + 1));
3289  Assert(face_no < GeometryInfo<dim>::faces_per_cell,
3290  ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
3291 
3292  const RefinementCase<dim - 1>
3293  ref_cases[RefinementCase<dim>::isotropic_refinement +
3294  1][GeometryInfo<dim>::faces_per_cell / 2] = {
3295  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3296  RefinementCase<dim - 1>::no_refinement},
3297 
3298  {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3299 
3300  {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3301 
3302  {RefinementCase<dim - 1>::cut_x, // cut_xy
3303  RefinementCase<dim - 1>::cut_x}};
3304 
3305  return ref_cases[cell_refinement_case][face_no / 2];
3306 }
3307 
3308 
3309 template <>
3310 inline RefinementCase<2>
3311 GeometryInfo<3>::face_refinement_case(
3312  const RefinementCase<3> &cell_refinement_case,
3313  const unsigned int face_no,
3314  const bool face_orientation,
3315  const bool /*face_flip*/,
3316  const bool face_rotation)
3317 {
3318  const unsigned int dim = 3;
3319  Assert(cell_refinement_case < RefinementCase<dim>::isotropic_refinement + 1,
3320  ExcIndexRange(cell_refinement_case,
3321  0,
3322  RefinementCase<dim>::isotropic_refinement + 1));
3323  Assert(face_no < GeometryInfo<dim>::faces_per_cell,
3324  ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
3325 
3326  const RefinementCase<dim - 1>
3327  ref_cases[RefinementCase<dim>::isotropic_refinement + 1]
3328  [GeometryInfo<dim>::faces_per_cell / 2] = {
3329  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3330  RefinementCase<dim - 1>::no_refinement,
3331  RefinementCase<dim - 1>::no_refinement},
3332 
3333  {RefinementCase<dim - 1>::no_refinement, // cut_x
3334  RefinementCase<dim - 1>::cut_y,
3335  RefinementCase<dim - 1>::cut_x},
3336 
3337  {RefinementCase<dim - 1>::cut_x, // cut_y
3338  RefinementCase<dim - 1>::no_refinement,
3339  RefinementCase<dim - 1>::cut_y},
3340 
3341  {RefinementCase<dim - 1>::cut_x, // cut_xy
3342  RefinementCase<dim - 1>::cut_y,
3343  RefinementCase<dim - 1>::cut_xy},
3344 
3345  {RefinementCase<dim - 1>::cut_y, // cut_z
3346  RefinementCase<dim - 1>::cut_x,
3347  RefinementCase<dim - 1>::no_refinement},
3348 
3349  {RefinementCase<dim - 1>::cut_y, // cut_xz
3350  RefinementCase<dim - 1>::cut_xy,
3351  RefinementCase<dim - 1>::cut_x},
3352 
3353  {RefinementCase<dim - 1>::cut_xy, // cut_yz
3354  RefinementCase<dim - 1>::cut_x,
3355  RefinementCase<dim - 1>::cut_y},
3356 
3357  {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3358  RefinementCase<dim - 1>::cut_xy,
3359  RefinementCase<dim - 1>::cut_xy},
3360  };
3361 
3362  const RefinementCase<dim - 1> ref_case =
3363  ref_cases[cell_refinement_case][face_no / 2];
3364 
3365  const RefinementCase<dim - 1> flip[4] = {
3366  RefinementCase<dim - 1>::no_refinement,
3367  RefinementCase<dim - 1>::cut_y,
3368  RefinementCase<dim - 1>::cut_x,
3369  RefinementCase<dim - 1>::cut_xy};
3370 
3371  // correct the ref_case for face_orientation
3372  // and face_rotation. for face_orientation,
3373  // 'true' is the default value whereas for
3374  // face_rotation, 'false' is standard. If
3375  // <tt>face_rotation==face_orientation</tt>,
3376  // then one of them is non-standard and we
3377  // have to swap cut_x and cut_y, otherwise no
3378  // change is necessary. face_flip has no
3379  // influence. however, in order to keep the
3380  // interface consistent with other functions,
3381  // we still include it as an argument to this
3382  // function
3383  return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3384 }
3385 
3386 
3387 
3388 template <int dim>
3389 inline RefinementCase<1>
3390 GeometryInfo<dim>::line_refinement_case(const RefinementCase<dim> &,
3391  const unsigned int)
3392 {
3393  Assert(false, ExcNotImplemented());
3394  return RefinementCase<1>::no_refinement;
3395 }
3396 
3397 template <>
3398 inline RefinementCase<1>
3399 GeometryInfo<1>::line_refinement_case(
3400  const RefinementCase<1> &cell_refinement_case,
3401  const unsigned int line_no)
3402 {
3403  (void)line_no;
3404  const unsigned int dim = 1;
3405  (void)dim;
3406  Assert(cell_refinement_case < RefinementCase<dim>::isotropic_refinement + 1,
3407  ExcIndexRange(cell_refinement_case,
3408  0,
3409  RefinementCase<dim>::isotropic_refinement + 1));
3410  Assert(line_no < GeometryInfo<dim>::lines_per_cell,
3411  ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
3412 
3413  return cell_refinement_case;
3414 }
3415 
3416 
3417 template <>
3418 inline RefinementCase<1>
3419 GeometryInfo<2>::line_refinement_case(
3420  const RefinementCase<2> &cell_refinement_case,
3421  const unsigned int line_no)
3422 {
3423  // Assertions are in face_refinement_case()
3424  return face_refinement_case(cell_refinement_case, line_no);
3425 }
3426 
3427 
3428 template <>
3429 inline RefinementCase<1>
3430 GeometryInfo<3>::line_refinement_case(
3431  const RefinementCase<3> &cell_refinement_case,
3432  const unsigned int line_no)
3433 {
3434  const unsigned int dim = 3;
3435  Assert(cell_refinement_case < RefinementCase<dim>::isotropic_refinement + 1,
3436  ExcIndexRange(cell_refinement_case,
3437  0,
3438  RefinementCase<dim>::isotropic_refinement + 1));
3439  Assert(line_no < GeometryInfo<dim>::lines_per_cell,
3440  ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
3441 
3442  // array indicating, which simple refine
3443  // case cuts a line in direction x, y or
3444  // z. For example, cut_y and everything
3445  // containing cut_y (cut_xy, cut_yz,
3446  // cut_xyz) cuts lines, which are in y
3447  // direction.
3448  const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3449  RefinementCase<dim>::cut_y,
3450  RefinementCase<dim>::cut_z};
3451 
3452  // order the direction of lines
3453  // 0->x, 1->y, 2->z
3454  const unsigned int direction[lines_per_cell] = {
3455  1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3456 
3457  return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3458  RefinementCase<1>::cut_x :
3459  RefinementCase<1>::no_refinement);
3460 }
3461 
3462 
3463 
3464 template <int dim>
3465 inline RefinementCase<dim>
3466 GeometryInfo<dim>::min_cell_refinement_case_for_face_refinement(
3467  const RefinementCase<dim - 1> &,
3468  const unsigned int,
3469  const bool,
3470  const bool,
3471  const bool)
3472 {
3473  Assert(false, ExcNotImplemented());
3474 
3475  return RefinementCase<dim>::no_refinement;
3476 }
3477 
3478 template <>
3479 inline RefinementCase<1>
3480 GeometryInfo<1>::min_cell_refinement_case_for_face_refinement(
3481  const RefinementCase<0> &,
3482  const unsigned int,
3483  const bool,
3484  const bool,
3485  const bool)
3486 {
3487  const unsigned int dim = 1;
3488  Assert(false, ExcImpossibleInDim(dim));
3489 
3490  return RefinementCase<dim>::no_refinement;
3491 }
3492 
3493 
3494 template <>
3495 inline RefinementCase<2>
3496 GeometryInfo<2>::min_cell_refinement_case_for_face_refinement(
3497  const RefinementCase<1> &face_refinement_case,
3498  const unsigned int face_no,
3499  const bool,
3500  const bool,
3501  const bool)
3502 {
3503  const unsigned int dim = 2;
3504  Assert(face_refinement_case <
3505  RefinementCase<dim - 1>::isotropic_refinement + 1,
3506  ExcIndexRange(face_refinement_case,
3507  0,
3508  RefinementCase<dim - 1>::isotropic_refinement + 1));
3509  Assert(face_no < GeometryInfo<dim>::faces_per_cell,
3510  ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
3511 
3512  if (face_refinement_case == RefinementCase<dim>::cut_x)
3513  return (face_no / 2) ? RefinementCase<dim>::cut_x :
3514  RefinementCase<dim>::cut_y;
3515  else
3516  return RefinementCase<dim>::no_refinement;
3517 }
3518 
3519 
3520 template <>
3521 inline RefinementCase<3>
3522 GeometryInfo<3>::min_cell_refinement_case_for_face_refinement(
3523  const RefinementCase<2> &face_refinement_case,
3524  const unsigned int face_no,
3525  const bool face_orientation,
3526  const bool /*face_flip*/,
3527  const bool face_rotation)
3528 {
3529  const unsigned int dim = 3;
3530  Assert(face_refinement_case <
3531  RefinementCase<dim - 1>::isotropic_refinement + 1,
3532  ExcIndexRange(face_refinement_case,
3533  0,
3534  RefinementCase<dim - 1>::isotropic_refinement + 1));
3535  Assert(face_no < GeometryInfo<dim>::faces_per_cell,
3536  ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
3537 
3538  const RefinementCase<2> flip[4] = {RefinementCase<2>::no_refinement,
3539  RefinementCase<2>::cut_y,
3540  RefinementCase<2>::cut_x,
3541  RefinementCase<2>::cut_xy};
3542 
3543  // correct the face_refinement_case for
3544  // face_orientation and face_rotation. for
3545  // face_orientation, 'true' is the default
3546  // value whereas for face_rotation, 'false'
3547  // is standard. If
3548  // <tt>face_rotation==face_orientation</tt>,
3549  // then one of them is non-standard and we
3550  // have to swap cut_x and cut_y, otherwise no
3551  // change is necessary. face_flip has no
3552  // influence. however, in order to keep the
3553  // interface consistent with other functions,
3554  // we still include it as an argument to this
3555  // function
3556  const RefinementCase<dim - 1> std_face_ref =
3557  (face_orientation == face_rotation) ? flip[face_refinement_case] :
3558  face_refinement_case;
3559 
3560  const RefinementCase<dim> face_to_cell[3][4] = {
3561  {RefinementCase<dim>::no_refinement, // faces 0 and 1
3562  RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3563  RefinementCase<dim>::cut_z,
3564  RefinementCase<dim>::cut_yz},
3565 
3566  {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3567  // "exchanged on faces 2 and 3")
3568  RefinementCase<dim>::cut_z,
3569  RefinementCase<dim>::cut_x,
3570  RefinementCase<dim>::cut_xz},
3571 
3572  {RefinementCase<dim>::no_refinement, // faces 4 and 5
3573  RefinementCase<dim>::cut_x,
3574  RefinementCase<dim>::cut_y,
3575  RefinementCase<dim>::cut_xy}};
3576 
3577  return face_to_cell[face_no / 2][std_face_ref];
3578 }
3579 
3580 
3581 
3582 template <int dim>
3583 inline RefinementCase<dim>
3584 GeometryInfo<dim>::min_cell_refinement_case_for_line_refinement(
3585  const unsigned int)
3586 {
3587  Assert(false, ExcNotImplemented());
3588 
3589  return RefinementCase<dim>::no_refinement;
3590 }
3591 
3592 template <>
3593 inline RefinementCase<1>
3594 GeometryInfo<1>::min_cell_refinement_case_for_line_refinement(
3595  const unsigned int line_no)
3596 {
3597  (void)line_no;
3598  Assert(line_no == 0, ExcIndexRange(line_no, 0, 1));
3599 
3600  return RefinementCase<1>::cut_x;
3601 }
3602 
3603 
3604 template <>
3605 inline RefinementCase<2>
3606 GeometryInfo<2>::min_cell_refinement_case_for_line_refinement(
3607  const unsigned int line_no)
3608 {
3609  const unsigned int dim = 2;
3610  (void)dim;
3611  Assert(line_no < GeometryInfo<dim>::lines_per_cell,
3612  ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
3613 
3614  return (line_no / 2) ? RefinementCase<2>::cut_x : RefinementCase<2>::cut_y;
3615 }
3616 
3617 
3618 template <>
3619 inline RefinementCase<3>
3620 GeometryInfo<3>::min_cell_refinement_case_for_line_refinement(
3621  const unsigned int line_no)
3622 {
3623  const unsigned int dim = 3;
3624  Assert(line_no < GeometryInfo<dim>::lines_per_cell,
3625  ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
3626 
3627  const RefinementCase<dim> ref_cases[6] = {
3628  RefinementCase<dim>::cut_y, // lines 0 and 1
3629  RefinementCase<dim>::cut_x, // lines 2 and 3
3630  RefinementCase<dim>::cut_y, // lines 4 and 5
3631  RefinementCase<dim>::cut_x, // lines 6 and 7
3632  RefinementCase<dim>::cut_z, // lines 8 and 9
3633  RefinementCase<dim>::cut_z}; // lines 10 and 11
3634 
3635  return ref_cases[line_no / 2];
3636 }
3637 
3638 
3639 
3640 template <>
3641 inline unsigned int
3642 GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
3643  const bool face_orientation,
3644  const bool face_flip,
3645  const bool face_rotation)
3646 {
3647  Assert(vertex < GeometryInfo<3>::vertices_per_face,
3648  ExcIndexRange(vertex, 0, GeometryInfo<3>::vertices_per_face));
3649 
3650  // set up a table to make sure that
3651  // we handle non-standard faces correctly
3652  //
3653  // so set up a table that for each vertex (of
3654  // a quad in standard position) describes
3655  // which vertex to take
3656  //
3657  // first index: four vertices 0...3
3658  //
3659  // second index: face_orientation; 0:
3660  // opposite normal, 1: standard
3661  //
3662  // third index: face_flip; 0: standard, 1:
3663  // face rotated by 180 degrees
3664  //
3665  // forth index: face_rotation: 0: standard,
3666  // 1: face rotated by 90 degrees
3667 
3668  const unsigned int vertex_translation[4][2][2][2] = {
3669  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3670  // face_rotation=false and true
3671  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3672  // face_rotation=false and true
3673  {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
3674  // face_rotation=false and true
3675  {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
3676  // face_rotation=false and true
3677 
3678  {{{2, 3}, // vertex 1 ...
3679  {1, 0}},
3680  {{1, 3}, {2, 0}}},
3681 
3682  {{{1, 0}, // vertex 2 ...
3683  {2, 3}},
3684  {{2, 0}, {1, 3}}},
3685 
3686  {{{3, 1}, // vertex 3 ...
3687  {0, 2}},
3688  {{3, 2}, {0, 1}}}};
3689 
3690  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3691 }
3692 
3693 
3694 
3695 template <int dim>
3696 inline unsigned int
3697 GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
3698  const bool,
3699  const bool,
3700  const bool)
3701 {
3702  Assert(dim > 1, ExcImpossibleInDim(dim));
3703  Assert(vertex < GeometryInfo<dim>::vertices_per_face,
3704  ExcIndexRange(vertex, 0, GeometryInfo<dim>::vertices_per_face));
3705  return vertex;
3706 }
3707 
3708 
3709 
3710 template <>
3711 inline unsigned int
3712 GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
3713  const bool face_orientation,
3714  const bool face_flip,
3715  const bool face_rotation)
3716 {
3717  Assert(line < GeometryInfo<3>::lines_per_face,
3718  ExcIndexRange(line, 0, GeometryInfo<3>::lines_per_face));
3719 
3720 
3721  // make sure we handle
3722  // non-standard faces correctly
3723  //
3724  // so set up a table that for each line (of a
3725  // quad) describes which line to take
3726  //
3727  // first index: four lines 0...3
3728  //
3729  // second index: face_orientation; 0:
3730  // opposite normal, 1: standard
3731  //
3732  // third index: face_flip; 0: standard, 1:
3733  // face rotated by 180 degrees
3734  //
3735  // forth index: face_rotation: 0: standard,
3736  // 1: face rotated by 90 degrees
3737 
3738  const unsigned int line_translation[4][2][2][2] = {
3739  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
3740  // face_rotation=false and true
3741  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
3742  // face_rotation=false and true
3743  {{0, 3}, // line 0, face_orientation=true, face_flip=false,
3744  // face_rotation=false and true
3745  {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
3746  // face_rotation=false and true
3747 
3748  {{{3, 1}, // line 1 ...
3749  {2, 0}},
3750  {{1, 2}, {0, 3}}},
3751 
3752  {{{0, 3}, // line 2 ...
3753  {1, 2}},
3754  {{2, 0}, {3, 1}}},
3755 
3756  {{{1, 2}, // line 3 ...
3757  {0, 3}},
3758  {{3, 1}, {2, 0}}}};
3759 
3760  return line_translation[line][face_orientation][face_flip][face_rotation];
3761 }
3762 
3763 
3764 
3765 template <int dim>
3766 inline unsigned int
3767 GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
3768  const bool,
3769  const bool,
3770  const bool)
3771 {
3772  Assert(false, ExcNotImplemented());
3773  return line;
3774 }
3775 
3776 
3777 
3778 template <>
3779 inline unsigned int
3780 GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
3781  const bool face_orientation,
3782  const bool face_flip,
3783  const bool face_rotation)
3784 {
3785  Assert(line < GeometryInfo<3>::lines_per_face,
3786  ExcIndexRange(line, 0, GeometryInfo<3>::lines_per_face));
3787 
3788 
3789  // make sure we handle
3790  // non-standard faces correctly
3791  //
3792  // so set up a table that for each line (of a
3793  // quad) describes which line to take
3794  //
3795  // first index: four lines 0...3
3796  //
3797  // second index: face_orientation; 0:
3798  // opposite normal, 1: standard
3799  //
3800  // third index: face_flip; 0: standard, 1:
3801  // face rotated by 180 degrees
3802  //
3803  // forth index: face_rotation: 0: standard,
3804  // 1: face rotated by 90 degrees
3805 
3806  const unsigned int line_translation[4][2][2][2] = {
3807  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
3808  // face_rotation=false and true
3809  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
3810  // face_rotation=false and true
3811  {{0, 2}, // line 0, face_orientation=true, face_flip=false,
3812  // face_rotation=false and true
3813  {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
3814  // face_rotation=false and true
3815 
3816  {{{3, 1}, // line 1 ...
3817  {2, 0}},
3818  {{1, 3}, {0, 2}}},
3819 
3820  {{{0, 3}, // line 2 ...
3821  {1, 2}},
3822  {{2, 1}, {3, 0}}},
3823 
3824  {{{1, 2}, // line 3 ...
3825  {0, 3}},
3826  {{3, 0}, {2, 1}}}};
3827 
3828  return line_translation[line][face_orientation][face_flip][face_rotation];
3829 }
3830 
3831 
3832 
3833 template <int dim>
3834 inline unsigned int
3835 GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
3836  const bool,
3837  const bool,
3838  const bool)
3839 {
3840  Assert(false, ExcNotImplemented());
3841  return line;
3842 }
3843 
3844 
3845 
3846 template <>
3847 inline unsigned int
3848 GeometryInfo<1>::child_cell_on_face(const RefinementCase<1> &,
3849  const unsigned int face,
3850  const unsigned int subface,
3851  const bool,
3852  const bool,
3853  const bool,
3854  const RefinementCase<0> &)
3855 {
3856  (void)subface;
3857  Assert(face < faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
3858  Assert(subface < max_children_per_face,
3859  ExcIndexRange(subface, 0, max_children_per_face));
3860 
3861  return face;
3862 }
3863 
3864 
3865 
3866 template <>
3867 inline unsigned int
3868 GeometryInfo<2>::child_cell_on_face(const RefinementCase<2> &ref_case,
3869  const unsigned int face,
3870  const unsigned int subface,
3871  const bool /*face_orientation*/,
3872  const bool face_flip,
3873  const bool /*face_rotation*/,
3874  const RefinementCase<1> &)
3875 {
3876  Assert(face < faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
3877  Assert(subface < max_children_per_face,
3878  ExcIndexRange(subface, 0, max_children_per_face));
3879 
3880  // always return the child adjacent to the specified
3881  // subface. if the face of a cell is not refined, don't
3882  // throw an assertion but deliver the child adjacent to
3883  // the face nevertheless, i.e. deliver the child of
3884  // this cell adjacent to the subface of a possibly
3885  // refined neighbor. this simplifies setting neighbor
3886  // information in execute_refinement.
3887  const unsigned int
3888  subcells[2][RefinementCase<2>::isotropic_refinement][faces_per_cell]
3889  [max_children_per_face] = {
3890  {
3891  // Normal orientation (face_flip = false)
3892  {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
3893  {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
3894  {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
3895  },
3896  {
3897  // Flipped orientation (face_flip = true)
3898  {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
3899  {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
3900  {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
3901  }};
3902 
3903  return subcells[face_flip][ref_case - 1][face][subface];
3904 }
3905 
3906 
3907 
3908 template <>
3909 inline unsigned int
3910 GeometryInfo<3>::child_cell_on_face(const RefinementCase<3> &ref_case,
3911  const unsigned int face,
3912  const unsigned int subface,
3913  const bool face_orientation,
3914  const bool face_flip,
3915  const bool face_rotation,
3916  const RefinementCase<2> &face_ref_case)
3917 {
3918  const unsigned int dim = 3;
3919 
3920  Assert(ref_case > RefinementCase<dim - 1>::no_refinement,
3921  ExcMessage("Cell has no children."));
3922  Assert(face < faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
3923  Assert(subface < GeometryInfo<dim - 1>::n_children(face_ref_case) ||
3924  (subface == 0 &&
3925  face_ref_case == RefinementCase<dim - 1>::no_refinement),
3926  ExcIndexRange(subface, 0, GeometryInfo<2>::n_children(face_ref_case)));
3927 
3928  // invalid number used for invalid cases,
3929  // e.g. when the children are more refined at
3930  // a given face than the face itself
3931  const unsigned int e = numbers::invalid_unsigned_int;
3932 
3933  // the whole process of finding a child cell
3934  // at a given subface considering the
3935  // possibly anisotropic refinement cases of
3936  // the cell and the face as well as
3937  // orientation, flip and rotation of the face
3938  // is quite complicated. thus, we break it
3939  // down into several steps.
3940 
3941  // first step: convert the given face refine
3942  // case to a face refine case concerning the
3943  // face in standard orientation (, flip and
3944  // rotation). This only affects cut_x and
3945  // cut_y
3946  const RefinementCase<dim - 1> flip[4] = {
3947  RefinementCase<dim - 1>::no_refinement,
3948  RefinementCase<dim - 1>::cut_y,
3949  RefinementCase<dim - 1>::cut_x,
3950  RefinementCase<dim - 1>::cut_xy};
3951  // for face_orientation, 'true' is the
3952  // default value whereas for face_rotation,
3953  // 'false' is standard. If
3954  // <tt>face_rotation==face_orientation</tt>,
3955  // then one of them is non-standard and we
3956  // have to swap cut_x and cut_y, otherwise no
3957  // change is necessary.
3958  const RefinementCase<dim - 1> std_face_ref =
3959  (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
3960 
3961  // second step: convert the given subface
3962  // index to the one for a standard face
3963  // respecting face_orientation, face_flip and
3964  // face_rotation
3965 
3966  // first index: face_ref_case
3967  // second index: face_orientation
3968  // third index: face_flip
3969  // forth index: face_rotation
3970  // fifth index: subface index
3971  const unsigned int subface_exchange[4][2][2][2][4] = {
3972  // no_refinement (subface 0 stays 0,
3973  // all others are invalid)
3974  {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
3975  {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
3976  // cut_x (here, if the face is only
3977  // rotated OR only falsely oriented,
3978  // then subface 0 of the non-standard
3979  // face does NOT correspond to one of
3980  // the subfaces of a standard
3981  // face. Thus we indicate the subface
3982  // which is located at the lower left
3983  // corner (the origin of the face's
3984  // local coordinate system) with
3985  // '0'. The rest of this issue is
3986  // taken care of using the above
3987  // conversion to a 'standard face
3988  // refine case')
3989  {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
3990  {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
3991  // cut_y (the same applies as for
3992  // cut_x)
3993  {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
3994  {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
3995  // cut_xyz: this information is
3996  // identical to the information
3997  // returned by
3998  // GeometryInfo<3>::real_to_standard_face_vertex()
3999  {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4000  // face_rotation=false, subfaces 0,1,2,3
4001  {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4002  // face_rotation=true, subfaces 0,1,2,3
4003  {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4004  // face_rotation=false, subfaces 0,1,2,3
4005  {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4006  // face_rotation=true, subfaces 0,1,2,3
4007  {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4008  // face_rotation=false, subfaces 0,1,2,3
4009  {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4010  // face_rotation=true, subfaces 0,1,2,3
4011  {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4012  // face_rotation=false, subfaces 0,1,2,3
4013  {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4014  // face_rotation=true, subfaces 0,1,2,3
4015 
4016  const unsigned int std_subface =
4017  subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4018  [subface];
4019  Assert(std_subface != e, ExcInternalError());
4020 
4021  // third step: these are the children, which
4022  // can be found at the given subfaces of an
4023  // isotropically refined (standard) face
4024  //
4025  // first index: (refinement_case-1)
4026  // second index: face_index
4027  // third index: subface_index (isotropic refinement)
4028  const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4029  [max_children_per_face] = {
4030  // cut_x
4031  {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4032  {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4033  {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4034  {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4035  {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4036  {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4037  // cut_y
4038  {{0, 1, 0, 1},
4039  {0, 1, 0, 1},
4040  {0, 0, 0, 0},
4041  {1, 1, 1, 1},
4042  {0, 0, 1, 1},
4043  {0, 0, 1, 1}},
4044  // cut_xy
4045  {{0, 2, 0, 2},
4046  {1, 3, 1, 3},
4047  {0, 0, 1, 1},
4048  {2, 2, 3, 3},
4049  {0, 1, 2, 3},
4050  {0, 1, 2, 3}},
4051  // cut_z
4052  {{0, 0, 1, 1},
4053  {0, 0, 1, 1},
4054  {0, 1, 0, 1},
4055  {0, 1, 0, 1},
4056  {0, 0, 0, 0},
4057  {1, 1, 1, 1}},
4058  // cut_xz
4059  {{0, 0, 1, 1},
4060  {2, 2, 3, 3},
4061  {0, 1, 2, 3},
4062  {0, 1, 2, 3},
4063  {0, 2, 0, 2},
4064  {1, 3, 1, 3}},
4065  // cut_yz
4066  {{0, 1, 2, 3},
4067  {0, 1, 2, 3},
4068  {0, 2, 0, 2},
4069  {1, 3, 1, 3},
4070  {0, 0, 1, 1},
4071  {2, 2, 3, 3}},
4072  // cut_xyz
4073  {{0, 2, 4, 6},
4074  {1, 3, 5, 7},
4075  {0, 4, 1, 5},
4076  {2, 6, 3, 7},
4077  {0, 1, 2, 3},
4078  {4, 5, 6, 7}}};
4079 
4080  // forth step: check, whether the given face
4081  // refine case is valid for the given cell
4082  // refine case. this is the case, if the
4083  // given face refine case is at least as
4084  // refined as the face is for the given cell
4085  // refine case
4086 
4087  // note, that we are considering standard
4088  // face refinement cases here and thus must
4089  // not pass the given orientation, flip and
4090  // rotation flags
4091  if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4092  face_refinement_case(ref_case, face))
4093  {
4094  // all is fine. for anisotropic face
4095  // refine cases, select one of the
4096  // isotropic subfaces which neighbors the
4097  // same child
4098 
4099  // first index: (standard) face refine case
4100  // second index: subface index
4101  const unsigned int equivalent_iso_subface[4][4] = {
4102  {0, e, e, e}, // no_refinement
4103  {0, 3, e, e}, // cut_x
4104  {0, 3, e, e}, // cut_y
4105  {0, 1, 2, 3}}; // cut_xy
4106 
4107  const unsigned int equ_std_subface =
4108  equivalent_iso_subface[std_face_ref][std_subface];
4109  Assert(equ_std_subface != e, ExcInternalError());
4110 
4111  return iso_children[ref_case - 1][face][equ_std_subface];
4112  }
4113  else
4114  {
4115  // the face_ref_case was too coarse,
4116  // throw an error
4117  Assert(false,
4118  ExcMessage("The face RefineCase is too coarse "
4119  "for the given cell RefineCase."));
4120  }
4121  // we only get here in case of an error
4122  return e;
4123 }
4124 
4125 
4126 
4127 template <>
4128 inline unsigned int
4129 GeometryInfo<4>::child_cell_on_face(const RefinementCase<4> &,
4130  const unsigned int,
4131  const unsigned int,
4132  const bool,
4133  const bool,
4134  const bool,
4135  const RefinementCase<3> &)
4136 {
4137  Assert(false, ExcNotImplemented());
4138  return numbers::invalid_unsigned_int;
4139 }
4140 
4141 
4142 
4143 template <>
4144 inline unsigned int
4145 GeometryInfo<1>::face_to_cell_lines(const unsigned int face,
4146  const unsigned int line,
4147  const bool,
4148  const bool,
4149  const bool)
4150 {
4151  (void)face;
4152  (void)line;
4153  Assert(face + 1 < faces_per_cell + 1, ExcIndexRange(face, 0, faces_per_cell));
4154  Assert(line + 1 < lines_per_face + 1, ExcIndexRange(line, 0, lines_per_face));
4155 
4156  // There is only a single line, so
4157  // it must be this.
4158  return 0;
4159 }
4160 
4161 
4162 
4163 template <>
4164 inline unsigned int
4165 GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4166  const unsigned int line,
4167  const bool,
4168  const bool,
4169  const bool)
4170 {
4171  (void)line;
4172  Assert(face < faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
4173  Assert(line < lines_per_face, ExcIndexRange(line, 0, lines_per_face));
4174 
4175  // The face is a line itself.
4176  return face;
4177 }
4178 
4179 
4180 
4181 template <>
4182 inline unsigned int
4183 GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4184  const unsigned int line,
4185  const bool face_orientation,
4186  const bool face_flip,
4187  const bool face_rotation)
4188 {
4189  Assert(face < faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
4190  Assert(line < lines_per_face, ExcIndexRange(line, 0, lines_per_face));
4191 
4192  const unsigned lines[faces_per_cell][lines_per_face] = {
4193  {8, 10, 0, 4}, // left face
4194  {9, 11, 1, 5}, // right face
4195  {2, 6, 8, 9}, // front face
4196  {3, 7, 10, 11}, // back face
4197  {0, 1, 2, 3}, // bottom face
4198  {4, 5, 6, 7}}; // top face
4199  return lines[face][real_to_standard_face_line(
4200  line, face_orientation, face_flip, face_rotation)];
4201 }
4202 
4203 
4204 
4205 template <int dim>
4206 inline unsigned int
4207 GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4208  const unsigned int,
4209  const bool,
4210  const bool,
4211  const bool)
4212 {
4213  Assert(false, ExcNotImplemented());
4214  return numbers::invalid_unsigned_int;
4215 }
4216 
4217 
4218 
4219 template <int dim>
4220 inline unsigned int
4221 GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4222  const unsigned int vertex,
4223  const bool face_orientation,
4224  const bool face_flip,
4225  const bool face_rotation)
4226 {
4227  return child_cell_on_face(RefinementCase<dim>::isotropic_refinement,
4228  face,
4229  vertex,
4230  face_orientation,
4231  face_flip,
4232  face_rotation);
4233 }
4234 
4235 
4236 
4237 template <int dim>
4238 inline Point<dim>
4239 GeometryInfo<dim>::project_to_unit_cell(const Point<dim> &q)
4240 {
4241  Point<dim> p = q;
4242  for (unsigned int i = 0; i < dim; i++)
4243  if (p[i] < 0.)
4244  p[i] = 0.;
4245  else if (p[i] > 1.)
4246  p[i] = 1.;
4247 
4248  return p;
4249 }
4250 
4251 
4252 
4253 template <int dim>
4254 inline double
4255 GeometryInfo<dim>::distance_to_unit_cell(const Point<dim> &p)
4256 {
4257  double result = 0.0;
4258 
4259  for (unsigned int i = 0; i < dim; i++)
4260  if ((-p[i]) > result)
4261  result = -p[i];
4262  else if ((p[i] - 1.) > result)
4263  result = (p[i] - 1.);
4264 
4265  return result;
4266 }
4267 
4268 
4269 
4270 template <int dim>
4271 inline double
4272 GeometryInfo<dim>::d_linear_shape_function(const Point<dim> & xi,
4273  const unsigned int i)
4274 {
4275  Assert(i < GeometryInfo<dim>::vertices_per_cell,
4276  ExcIndexRange(i, 0, GeometryInfo<dim>::vertices_per_cell));
4277 
4278  switch (dim)
4279  {
4280  case 1:
4281  {
4282  const double x = xi[0];
4283  switch (i)
4284  {
4285  case 0:
4286  return 1 - x;
4287  case 1:
4288  return x;
4289  }
4290  break;
4291  }
4292 
4293  case 2:
4294  {
4295  const double x = xi[0];
4296  const double y = xi[1];
4297  switch (i)
4298  {
4299  case 0:
4300  return (1 - x) * (1 - y);
4301  case 1:
4302  return x * (1 - y);
4303  case 2:
4304  return (1 - x) * y;
4305  case 3:
4306  return x * y;
4307  }
4308  break;
4309  }
4310 
4311  case 3:
4312  {
4313  const double x = xi[0];
4314  const double y = xi[1];
4315  const double z = xi[2];
4316  switch (i)
4317  {
4318  case 0:
4319  return (1 - x) * (1 - y) * (1 - z);
4320  case 1:
4321  return x * (1 - y) * (1 - z);
4322  case 2:
4323  return (1 - x) * y * (1 - z);
4324  case 3:
4325  return x * y * (1 - z);
4326  case 4:
4327  return (1 - x) * (1 - y) * z;
4328  case 5:
4329  return x * (1 - y) * z;
4330  case 6:
4331  return (1 - x) * y * z;
4332  case 7:
4333  return x * y * z;
4334  }
4335  break;
4336  }
4337 
4338  default:
4339  Assert(false, ExcNotImplemented());
4340  }
4341  return -1e9;
4342 }
4343 
4344 
4345 
4346 template <>
4347 Tensor<1, 1> inline GeometryInfo<1>::d_linear_shape_function_gradient(
4348  const Point<1> &,
4349  const unsigned int i)
4350 {
4351  Assert(i < GeometryInfo<1>::vertices_per_cell,
4352  ExcIndexRange(i, 0, GeometryInfo<1>::vertices_per_cell));
4353 
4354  switch (i)
4355  {
4356  case 0:
4357  return Point<1>(-1.);
4358  case 1:
4359  return Point<1>(1.);
4360  }
4361 
4362  return Point<1>(-1e9);
4363 }
4364 
4365 
4366 
4367 template <>
4368 Tensor<1, 2> inline GeometryInfo<2>::d_linear_shape_function_gradient(
4369  const Point<2> & xi,
4370  const unsigned int i)
4371 {
4372  Assert(i < GeometryInfo<2>::vertices_per_cell,
4373  ExcIndexRange(i, 0, GeometryInfo<2>::vertices_per_cell));
4374 
4375  const double x = xi[0];
4376  const double y = xi[1];
4377  switch (i)
4378  {
4379  case 0:
4380  return Point<2>(-(1 - y), -(1 - x));
4381  case 1:
4382  return Point<2>(1 - y, -x);
4383  case 2:
4384  return Point<2>(-y, 1 - x);
4385  case 3:
4386  return Point<2>(y, x);
4387  }
4388  return Point<2>(-1e9, -1e9);
4389 }
4390 
4391 
4392 
4393 template <>
4394 Tensor<1, 3> inline GeometryInfo<3>::d_linear_shape_function_gradient(
4395  const Point<3> & xi,
4396  const unsigned int i)
4397 {
4398  Assert(i < GeometryInfo<3>::vertices_per_cell,
4399  ExcIndexRange(i, 0, GeometryInfo<3>::vertices_per_cell));
4400 
4401  const double x = xi[0];
4402  const double y = xi[1];
4403  const double z = xi[2];
4404  switch (i)
4405  {
4406  case 0:
4407  return Point<3>(-(1 - y) * (1 - z),
4408  -(1 - x) * (1 - z),
4409  -(1 - x) * (1 - y));
4410  case 1:
4411  return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4412  case 2:
4413  return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4414  case 3:
4415  return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4416  case 4:
4417  return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4418  case 5:
4419  return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4420  case 6:
4421  return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4422  case 7:
4423  return Point<3>(y * z, x * z, x * y);
4424  }
4425 
4426  return Point<3>(-1e9, -1e9, -1e9);
4427 }
4428 
4429 
4430 
4431 template <int dim>
4432 inline Tensor<1, dim>
4433 GeometryInfo<dim>::d_linear_shape_function_gradient(const Point<dim> &,
4434  const unsigned int)
4435 {
4436  Assert(false, ExcNotImplemented());
4437  return Tensor<1, dim>();
4438 }
4439 
4440 
4441 unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4442 {
4443  return 0;
4444 }
4445 
4446 
4447 namespace internal
4448 {
4449  namespace GeometryInfoHelper
4450  {
4451  // wedge product of a single
4452  // vector in 2d: we just have to
4453  // rotate it by 90 degrees to the
4454  // right
4455  inline Tensor<1, 2>
4456  wedge_product(const Tensor<1, 2> (&derivative)[1])
4457  {
4458  Tensor<1, 2> result;
4459  result[0] = derivative[0][1];
4460  result[1] = -derivative[0][0];
4461 
4462  return result;
4463  }
4464 
4465 
4466  // wedge product of 2 vectors in
4467  // 3d is the cross product
4468  inline Tensor<1, 3>
4469  wedge_product(const Tensor<1, 3> (&derivative)[2])
4470  {
4471  return cross_product_3d(derivative[0], derivative[1]);
4472  }
4473 
4474 
4475  // wedge product of dim vectors
4476  // in dim-d: that's the
4477  // determinant of the matrix
4478  template <int dim>
4479  inline Tensor<0, dim>
4480  wedge_product(const Tensor<1, dim> (&derivative)[dim])
4481  {
4482  Tensor<2, dim> jacobian;
4483  for (unsigned int i = 0; i < dim; ++i)
4484  jacobian[i] = derivative[i];
4485 
4486  return determinant(jacobian);
4487  }
4488  } // namespace GeometryInfoHelper
4489 } // namespace internal
4490 
4491 
4492 template <int dim>
4493 template <int spacedim>
4494 inline void
4495 GeometryInfo<dim>::alternating_form_at_vertices
4496 # ifndef DEAL_II_CONSTEXPR_BUG
4497  (const Point<spacedim> (&vertices)[vertices_per_cell],
4498  Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell])
4499 # else
4500  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms)
4501 # endif
4502 {
4503  // for each of the vertices,
4504  // compute the alternating form
4505  // of the mapped unit
4506  // vectors. consider for
4507  // example the case of a quad
4508  // in spacedim==3: to do so, we
4509  // need to see how the
4510  // infinitesimal vectors
4511  // (d\xi_1,0) and (0,d\xi_2)
4512  // are transformed into
4513  // spacedim-dimensional space
4514  // and then form their cross
4515  // product (i.e. the wedge product
4516  // of two vectors). to this end, note
4517  // that
4518  // \vec x = sum_i \vec v_i phi_i(\vec xi)
4519  // so the transformed vectors are
4520  // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
4521  // and
4522  // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
4523  // which boils down to the columns
4524  // of the 3x2 matrix \grad_\xi x(\xi)
4525  //
4526  // a similar reasoning would
4527  // hold for all dim,spacedim
4528  // pairs -- we only have to
4529  // compute the wedge product of
4530  // the columns of the
4531  // derivatives
4532  for (unsigned int i = 0; i < vertices_per_cell; ++i)
4533  {
4534  Tensor<1, spacedim> derivatives[dim];
4535 
4536  for (unsigned int j = 0; j < vertices_per_cell; ++j)
4537  {
4538  const Tensor<1, dim> grad_phi_j =
4539  d_linear_shape_function_gradient(unit_cell_vertex(i), j);
4540  for (unsigned int l = 0; l < dim; ++l)
4541  derivatives[l] += vertices[j] * grad_phi_j[l];
4542  }
4543 
4544  forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
4545  }
4546 }
4547 
4548 #endif // DOXYGEN
4549 DEAL_II_NAMESPACE_CLOSE
4550 
4551 #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:2688
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:1727
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)
RefinementCase::operator &
RefinementCase operator &(const RefinementCase &r) const
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:1159
GeometryInfo::vertices_per_cell
static constexpr unsigned int vertices_per_cell
Definition: geometry_info.h:1733
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)
RefinementCase::operator|
RefinementCase operator|(const RefinementCase &r) const
operator &
UpdateFlags operator &(const UpdateFlags f1, const UpdateFlags f2)
Definition: fe_update_flags.h:339
DeclException1
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:518
GeometryInfo::lines_per_cell
static constexpr unsigned int lines_per_cell
Definition: geometry_info.h:1762
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1407
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:1718
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)
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:1174
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:1744
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)
GeometryInfo::n_subfaces
static unsigned int n_subfaces(const internal::SubfaceCase< dim > &subface_case)
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)
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
RefinementCase::value
std::uint8_t value
Definition: geometry_info.h:716
DeclException3
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:564
GeometryInfo::subface_ratio
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
RefinementCase::cut_axis
static RefinementCase cut_axis(const unsigned int i)
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
Generated by   doxygen 1.8.13