Reference documentation for deal.II version 9.0.0
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 at
12 // the top level of the deal.II distribution.
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 #include <deal.II/base/exceptions.h>
22 #include <deal.II/base/point.h>
23 
24 #include <cstdint>
25 
26 
27 
28 DEAL_II_NAMESPACE_OPEN
29 
30 #ifndef DOXYGEN
31 namespace internal
32 {
33  namespace GeometryInfoHelper
34  {
35  // A struct that holds the values for all the arrays we want to initialize in GeometryInfo
36  template <int dim>
37  struct Initializers;
38 
39  template<>
40  struct Initializers<1>
41  {
42  static constexpr std::array<unsigned int, 2> ucd_to_deal
43  {{0, 1}};
44 
45  static constexpr std::array<unsigned int, 2> unit_normal_direction
46  {{ 0, 0 }};
47 
48  static constexpr std::array<int, 2> unit_normal_orientation
49  {{ -1, 1 }};
50 
51  static constexpr std::array<unsigned int, 2> opposite_face
52  {{ 1, 0 }};
53 
54  static constexpr std::array<unsigned int, 2> dx_to_deal
55  {{ 0, 1}};
56 
57  static constexpr std::array<std::array<unsigned int, 1>, 2> vertex_to_face
58  {
59  { {{ 0 }},
60  {{ 1 }}
61  }
62  };
63  };
64 
65  template<>
66  struct Initializers<2>
67  {
68  static constexpr std::array<unsigned int, 4> ucd_to_deal
69  {{ 0, 1, 3, 2}};
70 
71  static constexpr std::array<unsigned int, 4> unit_normal_direction
72  {{ 0, 0, 1, 1 }};
73 
74  static constexpr std::array<int, 4> unit_normal_orientation
75  {{-1, 1, -1, 1}};
76 
77  static constexpr std::array<unsigned int, 4> opposite_face
78  {{1,0,3,2}};
79 
80  static constexpr std::array<unsigned int, 4> dx_to_deal
81  {{ 0, 2, 1, 3}};
82 
83  static constexpr std::array<std::array<unsigned int, 2>, 4> vertex_to_face
84  {
85  { {{ 0, 2 }},
86  {{ 1, 2 }},
87  {{ 0, 3 }},
88  {{ 1, 3 }}
89  }
90  };
91  };
92 
93  template<>
94  struct Initializers<3>
95  {
96  static constexpr std::array<unsigned int, 8> ucd_to_deal
97  {{ 0, 1, 5, 4, 2, 3, 7, 6}};
98 
99  static constexpr std::array<unsigned int, 6> unit_normal_direction
100  {{ 0, 0, 1, 1, 2, 2 }};
101 
102  static constexpr std::array<int, 6> unit_normal_orientation
103  {{ -1, 1, -1, 1, -1, 1 }};
104 
105  static constexpr std::array<unsigned int, 6> opposite_face
106  {{ 1, 0, 3, 2, 5, 4 }};
107 
108  static constexpr std::array<unsigned int, 8> dx_to_deal
109  {{ 0, 4, 2, 6, 1, 5, 3, 7}};
110 
111  static constexpr std::array<std::array<unsigned int, 3>, 8> vertex_to_face
112  {
113  { {{ 0, 2, 4 }},
114  {{ 1, 2, 4 }},
115  {{ 0, 3, 4 }},
116  {{ 1, 3, 4 }},
117  {{ 0, 2, 5 }},
118  {{ 1, 2, 5 }},
119  {{ 0, 3, 5 }},
120  {{ 1, 3, 5 }}
121  }
122  };
123  };
124 
125  template<>
126  struct Initializers<4>
127  {
128  static constexpr std::array<unsigned int, 16> ucd_to_deal
129  {
130  {
131  numbers::invalid_unsigned_int,
132  numbers::invalid_unsigned_int,
133  numbers::invalid_unsigned_int,
134  numbers::invalid_unsigned_int,
135  numbers::invalid_unsigned_int,
136  numbers::invalid_unsigned_int,
137  numbers::invalid_unsigned_int,
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  }
148  };
149 
150  static constexpr std::array<unsigned int, 8> unit_normal_direction
151  {{ 0, 0, 1, 1, 2, 2, 3, 3 }};
152 
153  static constexpr std::array<int, 8> unit_normal_orientation
154  {{ -1, 1, -1, 1, -1, 1, -1, 1 }};
155 
156  static constexpr std::array<unsigned int, 8> opposite_face
157  {{ 1, 0, 3, 2, 5, 4, 7, 6 }};
158 
159  static constexpr std::array<unsigned int, 16> dx_to_deal
160  {
161  {
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  }
179  };
180 
181  static constexpr std::array<std::array<unsigned int, 4>, 16> vertex_to_face
182  {
183  { {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}},
184  {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}},
185  {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}},
186  {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}},
187  {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}},
188  {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}},
189  {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}},
190  {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}},
191  {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}},
192  {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}},
193  {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}},
194  {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}},
195  {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}},
196  {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}},
197  {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}},
198  {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}}
199  }
200  };
201  };
202  } // GeometryInfoHelper
203 } // internal
204 #endif //DOXYGEN
205 
206 
224 class GeometryPrimitive
225 {
226 public:
233  enum Object
234  {
238  vertex = 0,
242  line = 1,
246  quad = 2,
250  hex = 3
251  };
252 
257  GeometryPrimitive (const Object object);
258 
264  GeometryPrimitive (const unsigned int object_dimension);
265 
270  operator unsigned int () const;
271 
272 private:
276  Object object;
277 };
278 
279 
280 
281 
295 template <int dim>
296 struct RefinementPossibilities
297 {
332  enum Possibilities
333  {
337  no_refinement = 0,
338 
342  isotropic_refinement = static_cast<std::uint8_t>(-1)
343  };
344 };
345 
346 
347 
358 template <>
359 struct RefinementPossibilities<1>
360 {
395  enum Possibilities
396  {
400  no_refinement = 0,
404  cut_x = 1,
408  isotropic_refinement = cut_x
409  };
410 };
411 
412 
413 
425 template <>
426 struct RefinementPossibilities<2>
427 {
462  enum Possibilities
463  {
467  no_refinement = 0,
471  cut_x = 1,
475  cut_y = 2,
479  cut_xy = cut_x | cut_y,
480 
484  isotropic_refinement = cut_xy
485  };
486 };
487 
488 
489 
501 template <>
502 struct RefinementPossibilities<3>
503 {
538  enum Possibilities
539  {
543  no_refinement = 0,
547  cut_x = 1,
551  cut_y = 2,
555  cut_xy = cut_x | cut_y,
559  cut_z = 4,
563  cut_xz = cut_x | cut_z,
567  cut_yz = cut_y | cut_z,
571  cut_xyz = cut_x | cut_y | cut_z,
572 
576  isotropic_refinement = cut_xyz
577  };
578 };
579 
580 
581 
593 template <int dim>
594 class RefinementCase : public RefinementPossibilities<dim>
595 {
596 public:
600  RefinementCase ();
601 
606  RefinementCase (const typename RefinementPossibilities<dim>::Possibilities refinement_case);
607 
613  explicit RefinementCase (const std::uint8_t refinement_case);
614 
626  operator std::uint8_t () const;
627 
632  RefinementCase operator | (const RefinementCase &r) const;
633 
638  RefinementCase operator & (const RefinementCase &r) const;
639 
647  RefinementCase operator ~ () const;
648 
649 
655  static
656  RefinementCase cut_axis (const unsigned int i);
657 
661  static std::size_t memory_consumption ();
662 
667  template <class Archive>
668  void serialize(Archive &ar,
669  const unsigned int version);
670 
674  DeclException1 (ExcInvalidRefinementCase,
675  int,
676  << "The refinement flags given (" << arg1 << ") contain set bits that do not "
677  << "make sense for the space dimension of the object to which they are applied.");
678 
679 private:
684 std::uint8_t value :
685  (dim > 0 ? dim : 1);
686 };
687 
688 
689 namespace internal
690 {
691 
712  template <int dim>
713  struct SubfacePossibilities
714  {
718  enum Possibilities
719  {
723  case_none = 0,
724 
728  case_isotropic = static_cast<std::uint8_t>(-1)
729  };
730  };
731 
732 
742  template <>
743  struct SubfacePossibilities<0>
744  {
750  enum Possibilities
751  {
755  case_none = 0,
756 
760  case_isotropic = case_none
761  };
762  };
763 
764 
765 
776  template <>
777  struct SubfacePossibilities<1>
778  {
784  enum Possibilities
785  {
789  case_none = 0,
790 
794  case_isotropic = case_none
795  };
796  };
797 
798 
799 
811  template <>
812  struct SubfacePossibilities<2>
813  {
820  enum Possibilities
821  {
825  case_none = 0,
829  case_x = 1,
833  case_isotropic = case_x
834  };
835  };
836 
837 
838 
931  template <>
932  struct SubfacePossibilities<3>
933  {
940  enum Possibilities
941  {
942  case_none = 0,
943  case_x = 1,
944  case_x1y = 2,
945  case_x2y = 3,
946  case_x1y2y = 4,
947  case_y = 5,
948  case_y1x = 6,
949  case_y2x = 7,
950  case_y1x2x = 8,
951  case_xy = 9,
952 
953  case_isotropic = case_xy
954  };
955  };
956 
957 
958 
959 
967  template <int dim>
968  class SubfaceCase : public SubfacePossibilities<dim>
969  {
970  public:
976  SubfaceCase (const typename SubfacePossibilities<dim>::Possibilities subface_possibility);
977 
989  operator std::uint8_t () const;
990 
994  static constexpr std::size_t memory_consumption ();
995 
999  DeclException1 (ExcInvalidSubfaceCase,
1000  int,
1001  << "The subface case given (" << arg1 << ") does not make sense "
1002  << "for the space dimension of the object to which they are applied.");
1003 
1004  private:
1009  std::uint8_t value :
1010  (dim == 3 ? 4 : 1);
1011  };
1012 
1013 } // namespace internal
1014 
1015 
1016 
1017 template <int dim> struct GeometryInfo;
1018 
1019 
1020 
1021 
1039 template <>
1040 struct GeometryInfo<0>
1041 {
1042 
1050  static constexpr unsigned int max_children_per_cell = 1;
1051 
1055  static constexpr unsigned int faces_per_cell = 0;
1056 
1064  static constexpr unsigned int max_children_per_face = 0;
1065 
1071  static unsigned int n_children(const RefinementCase<0> &refinement_case);
1072 
1076  static constexpr unsigned int vertices_per_cell = 1;
1077 
1084  static constexpr unsigned int vertices_per_face = 0;
1085 
1089  static constexpr unsigned int lines_per_face = 0;
1090 
1094  static constexpr unsigned int quads_per_face = 0;
1095 
1099  static constexpr unsigned int lines_per_cell = 0;
1100 
1104  static constexpr unsigned int quads_per_cell = 0;
1105 
1109  static constexpr unsigned int hexes_per_cell = 0;
1110 
1128  static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1129 
1143  static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1144 };
1145 
1146 
1147 
1148 
1149 
1674 template <int dim>
1675 struct GeometryInfo
1676 {
1677 
1685  static constexpr unsigned int max_children_per_cell = 1 << dim;
1686 
1690  static constexpr unsigned int faces_per_cell = 2 * dim;
1691 
1699  static constexpr unsigned int max_children_per_face = GeometryInfo<dim-1>::max_children_per_cell;
1700 
1704  static constexpr unsigned int vertices_per_cell = 1 << dim;
1705 
1709  static constexpr unsigned int vertices_per_face = GeometryInfo<dim-1>::vertices_per_cell;
1710 
1714  static constexpr unsigned int lines_per_face
1715  = GeometryInfo<dim-1>::lines_per_cell;
1716 
1720  static constexpr unsigned int quads_per_face
1721  = GeometryInfo<dim-1>::quads_per_cell;
1722 
1732  static constexpr unsigned int lines_per_cell
1733  = (2*GeometryInfo<dim-1>::lines_per_cell +
1734  GeometryInfo<dim-1>::vertices_per_cell);
1735 
1743  static constexpr unsigned int quads_per_cell
1744  = (2*GeometryInfo<dim-1>::quads_per_cell +
1745  GeometryInfo<dim-1>::lines_per_cell);
1746 
1750  static constexpr unsigned int hexes_per_cell
1751  = (2*GeometryInfo<dim-1>::hexes_per_cell +
1752  GeometryInfo<dim-1>::quads_per_cell);
1753 
1771  static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal
1772  = internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal;
1773 
1787  static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal
1788  = internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal;
1789 
1800  static constexpr std::array<std::array<unsigned int, dim>, vertices_per_cell> vertex_to_face
1801  = internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face;
1802 
1807  static
1808  unsigned int
1809  n_children(const RefinementCase<dim> &refinement_case);
1810 
1815  static
1816  unsigned int
1817  n_subfaces(const internal::SubfaceCase<dim> &subface_case);
1818 
1828  static
1829  double
1830  subface_ratio(const internal::SubfaceCase<dim> &subface_case,
1831  const unsigned int subface_no);
1832 
1838  static
1839  RefinementCase<dim-1>
1840  face_refinement_case (const RefinementCase<dim> &cell_refinement_case,
1841  const unsigned int face_no,
1842  const bool face_orientation = true,
1843  const bool face_flip = false,
1844  const bool face_rotation = false);
1845 
1851  static
1852  RefinementCase<dim>
1853  min_cell_refinement_case_for_face_refinement
1854  (const RefinementCase<dim-1> &face_refinement_case,
1855  const unsigned int face_no,
1856  const bool face_orientation = true,
1857  const bool face_flip = false,
1858  const bool face_rotation = false);
1859 
1864  static
1865  RefinementCase<1>
1866  line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
1867  const unsigned int line_no);
1868 
1873  static
1874  RefinementCase<dim>
1875  min_cell_refinement_case_for_line_refinement(const unsigned int line_no);
1876 
1923  static
1924  unsigned int
1925  child_cell_on_face (const RefinementCase<dim> &ref_case,
1926  const unsigned int face,
1927  const unsigned int subface,
1928  const bool face_orientation = true,
1929  const bool face_flip = false,
1930  const bool face_rotation = false,
1931  const RefinementCase<dim-1> &face_refinement_case
1932  = RefinementCase<dim-1>::isotropic_refinement);
1933 
1947  static
1948  unsigned int
1949  line_to_cell_vertices (const unsigned int line,
1950  const unsigned int vertex);
1951 
1972  static
1973  unsigned int
1974  face_to_cell_vertices (const unsigned int face,
1975  const unsigned int vertex,
1976  const bool face_orientation = true,
1977  const bool face_flip = false,
1978  const bool face_rotation = false);
1979 
1991  static
1992  unsigned int
1993  face_to_cell_lines (const unsigned int face,
1994  const unsigned int line,
1995  const bool face_orientation = true,
1996  const bool face_flip = false,
1997  const bool face_rotation = false);
1998 
2008  static
2009  unsigned int
2010  standard_to_real_face_vertex (const unsigned int vertex,
2011  const bool face_orientation = true,
2012  const bool face_flip = false,
2013  const bool face_rotation = false);
2014 
2024  static
2025  unsigned int
2026  real_to_standard_face_vertex (const unsigned int vertex,
2027  const bool face_orientation = true,
2028  const bool face_flip = false,
2029  const bool face_rotation = false);
2030 
2040  static
2041  unsigned int
2042  standard_to_real_face_line (const unsigned int line,
2043  const bool face_orientation = true,
2044  const bool face_flip = false,
2045  const bool face_rotation = false);
2046 
2056  static
2057  unsigned int
2058  real_to_standard_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 
2068  static
2069  Point<dim>
2070  unit_cell_vertex (const unsigned int vertex);
2071 
2081  static
2082  unsigned int
2083  child_cell_from_point (const Point<dim> &p);
2084 
2092  static
2093  Point<dim>
2094  cell_to_child_coordinates (const Point<dim> &p,
2095  const unsigned int child_index,
2096  const RefinementCase<dim> refine_case
2097  = RefinementCase<dim>::isotropic_refinement);
2098 
2104  static
2105  Point<dim>
2106  child_to_cell_coordinates (const Point<dim> &p,
2107  const unsigned int child_index,
2108  const RefinementCase<dim> refine_case
2109  = RefinementCase<dim>::isotropic_refinement);
2110 
2115  static
2116  bool
2117  is_inside_unit_cell (const Point<dim> &p);
2118 
2130  static
2131  bool
2132  is_inside_unit_cell (const Point<dim> &p,
2133  const double eps);
2134 
2139  static
2140  Point<dim>
2141  project_to_unit_cell (const Point<dim> &p);
2142 
2148  static
2149  double
2150  distance_to_unit_cell (const Point<dim> &p);
2151 
2156  static
2157  double
2158  d_linear_shape_function (const Point<dim> &xi,
2159  const unsigned int i);
2160 
2165  static
2166  Tensor<1,dim>
2167  d_linear_shape_function_gradient (const Point<dim> &xi,
2168  const unsigned int i);
2169 
2221  template <int spacedim>
2222  static void
2223  alternating_form_at_vertices
2224 #ifndef DEAL_II_CONSTEXPR_BUG
2225  (const Point<spacedim> (&vertices)[vertices_per_cell],
2226  Tensor<spacedim-dim,spacedim> (&forms)[vertices_per_cell]);
2227 #else
2228  (const Point<spacedim> *vertices,
2229  Tensor<spacedim-dim,spacedim> *forms);
2230 #endif
2231 
2241  static constexpr std::array<unsigned int, faces_per_cell> unit_normal_direction
2242  = internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction;
2243 
2260  static constexpr std::array<int, faces_per_cell> unit_normal_orientation
2261  = internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation;
2262 
2268  static constexpr std::array<unsigned int, faces_per_cell> opposite_face
2269  = internal::GeometryInfoHelper::Initializers<dim>::opposite_face;
2270 
2271 
2275  DeclException1 (ExcInvalidCoordinate,
2276  double,
2277  << "The coordinates must satisfy 0 <= x_i <= 1, "
2278  << "but here we have x_i=" << arg1);
2279 
2283  DeclException3 (ExcInvalidSubface,
2284  int, int, int,
2285  << "RefinementCase<dim> " << arg1 << ": face " << arg2
2286  << " has no subface " << arg3);
2287 };
2288 
2289 
2290 
2291 
2292 #ifndef DOXYGEN
2293 
2294 
2295 /* -------------- declaration of explicit specializations ------------- */
2296 
2297 
2298 template <>
2299 Tensor<1,1>
2300 GeometryInfo<1>::
2301 d_linear_shape_function_gradient (const Point<1> &xi,
2302  const unsigned int i);
2303 template <>
2304 Tensor<1,2>
2305 GeometryInfo<2>::
2306 d_linear_shape_function_gradient (const Point<2> &xi,
2307  const unsigned int i);
2308 template <>
2309 Tensor<1,3>
2310 GeometryInfo<3>::
2311 d_linear_shape_function_gradient (const Point<3> &xi,
2312  const unsigned int i);
2313 
2314 
2315 
2316 
2317 /* -------------- inline functions ------------- */
2318 
2319 
2320 inline
2321 GeometryPrimitive::GeometryPrimitive (const Object object)
2322  :
2323  object (object)
2324 {}
2325 
2326 
2327 
2328 inline
2329 GeometryPrimitive::GeometryPrimitive (const unsigned int object_dimension)
2330  :
2331  object (static_cast<Object>(object_dimension))
2332 {}
2333 
2334 
2335 inline
2336 GeometryPrimitive::operator unsigned int () const
2337 {
2338  return static_cast<unsigned int>(object);
2339 }
2340 
2341 
2342 
2343 namespace internal
2344 {
2345 
2346  template <int dim>
2347  inline
2348  SubfaceCase<dim>::SubfaceCase (const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2349  :
2350  value (subface_possibility)
2351  {}
2352 
2353 
2354  template <int dim>
2355  inline
2356  SubfaceCase<dim>::operator std::uint8_t () const
2357  {
2358  return value;
2359  }
2360 
2361 
2362 } // namespace internal
2363 
2364 
2365 template <int dim>
2366 inline
2367 RefinementCase<dim>
2368 RefinementCase<dim>::cut_axis (const unsigned int)
2369 {
2370  Assert (false, ExcInternalError());
2371  return static_cast<std::uint8_t>(-1);
2372 }
2373 
2374 
2375 template <>
2376 inline
2377 RefinementCase<1>
2378 RefinementCase<1>::cut_axis (const unsigned int i)
2379 {
2380  Assert (i < 1, ExcIndexRange(i, 0, 1));
2381 
2382  const RefinementCase options[1] = { RefinementPossibilities<1>::cut_x };
2383  return options[i];
2384 }
2385 
2386 
2387 
2388 template <>
2389 inline
2390 RefinementCase<2>
2391 RefinementCase<2>::cut_axis (const unsigned int i)
2392 {
2393  Assert (i < 2, ExcIndexRange(i, 0, 2));
2394 
2395  const RefinementCase options[2] = { RefinementPossibilities<2>::cut_x,
2396  RefinementPossibilities<2>::cut_y
2397  };
2398  return options[i];
2399 }
2400 
2401 
2402 
2403 template <>
2404 inline
2405 RefinementCase<3>
2406 RefinementCase<3>::cut_axis (const unsigned int i)
2407 {
2408  Assert (i < 3, ExcIndexRange(i, 0, 3));
2409 
2410  const RefinementCase options[3] = { RefinementPossibilities<3>::cut_x,
2411  RefinementPossibilities<3>::cut_y,
2412  RefinementPossibilities<3>::cut_z
2413  };
2414  return options[i];
2415 }
2416 
2417 
2418 
2419 template <int dim>
2420 inline
2421 RefinementCase<dim>::RefinementCase ()
2422  :
2423  value (RefinementPossibilities<dim>::no_refinement)
2424 {}
2425 
2426 
2427 
2428 template <int dim>
2429 inline
2430 RefinementCase<dim>::
2431 RefinementCase (const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2432  :
2433  value (refinement_case)
2434 {
2435  // check that only those bits of
2436  // the given argument are set that
2437  // make sense for a given space
2438  // dimension
2439  Assert ((refinement_case & RefinementPossibilities<dim>::isotropic_refinement) ==
2440  refinement_case,
2441  ExcInvalidRefinementCase (refinement_case));
2442 }
2443 
2444 
2445 
2446 template <int dim>
2447 inline
2448 RefinementCase<dim>::RefinementCase (const std::uint8_t refinement_case)
2449  :
2450  value (refinement_case)
2451 {
2452  // check that only those bits of
2453  // the given argument are set that
2454  // make sense for a given space
2455  // dimension
2456  Assert ((refinement_case & RefinementPossibilities<dim>::isotropic_refinement) ==
2457  refinement_case,
2458  ExcInvalidRefinementCase (refinement_case));
2459 }
2460 
2461 
2462 
2463 template <int dim>
2464 inline
2465 RefinementCase<dim>::operator std::uint8_t () const
2466 {
2467  return value;
2468 }
2469 
2470 
2471 
2472 template <int dim>
2473 inline
2474 RefinementCase<dim>
2475 RefinementCase<dim>::operator | (const RefinementCase<dim> &r) const
2476 {
2477  return RefinementCase<dim>(static_cast<std::uint8_t> (value | r.value));
2478 }
2479 
2480 
2481 
2482 template <int dim>
2483 inline
2484 RefinementCase<dim>
2485 RefinementCase<dim>::operator & (const RefinementCase<dim> &r) const
2486 {
2487  return RefinementCase<dim>(static_cast<std::uint8_t> (value & r.value));
2488 }
2489 
2490 
2491 
2492 template <int dim>
2493 inline
2494 RefinementCase<dim>
2495 RefinementCase<dim>::operator ~ () const
2496 {
2497  return RefinementCase<dim>(static_cast<std::uint8_t> (
2498  (~value) & RefinementPossibilities<dim>::isotropic_refinement));
2499 }
2500 
2501 
2502 
2503 
2504 template <int dim>
2505 inline
2506 std::size_t
2507 RefinementCase<dim>::memory_consumption ()
2508 {
2509  return sizeof(RefinementCase<dim>);
2510 }
2511 
2512 
2513 
2514 template <int dim>
2515 template <class Archive>
2516 inline
2517 void RefinementCase<dim>::serialize (Archive &ar,
2518  const unsigned int)
2519 {
2520  // serialization can't deal with bitfields, so copy from/to a full sized
2521  // std::uint8_t
2522  std::uint8_t uchar_value = value;
2523  ar &uchar_value;
2524  value = uchar_value;
2525 }
2526 
2527 
2528 
2529 
2530 template <>
2531 inline
2532 Point<1>
2533 GeometryInfo<1>::unit_cell_vertex (const unsigned int vertex)
2534 {
2535  Assert (vertex < vertices_per_cell,
2536  ExcIndexRange (vertex, 0, vertices_per_cell));
2537 
2538  return Point<1>(static_cast<double>(vertex));
2539 }
2540 
2541 
2542 
2543 template <>
2544 inline
2545 Point<2>
2546 GeometryInfo<2>::unit_cell_vertex (const unsigned int vertex)
2547 {
2548  Assert (vertex < vertices_per_cell,
2549  ExcIndexRange (vertex, 0, vertices_per_cell));
2550 
2551  return Point<2>(vertex%2, vertex/2);
2552 }
2553 
2554 
2555 
2556 template <>
2557 inline
2558 Point<3>
2559 GeometryInfo<3>::unit_cell_vertex (const unsigned int vertex)
2560 {
2561  Assert (vertex < vertices_per_cell,
2562  ExcIndexRange (vertex, 0, vertices_per_cell));
2563 
2564  return Point<3>(vertex%2, vertex/2%2, vertex/4);
2565 }
2566 
2567 
2568 
2569 template <int dim>
2570 inline
2571 Point<dim>
2572 GeometryInfo<dim>::unit_cell_vertex (const unsigned int)
2573 {
2574  Assert(false, ExcNotImplemented());
2575 
2576  return Point<dim> ();
2577 }
2578 
2579 
2580 
2581 template <>
2582 inline
2583 unsigned int
2584 GeometryInfo<1>::child_cell_from_point (const Point<1> &p)
2585 {
2586  Assert ((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2587 
2588  return (p[0] <= 0.5 ? 0 : 1);
2589 }
2590 
2591 
2592 
2593 template <>
2594 inline
2595 unsigned int
2596 GeometryInfo<2>::child_cell_from_point (const Point<2> &p)
2597 {
2598  Assert ((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2599  Assert ((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2600 
2601  return (p[0] <= 0.5 ?
2602  (p[1] <= 0.5 ? 0 : 2) :
2603  (p[1] <= 0.5 ? 1 : 3));
2604 }
2605 
2606 
2607 
2608 template <>
2609 inline
2610 unsigned int
2611 GeometryInfo<3>::child_cell_from_point (const Point<3> &p)
2612 {
2613  Assert ((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2614  Assert ((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2615  Assert ((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
2616 
2617  return (p[0] <= 0.5 ?
2618  (p[1] <= 0.5 ?
2619  (p[2] <= 0.5 ? 0 : 4) :
2620  (p[2] <= 0.5 ? 2 : 6)) :
2621  (p[1] <= 0.5 ?
2622  (p[2] <= 0.5 ? 1 : 5) :
2623  (p[2] <= 0.5 ? 3 : 7)));
2624 }
2625 
2626 
2627 template <int dim>
2628 inline
2629 unsigned int
2630 GeometryInfo<dim>::child_cell_from_point (const Point<dim> &)
2631 {
2632  Assert(false, ExcNotImplemented());
2633 
2634  return 0;
2635 }
2636 
2637 
2638 
2639 template <>
2640 inline
2641 Point<1>
2642 GeometryInfo<1>::cell_to_child_coordinates (const Point<1> &p,
2643  const unsigned int child_index,
2644  const RefinementCase<1> refine_case)
2645 
2646 {
2647  Assert (child_index < 2,
2648  ExcIndexRange (child_index, 0, 2));
2649  Assert (refine_case==RefinementCase<1>::cut_x,
2650  ExcInternalError());
2651  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
2652 
2653  return Point<1>(p*2.0-unit_cell_vertex(child_index));
2654 }
2655 
2656 
2657 
2658 template <>
2659 inline
2660 Point<2>
2661 GeometryInfo<2>::cell_to_child_coordinates (const Point<2> &p,
2662  const unsigned int child_index,
2663  const RefinementCase<2> refine_case)
2664 
2665 {
2666  Assert (child_index < GeometryInfo<2>::n_children(refine_case),
2667  ExcIndexRange (child_index, 0, GeometryInfo<2>::n_children(refine_case)));
2668 
2669  Point<2> point=p;
2670  switch (refine_case)
2671  {
2672  case RefinementCase<2>::cut_x:
2673  point[0]*=2.0;
2674  if (child_index==1)
2675  point[0]-=1.0;
2676  break;
2677  case RefinementCase<2>::cut_y:
2678  point[1]*=2.0;
2679  if (child_index==1)
2680  point[1]-=1.0;
2681  break;
2682  case RefinementCase<2>::cut_xy:
2683  point*=2.0;
2684  point-=unit_cell_vertex(child_index);
2685  break;
2686  default:
2687  Assert(false, ExcInternalError());
2688  }
2689 
2690  return point;
2691 }
2692 
2693 
2694 
2695 template <>
2696 inline
2697 Point<3>
2698 GeometryInfo<3>::cell_to_child_coordinates (const Point<3> &p,
2699  const unsigned int child_index,
2700  const RefinementCase<3> refine_case)
2701 
2702 {
2703  Assert (child_index < GeometryInfo<3>::n_children(refine_case),
2704  ExcIndexRange (child_index, 0, GeometryInfo<3>::n_children(refine_case)));
2705 
2706  Point<3> point=p;
2707  // there might be a cleverer way to do
2708  // this, but since this function is called
2709  // in very few places for initialization
2710  // purposes only, I don't care at the
2711  // moment
2712  switch (refine_case)
2713  {
2714  case RefinementCase<3>::cut_x:
2715  point[0]*=2.0;
2716  if (child_index==1)
2717  point[0]-=1.0;
2718  break;
2719  case RefinementCase<3>::cut_y:
2720  point[1]*=2.0;
2721  if (child_index==1)
2722  point[1]-=1.0;
2723  break;
2724  case RefinementCase<3>::cut_z:
2725  point[2]*=2.0;
2726  if (child_index==1)
2727  point[2]-=1.0;
2728  break;
2729  case RefinementCase<3>::cut_xy:
2730  point[0]*=2.0;
2731  point[1]*=2.0;
2732  if (child_index%2==1)
2733  point[0]-=1.0;
2734  if (child_index/2==1)
2735  point[1]-=1.0;
2736  break;
2737  case RefinementCase<3>::cut_xz:
2738  // careful, this is slightly
2739  // different from xy and yz due to
2740  // different internal numbering of
2741  // children!
2742  point[0]*=2.0;
2743  point[2]*=2.0;
2744  if (child_index/2==1)
2745  point[0]-=1.0;
2746  if (child_index%2==1)
2747  point[2]-=1.0;
2748  break;
2749  case RefinementCase<3>::cut_yz:
2750  point[1]*=2.0;
2751  point[2]*=2.0;
2752  if (child_index%2==1)
2753  point[1]-=1.0;
2754  if (child_index/2==1)
2755  point[2]-=1.0;
2756  break;
2757  case RefinementCase<3>::cut_xyz:
2758  point*=2.0;
2759  point-=unit_cell_vertex(child_index);
2760  break;
2761  default:
2762  Assert(false, ExcInternalError());
2763  }
2764 
2765  return point;
2766 }
2767 
2768 
2769 
2770 template <int dim>
2771 inline
2772 Point<dim>
2773 GeometryInfo<dim>::cell_to_child_coordinates (const Point<dim> &/*p*/,
2774  const unsigned int /*child_index*/,
2775  const RefinementCase<dim> /*refine_case*/)
2776 
2777 {
2778  Assert (false, ExcNotImplemented());
2779  return Point<dim>();
2780 }
2781 
2782 
2783 
2784 template <>
2785 inline
2786 Point<1>
2787 GeometryInfo<1>::child_to_cell_coordinates (const Point<1> &p,
2788  const unsigned int child_index,
2789  const RefinementCase<1> refine_case)
2790 
2791 {
2792  Assert (child_index < 2,
2793  ExcIndexRange (child_index, 0, 2));
2794  Assert (refine_case==RefinementCase<1>::cut_x,
2795  ExcInternalError());
2796  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
2797 
2798  return (p+unit_cell_vertex(child_index))*0.5;
2799 }
2800 
2801 
2802 
2803 template <>
2804 inline
2805 Point<3>
2806 GeometryInfo<3>::child_to_cell_coordinates (const Point<3> &p,
2807  const unsigned int child_index,
2808  const RefinementCase<3> refine_case)
2809 
2810 {
2811  Assert (child_index < GeometryInfo<3>::n_children(refine_case),
2812  ExcIndexRange (child_index, 0, GeometryInfo<3>::n_children(refine_case)));
2813 
2814  Point<3> point=p;
2815  // there might be a cleverer way to do
2816  // this, but since this function is called
2817  // in very few places for initialization
2818  // purposes only, I don't care at the
2819  // moment
2820  switch (refine_case)
2821  {
2822  case RefinementCase<3>::cut_x:
2823  if (child_index==1)
2824  point[0]+=1.0;
2825  point[0]*=0.5;
2826  break;
2827  case RefinementCase<3>::cut_y:
2828  if (child_index==1)
2829  point[1]+=1.0;
2830  point[1]*=0.5;
2831  break;
2832  case RefinementCase<3>::cut_z:
2833  if (child_index==1)
2834  point[2]+=1.0;
2835  point[2]*=0.5;
2836  break;
2837  case RefinementCase<3>::cut_xy:
2838  if (child_index%2==1)
2839  point[0]+=1.0;
2840  if (child_index/2==1)
2841  point[1]+=1.0;
2842  point[0]*=0.5;
2843  point[1]*=0.5;
2844  break;
2845  case RefinementCase<3>::cut_xz:
2846  // careful, this is slightly
2847  // different from xy and yz due to
2848  // different internal numbering of
2849  // children!
2850  if (child_index/2==1)
2851  point[0]+=1.0;
2852  if (child_index%2==1)
2853  point[2]+=1.0;
2854  point[0]*=0.5;
2855  point[2]*=0.5;
2856  break;
2857  case RefinementCase<3>::cut_yz:
2858  if (child_index%2==1)
2859  point[1]+=1.0;
2860  if (child_index/2==1)
2861  point[2]+=1.0;
2862  point[1]*=0.5;
2863  point[2]*=0.5;
2864  break;
2865  case RefinementCase<3>::cut_xyz:
2866  point+=unit_cell_vertex(child_index);
2867  point*=0.5;
2868  break;
2869  default:
2870  Assert(false, ExcInternalError());
2871  }
2872 
2873  return point;
2874 }
2875 
2876 
2877 
2878 template <>
2879 inline
2880 Point<2>
2881 GeometryInfo<2>::child_to_cell_coordinates (const Point<2> &p,
2882  const unsigned int child_index,
2883  const RefinementCase<2> refine_case)
2884 {
2885  Assert (child_index < GeometryInfo<2>::n_children(refine_case),
2886  ExcIndexRange (child_index, 0, GeometryInfo<2>::n_children(refine_case)));
2887 
2888  Point<2> point=p;
2889  switch (refine_case)
2890  {
2891  case RefinementCase<2>::cut_x:
2892  if (child_index==1)
2893  point[0]+=1.0;
2894  point[0]*=0.5;
2895  break;
2896  case RefinementCase<2>::cut_y:
2897  if (child_index==1)
2898  point[1]+=1.0;
2899  point[1]*=0.5;
2900  break;
2901  case RefinementCase<2>::cut_xy:
2902  point+=unit_cell_vertex(child_index);
2903  point*=0.5;
2904  break;
2905  default:
2906  Assert(false, ExcInternalError());
2907  }
2908 
2909  return point;
2910 }
2911 
2912 
2913 
2914 template <int dim>
2915 inline
2916 Point<dim>
2917 GeometryInfo<dim>::child_to_cell_coordinates (const Point<dim> &/*p*/,
2918  const unsigned int /*child_index*/,
2919  const RefinementCase<dim> /*refine_case*/)
2920 {
2921  Assert (false, ExcNotImplemented());
2922  return Point<dim>();
2923 }
2924 
2925 
2926 
2927 template <int dim>
2928 inline
2929 bool
2930 GeometryInfo<dim>::is_inside_unit_cell (const Point<dim> &)
2931 {
2932  Assert(false, ExcNotImplemented());
2933  return false;
2934 }
2935 
2936 template <>
2937 inline
2938 bool
2939 GeometryInfo<1>::is_inside_unit_cell (const Point<1> &p)
2940 {
2941  return (p[0] >= 0.) && (p[0] <= 1.);
2942 }
2943 
2944 
2945 
2946 template <>
2947 inline
2948 bool
2949 GeometryInfo<2>::is_inside_unit_cell (const Point<2> &p)
2950 {
2951  return (p[0] >= 0.) && (p[0] <= 1.) &&
2952  (p[1] >= 0.) && (p[1] <= 1.);
2953 }
2954 
2955 
2956 
2957 template <>
2958 inline
2959 bool
2960 GeometryInfo<3>::is_inside_unit_cell (const Point<3> &p)
2961 {
2962  return (p[0] >= 0.) && (p[0] <= 1.) &&
2963  (p[1] >= 0.) && (p[1] <= 1.) &&
2964  (p[2] >= 0.) && (p[2] <= 1.);
2965 }
2966 
2967 
2968 
2969 template <int dim>
2970 inline
2971 bool
2972 GeometryInfo<dim>::is_inside_unit_cell (const Point<dim> &,
2973  const double)
2974 {
2975  Assert(false, ExcNotImplemented());
2976  return false;
2977 }
2978 
2979 template <>
2980 inline
2981 bool
2982 GeometryInfo<1>::is_inside_unit_cell (const Point<1> &p,
2983  const double eps)
2984 {
2985  return (p[0] >= -eps) && (p[0] <= 1.+eps);
2986 }
2987 
2988 
2989 
2990 template <>
2991 inline
2992 bool
2993 GeometryInfo<2>::is_inside_unit_cell (const Point<2> &p,
2994  const double eps)
2995 {
2996  const double l = -eps, u = 1+eps;
2997  return (p[0] >= l) && (p[0] <= u) &&
2998  (p[1] >= l) && (p[1] <= u);
2999 }
3000 
3001 
3002 
3003 template <>
3004 inline
3005 bool
3006 GeometryInfo<3>::is_inside_unit_cell (const Point<3> &p,
3007  const double eps)
3008 {
3009  const double l = -eps, u = 1.0+eps;
3010  return (p[0] >= l) && (p[0] <= u) &&
3011  (p[1] >= l) && (p[1] <= u) &&
3012  (p[2] >= l) && (p[2] <= u);
3013 }
3014 
3015 
3016 
3017 template <>
3018 inline
3019 unsigned int
3020 GeometryInfo<1>::line_to_cell_vertices (const unsigned int line,
3021  const unsigned int vertex)
3022 {
3023  (void)line;
3024  Assert (line<lines_per_cell, ExcIndexRange(line, 0, lines_per_cell));
3025  Assert (vertex<2, ExcIndexRange(vertex, 0, 2));
3026 
3027  return vertex;
3028 }
3029 
3030 
3031 template <>
3032 inline
3033 unsigned int
3034 GeometryInfo<2>::line_to_cell_vertices (const unsigned int line,
3035  const unsigned int vertex)
3036 {
3037  constexpr unsigned int cell_vertices[4][2] = {{0,2},{1,3},{0,1},{2,3}};
3038  return cell_vertices[line][vertex];
3039 }
3040 
3041 
3042 
3043 template <>
3044 inline
3045 unsigned int
3046 GeometryInfo<3>::line_to_cell_vertices (const unsigned int line,
3047  const unsigned int vertex)
3048 {
3049  Assert (line<lines_per_cell, ExcIndexRange(line, 0, lines_per_cell));
3050  Assert (vertex<2, ExcIndexRange(vertex, 0, 2));
3051 
3052  constexpr unsigned
3053  vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3054  {1, 3},
3055  {0, 1},
3056  {2, 3},
3057  {4, 6}, // top face
3058  {5, 7},
3059  {4, 5},
3060  {6, 7},
3061  {0, 4}, // connects of bottom
3062  {1, 5}, // top face
3063  {2, 6},
3064  {3, 7}
3065  };
3066  return vertices[line][vertex];
3067 }
3068 
3069 
3070 template <>
3071 inline
3072 unsigned int
3073 GeometryInfo<4>::line_to_cell_vertices (const unsigned int,
3074  const unsigned int)
3075 {
3076  Assert(false, ExcNotImplemented());
3077  return numbers::invalid_unsigned_int;
3078 }
3079 
3080 template <>
3081 inline
3082 unsigned int
3083 GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3084  const bool face_orientation,
3085  const bool face_flip,
3086  const bool face_rotation)
3087 {
3088  Assert(vertex<GeometryInfo<3>::vertices_per_face,
3089  ExcIndexRange(vertex,0,GeometryInfo<3>::vertices_per_face));
3090 
3091  // set up a table to make sure that
3092  // we handle non-standard faces correctly
3093  //
3094  // so set up a table that for each vertex (of
3095  // a quad in standard position) describes
3096  // which vertex to take
3097  //
3098  // first index: four vertices 0...3
3099  //
3100  // second index: face_orientation; 0:
3101  // opposite normal, 1: standard
3102  //
3103  // third index: face_flip; 0: standard, 1:
3104  // face rotated by 180 degrees
3105  //
3106  // forth index: face_rotation: 0: standard,
3107  // 1: face rotated by 90 degrees
3108 
3109  constexpr unsigned int vertex_translation[4][2][2][2] =
3110  {
3111  { { { 0, 2 }, // vertex 0, face_orientation=false, face_flip=false, face_rotation=false and true
3112  { 3, 1 }
3113  }, // vertex 0, face_orientation=false, face_flip=true, face_rotation=false and true
3114  { { 0, 2 }, // vertex 0, face_orientation=true, face_flip=false, face_rotation=false and true
3115  { 3, 1 }
3116  }
3117  },// vertex 0, face_orientation=true, face_flip=true, face_rotation=false and true
3118 
3119  { { { 2, 3 }, // vertex 1 ...
3120  { 1, 0 }
3121  },
3122  { { 1, 0 },
3123  { 2, 3 }
3124  }
3125  },
3126 
3127  { { { 1, 0 }, // vertex 2 ...
3128  { 2, 3 }
3129  },
3130  { { 2, 3 },
3131  { 1, 0 }
3132  }
3133  },
3134 
3135  { { { 3, 1 }, // vertex 3 ...
3136  { 0, 2 }
3137  },
3138  { { 3, 1 },
3139  { 0, 2 }
3140  }
3141  }
3142  };
3143 
3144  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3145 }
3146 
3147 
3148 
3149 template <int dim>
3150 inline
3151 unsigned int
3152 GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3153  const bool,
3154  const bool,
3155  const bool)
3156 {
3157  Assert(dim>1, ExcImpossibleInDim(dim));
3158  Assert(vertex<GeometryInfo<dim>::vertices_per_face,
3159  ExcIndexRange(vertex,0,GeometryInfo<dim>::vertices_per_face));
3160  return vertex;
3161 }
3162 
3163 template <int dim>
3164 inline
3165 unsigned int
3166 GeometryInfo<dim>::n_children(const RefinementCase<dim> &ref_case)
3167 {
3168  constexpr unsigned int n_children[RefinementCase<3>::cut_xyz+1]=
3169  {0, 2, 2, 4, 2, 4, 4, 8};
3170 
3171  return n_children[ref_case];
3172 }
3173 
3174 
3175 
3176 template <int dim>
3177 inline
3178 unsigned int
3179 GeometryInfo<dim>::n_subfaces(const internal::SubfaceCase<dim> &)
3180 {
3181  Assert(false, ExcNotImplemented());
3182  return 0;
3183 }
3184 
3185 template <>
3186 inline
3187 unsigned int
3188 GeometryInfo<1>::n_subfaces(const internal::SubfaceCase<1> &)
3189 {
3190  Assert(false, ExcImpossibleInDim(1));
3191  return 0;
3192 }
3193 
3194 template <>
3195 inline
3196 unsigned int
3197 GeometryInfo<2>::n_subfaces(const internal::SubfaceCase<2> &subface_case)
3198 {
3199  return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3200 }
3201 
3202 
3203 
3204 template <>
3205 inline
3206 unsigned int
3207 GeometryInfo<3>::n_subfaces(const internal::SubfaceCase<3> &subface_case)
3208 {
3209  const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic+1]=
3210  {0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3211  return nsubs[subface_case];
3212 }
3213 
3214 
3215 
3216 template <int dim>
3217 inline
3218 double
3219 GeometryInfo<dim>::subface_ratio(const internal::SubfaceCase<dim> &,
3220  const unsigned int)
3221 {
3222  Assert(false, ExcNotImplemented());
3223  return 0.;
3224 }
3225 
3226 template <>
3227 inline
3228 double
3229 GeometryInfo<1>::subface_ratio(const internal::SubfaceCase<1> &,
3230  const unsigned int)
3231 {
3232  return 1;
3233 }
3234 
3235 
3236 template <>
3237 inline
3238 double
3239 GeometryInfo<2>::subface_ratio(const internal::SubfaceCase<2> &subface_case,
3240  const unsigned int)
3241 {
3242  double ratio=1;
3243  switch (subface_case)
3244  {
3245  case internal::SubfaceCase<2>::case_none:
3246  // Here, an
3247  // Assert(false,ExcInternalError())
3248  // would be the right
3249  // choice, but
3250  // unfortunately the
3251  // current function is
3252  // also called for faces
3253  // without children (see
3254  // tests/fe/mapping.cc).
3255 // Assert(false, ExcMessage("Face has no subfaces."));
3256  // Furthermore, assign
3257  // following value as
3258  // otherwise the
3259  // bits/volume_x tests
3260  // break
3261  ratio=1./GeometryInfo<2>::max_children_per_face;
3262  break;
3263  case internal::SubfaceCase<2>::case_x:
3264  ratio=0.5;
3265  break;
3266  default:
3267  // there should be no
3268  // cases left
3269  Assert(false, ExcInternalError());
3270  break;
3271  }
3272 
3273  return ratio;
3274 }
3275 
3276 
3277 template <>
3278 inline
3279 double
3280 GeometryInfo<3>::subface_ratio(const internal::SubfaceCase<3> &subface_case,
3281  const unsigned int subface_no)
3282 {
3283  double ratio=1;
3284  switch (subface_case)
3285  {
3286  case internal::SubfaceCase<3>::case_none:
3287  // Here, an
3288  // Assert(false,ExcInternalError())
3289  // would be the right
3290  // choice, but
3291  // unfortunately the
3292  // current function is
3293  // also called for faces
3294  // without children (see
3295  // tests/bits/mesh_3d_16.cc). Add
3296  // following switch to
3297  // avoid diffs in
3298  // tests/bits/mesh_3d_16
3299  ratio=1./GeometryInfo<3>::max_children_per_face;
3300  break;
3301  case internal::SubfaceCase<3>::case_x:
3302  case internal::SubfaceCase<3>::case_y:
3303  ratio=0.5;
3304  break;
3305  case internal::SubfaceCase<3>::case_xy:
3306  case internal::SubfaceCase<3>::case_x1y2y:
3307  case internal::SubfaceCase<3>::case_y1x2x:
3308  ratio=0.25;
3309  break;
3310  case internal::SubfaceCase<3>::case_x1y:
3311  case internal::SubfaceCase<3>::case_y1x:
3312  if (subface_no<2)
3313  ratio=0.25;
3314  else
3315  ratio=0.5;
3316  break;
3317  case internal::SubfaceCase<3>::case_x2y:
3318  case internal::SubfaceCase<3>::case_y2x:
3319  if (subface_no==0)
3320  ratio=0.5;
3321  else
3322  ratio=0.25;
3323  break;
3324  default:
3325  // there should be no
3326  // cases left
3327  Assert(false, ExcInternalError());
3328  break;
3329  }
3330 
3331  return ratio;
3332 }
3333 
3334 
3335 
3336 template <int dim>
3337 RefinementCase<dim-1>
3338 inline
3339 GeometryInfo<dim>::face_refinement_case(const RefinementCase<dim> &,
3340  const unsigned int,
3341  const bool,
3342  const bool,
3343  const bool)
3344 {
3345  Assert (false, ExcNotImplemented());
3346  return RefinementCase<dim-1>::no_refinement;
3347 }
3348 
3349 template <>
3350 RefinementCase<0>
3351 inline
3352 GeometryInfo<1>::face_refinement_case(const RefinementCase<1> &,
3353  const unsigned int,
3354  const bool,
3355  const bool,
3356  const bool)
3357 {
3358  Assert(false, ExcImpossibleInDim(1));
3359 
3360  return RefinementCase<0>::no_refinement;
3361 }
3362 
3363 
3364 template <>
3365 inline
3366 RefinementCase<1>
3367 GeometryInfo<2>::face_refinement_case(const RefinementCase<2> &cell_refinement_case,
3368  const unsigned int face_no,
3369  const bool,
3370  const bool,
3371  const bool)
3372 {
3373  const unsigned int dim=2;
3374  Assert(cell_refinement_case<RefinementCase<dim>::isotropic_refinement+1,
3375  ExcIndexRange(cell_refinement_case, 0, RefinementCase<dim>::isotropic_refinement+1));
3376  Assert(face_no<GeometryInfo<dim>::faces_per_cell,
3377  ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
3378 
3379  const RefinementCase<dim-1>
3380  ref_cases[RefinementCase<dim>::isotropic_refinement+1][GeometryInfo<dim>::faces_per_cell/2]=
3381  {
3382  {
3383  RefinementCase<dim-1>::no_refinement, // no_refinement
3384  RefinementCase<dim-1>::no_refinement
3385  },
3386 
3387  {
3388  RefinementCase<dim-1>::no_refinement,
3389  RefinementCase<dim-1>::cut_x
3390  },
3391 
3392  {
3393  RefinementCase<dim-1>::cut_x,
3394  RefinementCase<dim-1>::no_refinement
3395  },
3396 
3397  {
3398  RefinementCase<dim-1>::cut_x, // cut_xy
3399  RefinementCase<dim-1>::cut_x
3400  }
3401  };
3402 
3403  return ref_cases[cell_refinement_case][face_no/2];
3404 }
3405 
3406 
3407 template <>
3408 inline
3409 RefinementCase<2>
3410 GeometryInfo<3>::face_refinement_case(const RefinementCase<3> &cell_refinement_case,
3411  const unsigned int face_no,
3412  const bool face_orientation,
3413  const bool /*face_flip*/,
3414  const bool face_rotation)
3415 {
3416  const unsigned int dim=3;
3417  Assert(cell_refinement_case<RefinementCase<dim>::isotropic_refinement+1,
3418  ExcIndexRange(cell_refinement_case, 0, RefinementCase<dim>::isotropic_refinement+1));
3419  Assert(face_no<GeometryInfo<dim>::faces_per_cell,
3420  ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
3421 
3422  const RefinementCase<dim-1>
3423  ref_cases[RefinementCase<dim>::isotropic_refinement+1][GeometryInfo<dim>::faces_per_cell/2]=
3424  {
3425  {
3426  RefinementCase<dim-1>::no_refinement, // no_refinement
3427  RefinementCase<dim-1>::no_refinement,
3428  RefinementCase<dim-1>::no_refinement
3429  },
3430 
3431  {
3432  RefinementCase<dim-1>::no_refinement, // cut_x
3433  RefinementCase<dim-1>::cut_y,
3434  RefinementCase<dim-1>::cut_x
3435  },
3436 
3437  {
3438  RefinementCase<dim-1>::cut_x, // cut_y
3439  RefinementCase<dim-1>::no_refinement,
3440  RefinementCase<dim-1>::cut_y
3441  },
3442 
3443  {
3444  RefinementCase<dim-1>::cut_x, // cut_xy
3445  RefinementCase<dim-1>::cut_y,
3446  RefinementCase<dim-1>::cut_xy
3447  },
3448 
3449  {
3450  RefinementCase<dim-1>::cut_y, // cut_z
3451  RefinementCase<dim-1>::cut_x,
3452  RefinementCase<dim-1>::no_refinement
3453  },
3454 
3455  {
3456  RefinementCase<dim-1>::cut_y, // cut_xz
3457  RefinementCase<dim-1>::cut_xy,
3458  RefinementCase<dim-1>::cut_x
3459  },
3460 
3461  {
3462  RefinementCase<dim-1>::cut_xy, // cut_yz
3463  RefinementCase<dim-1>::cut_x,
3464  RefinementCase<dim-1>::cut_y
3465  },
3466 
3467  {
3468  RefinementCase<dim-1>::cut_xy, // cut_xyz
3469  RefinementCase<dim-1>::cut_xy,
3470  RefinementCase<dim-1>::cut_xy
3471  },
3472  };
3473 
3474  const RefinementCase<dim-1> ref_case=ref_cases[cell_refinement_case][face_no/2];
3475 
3476  const RefinementCase<dim-1> flip[4]=
3477  {
3478  RefinementCase<dim-1>::no_refinement,
3479  RefinementCase<dim-1>::cut_y,
3480  RefinementCase<dim-1>::cut_x,
3481  RefinementCase<dim-1>::cut_xy
3482  };
3483 
3484  // correct the ref_case for face_orientation
3485  // and face_rotation. for face_orientation,
3486  // 'true' is the default value whereas for
3487  // face_rotation, 'false' is standard. If
3488  // <tt>face_rotation==face_orientation</tt>,
3489  // then one of them is non-standard and we
3490  // have to swap cut_x and cut_y, otherwise no
3491  // change is necessary. face_flip has no
3492  // influence. however, in order to keep the
3493  // interface consistent with other functions,
3494  // we still include it as an argument to this
3495  // function
3496  return (face_orientation==face_rotation) ? flip[ref_case] : ref_case;
3497 }
3498 
3499 
3500 
3501 template <int dim>
3502 inline
3503 RefinementCase<1>
3504 GeometryInfo<dim>::line_refinement_case(const RefinementCase<dim> &,
3505  const unsigned int)
3506 {
3507  Assert(false, ExcNotImplemented());
3508  return RefinementCase<1>::no_refinement;
3509 }
3510 
3511 template <>
3512 inline
3513 RefinementCase<1>
3514 GeometryInfo<1>::line_refinement_case(const RefinementCase<1> &cell_refinement_case,
3515  const unsigned int line_no)
3516 {
3517  (void)line_no;
3518  const unsigned int dim = 1;
3519  (void)dim;
3520  Assert(cell_refinement_case<RefinementCase<dim>::isotropic_refinement+1,
3521  ExcIndexRange(cell_refinement_case, 0, RefinementCase<dim>::isotropic_refinement+1));
3522  Assert(line_no<GeometryInfo<dim>::lines_per_cell,
3523  ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
3524 
3525  return cell_refinement_case;
3526 }
3527 
3528 
3529 template <>
3530 inline
3531 RefinementCase<1>
3532 GeometryInfo<2>::line_refinement_case(const RefinementCase<2> &cell_refinement_case,
3533  const unsigned int line_no)
3534 {
3535  // Assertions are in face_refinement_case()
3536  return face_refinement_case(cell_refinement_case, line_no);
3537 }
3538 
3539 
3540 template <>
3541 inline
3542 RefinementCase<1>
3543 GeometryInfo<3>::line_refinement_case(const RefinementCase<3> &cell_refinement_case,
3544  const unsigned int line_no)
3545 {
3546  const unsigned int dim=3;
3547  Assert(cell_refinement_case<RefinementCase<dim>::isotropic_refinement+1,
3548  ExcIndexRange(cell_refinement_case, 0, RefinementCase<dim>::isotropic_refinement+1));
3549  Assert(line_no<GeometryInfo<dim>::lines_per_cell,
3550  ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
3551 
3552  // array indicating, which simple refine
3553  // case cuts a line in direction x, y or
3554  // z. For example, cut_y and everything
3555  // containing cut_y (cut_xy, cut_yz,
3556  // cut_xyz) cuts lines, which are in y
3557  // direction.
3558  const RefinementCase<dim>
3559  cut_one[dim] =
3560  {
3561  RefinementCase<dim>::cut_x,
3562  RefinementCase<dim>::cut_y,
3563  RefinementCase<dim>::cut_z
3564  };
3565 
3566  // order the direction of lines
3567  // 0->x, 1->y, 2->z
3568  const unsigned int direction[lines_per_cell]=
3569  {1,1,0,0,1,1,0,0,2,2,2,2};
3570 
3571  return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3572  RefinementCase<1>::cut_x : RefinementCase<1>::no_refinement);
3573 }
3574 
3575 
3576 
3577 template <int dim>
3578 inline
3579 RefinementCase<dim>
3580 GeometryInfo<dim>::min_cell_refinement_case_for_face_refinement(const RefinementCase<dim-1> &,
3581  const unsigned int,
3582  const bool,
3583  const bool,
3584  const bool)
3585 {
3586  Assert(false, ExcNotImplemented());
3587 
3588  return RefinementCase<dim>::no_refinement;
3589 }
3590 
3591 template <>
3592 inline
3593 RefinementCase<1>
3594 GeometryInfo<1>::min_cell_refinement_case_for_face_refinement(const RefinementCase<0> &,
3595  const unsigned int,
3596  const bool,
3597  const bool,
3598  const bool)
3599 {
3600  const unsigned int dim = 1;
3601  Assert(false, ExcImpossibleInDim(dim));
3602 
3603  return RefinementCase<dim>::no_refinement;
3604 }
3605 
3606 
3607 template <>
3608 inline
3609 RefinementCase<2>
3610 GeometryInfo<2>::min_cell_refinement_case_for_face_refinement(const RefinementCase<1> &face_refinement_case,
3611  const unsigned int face_no,
3612  const bool,
3613  const bool,
3614  const bool)
3615 {
3616  const unsigned int dim = 2;
3617  Assert(face_refinement_case<RefinementCase<dim-1>::isotropic_refinement+1,
3618  ExcIndexRange(face_refinement_case, 0, RefinementCase<dim-1>::isotropic_refinement+1));
3619  Assert(face_no<GeometryInfo<dim>::faces_per_cell,
3620  ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
3621 
3622  if (face_refinement_case==RefinementCase<dim>::cut_x)
3623  return (face_no/2) ? RefinementCase<dim>::cut_x : RefinementCase<dim>::cut_y;
3624  else
3625  return RefinementCase<dim>::no_refinement;
3626 }
3627 
3628 
3629 template <>
3630 inline
3631 RefinementCase<3>
3632 GeometryInfo<3>::min_cell_refinement_case_for_face_refinement(const RefinementCase<2> &face_refinement_case,
3633  const unsigned int face_no,
3634  const bool face_orientation,
3635  const bool /*face_flip*/,
3636  const bool face_rotation)
3637 {
3638  const unsigned int dim=3;
3639  Assert(face_refinement_case<RefinementCase<dim-1>::isotropic_refinement+1,
3640  ExcIndexRange(face_refinement_case, 0, RefinementCase<dim-1>::isotropic_refinement+1));
3641  Assert(face_no<GeometryInfo<dim>::faces_per_cell,
3642  ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
3643 
3644  const RefinementCase<2> flip[4]=
3645  {
3646  RefinementCase<2>::no_refinement,
3647  RefinementCase<2>::cut_y,
3648  RefinementCase<2>::cut_x,
3649  RefinementCase<2>::cut_xy
3650  };
3651 
3652  // correct the face_refinement_case for
3653  // face_orientation and face_rotation. for
3654  // face_orientation, 'true' is the default
3655  // value whereas for face_rotation, 'false'
3656  // is standard. If
3657  // <tt>face_rotation==face_orientation</tt>,
3658  // then one of them is non-standard and we
3659  // have to swap cut_x and cut_y, otherwise no
3660  // change is necessary. face_flip has no
3661  // influence. however, in order to keep the
3662  // interface consistent with other functions,
3663  // we still include it as an argument to this
3664  // function
3665  const RefinementCase<dim-1> std_face_ref = (face_orientation==face_rotation) ? flip[face_refinement_case] : face_refinement_case;
3666 
3667  const RefinementCase<dim> face_to_cell[3][4]=
3668  {
3669  {
3670  RefinementCase<dim>::no_refinement, // faces 0 and 1
3671  RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3672  RefinementCase<dim>::cut_z,
3673  RefinementCase<dim>::cut_yz
3674  },
3675 
3676  {
3677  RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are "exchanged on faces 2 and 3")
3678  RefinementCase<dim>::cut_z,
3679  RefinementCase<dim>::cut_x,
3680  RefinementCase<dim>::cut_xz
3681  },
3682 
3683  {
3684  RefinementCase<dim>::no_refinement, // faces 4 and 5
3685  RefinementCase<dim>::cut_x,
3686  RefinementCase<dim>::cut_y,
3687  RefinementCase<dim>::cut_xy
3688  }
3689  };
3690 
3691  return face_to_cell[face_no/2][std_face_ref];
3692 }
3693 
3694 
3695 
3696 template <int dim>
3697 inline
3698 RefinementCase<dim>
3699 GeometryInfo<dim>::min_cell_refinement_case_for_line_refinement(const unsigned int)
3700 {
3701  Assert(false, ExcNotImplemented());
3702 
3703  return RefinementCase<dim>::no_refinement;
3704 }
3705 
3706 template <>
3707 inline
3708 RefinementCase<1>
3709 GeometryInfo<1>::min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
3710 {
3711  (void)line_no;
3712  Assert(line_no==0, ExcIndexRange(line_no,0,1));
3713 
3714  return RefinementCase<1>::cut_x;
3715 }
3716 
3717 
3718 template <>
3719 inline
3720 RefinementCase<2>
3721 GeometryInfo<2>::min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
3722 {
3723  const unsigned int dim = 2;
3724  (void)dim;
3725  Assert(line_no<GeometryInfo<dim>::lines_per_cell,
3726  ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
3727 
3728  return (line_no/2) ? RefinementCase<2>::cut_x : RefinementCase<2>::cut_y;
3729 }
3730 
3731 
3732 template <>
3733 inline
3734 RefinementCase<3>
3735 GeometryInfo<3>::min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
3736 {
3737  const unsigned int dim=3;
3738  Assert(line_no<GeometryInfo<dim>::lines_per_cell,
3739  ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
3740 
3741  const RefinementCase<dim> ref_cases[6]=
3742  {
3743  RefinementCase<dim>::cut_y, // lines 0 and 1
3744  RefinementCase<dim>::cut_x, // lines 2 and 3
3745  RefinementCase<dim>::cut_y, // lines 4 and 5
3746  RefinementCase<dim>::cut_x, // lines 6 and 7
3747  RefinementCase<dim>::cut_z, // lines 8 and 9
3748  RefinementCase<dim>::cut_z
3749  }; // lines 10 and 11
3750 
3751  return ref_cases[line_no/2];
3752 }
3753 
3754 
3755 
3756 template <>
3757 inline
3758 unsigned int
3759 GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
3760  const bool face_orientation,
3761  const bool face_flip,
3762  const bool face_rotation)
3763 {
3764  Assert(vertex<GeometryInfo<3>::vertices_per_face,
3765  ExcIndexRange(vertex,0,GeometryInfo<3>::vertices_per_face));
3766 
3767  // set up a table to make sure that
3768  // we handle non-standard faces correctly
3769  //
3770  // so set up a table that for each vertex (of
3771  // a quad in standard position) describes
3772  // which vertex to take
3773  //
3774  // first index: four vertices 0...3
3775  //
3776  // second index: face_orientation; 0:
3777  // opposite normal, 1: standard
3778  //
3779  // third index: face_flip; 0: standard, 1:
3780  // face rotated by 180 degrees
3781  //
3782  // forth index: face_rotation: 0: standard,
3783  // 1: face rotated by 90 degrees
3784 
3785  const unsigned int vertex_translation[4][2][2][2] =
3786  {
3787  { { { 0, 2 }, // vertex 0, face_orientation=false, face_flip=false, face_rotation=false and true
3788  { 3, 1 }
3789  }, // vertex 0, face_orientation=false, face_flip=true, face_rotation=false and true
3790  { { 0, 1 }, // vertex 0, face_orientation=true, face_flip=false, face_rotation=false and true
3791  { 3, 2 }
3792  }
3793  },// vertex 0, face_orientation=true, face_flip=true, face_rotation=false and true
3794 
3795  { { { 2, 3 }, // vertex 1 ...
3796  { 1, 0 }
3797  },
3798  { { 1, 3 },
3799  { 2, 0 }
3800  }
3801  },
3802 
3803  { { { 1, 0 }, // vertex 2 ...
3804  { 2, 3 }
3805  },
3806  { { 2, 0 },
3807  { 1, 3 }
3808  }
3809  },
3810 
3811  { { { 3, 1 }, // vertex 3 ...
3812  { 0, 2 }
3813  },
3814  { { 3, 2 },
3815  { 0, 1 }
3816  }
3817  }
3818  };
3819 
3820  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3821 }
3822 
3823 
3824 
3825 template <int dim>
3826 inline
3827 unsigned int
3828 GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
3829  const bool,
3830  const bool,
3831  const bool)
3832 {
3833  Assert(dim>1, ExcImpossibleInDim(dim));
3834  Assert(vertex<GeometryInfo<dim>::vertices_per_face,
3835  ExcIndexRange(vertex,0,GeometryInfo<dim>::vertices_per_face));
3836  return vertex;
3837 }
3838 
3839 
3840 
3841 template <>
3842 inline
3843 unsigned int
3844 GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
3845  const bool face_orientation,
3846  const bool face_flip,
3847  const bool face_rotation)
3848 {
3849  Assert(line<GeometryInfo<3>::lines_per_face,
3850  ExcIndexRange(line,0,GeometryInfo<3>::lines_per_face));
3851 
3852 
3853  // make sure we handle
3854  // non-standard faces correctly
3855  //
3856  // so set up a table that for each line (of a
3857  // quad) describes which line to take
3858  //
3859  // first index: four lines 0...3
3860  //
3861  // second index: face_orientation; 0:
3862  // opposite normal, 1: standard
3863  //
3864  // third index: face_flip; 0: standard, 1:
3865  // face rotated by 180 degrees
3866  //
3867  // forth index: face_rotation: 0: standard,
3868  // 1: face rotated by 90 degrees
3869 
3870  const unsigned int line_translation[4][2][2][2] =
3871  {
3872  { { { 2, 0 }, // line 0, face_orientation=false, face_flip=false, face_rotation=false and true
3873  { 3, 1 }
3874  }, // line 0, face_orientation=false, face_flip=true, face_rotation=false and true
3875  { { 0, 3 }, // line 0, face_orientation=true, face_flip=false, face_rotation=false and true
3876  { 1, 2 }
3877  }
3878  },// line 0, face_orientation=true, face_flip=true, face_rotation=false and true
3879 
3880  { { { 3, 1 }, // line 1 ...
3881  { 2, 0 }
3882  },
3883  { { 1, 2 },
3884  { 0, 3 }
3885  }
3886  },
3887 
3888  { { { 0, 3 }, // line 2 ...
3889  { 1, 2 }
3890  },
3891  { { 2, 0 },
3892  { 3, 1 }
3893  }
3894  },
3895 
3896  { { { 1, 2 }, // line 3 ...
3897  { 0, 3 }
3898  },
3899  { { 3, 1 },
3900  { 2, 0 }
3901  }
3902  }
3903  };
3904 
3905  return line_translation[line][face_orientation][face_flip][face_rotation];
3906 }
3907 
3908 
3909 
3910 template <int dim>
3911 inline
3912 unsigned int
3913 GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
3914  const bool,
3915  const bool,
3916  const bool)
3917 {
3918  Assert(false, ExcNotImplemented());
3919  return line;
3920 }
3921 
3922 
3923 
3924 template <>
3925 inline
3926 unsigned int
3927 GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
3928  const bool face_orientation,
3929  const bool face_flip,
3930  const bool face_rotation)
3931 {
3932  Assert(line<GeometryInfo<3>::lines_per_face,
3933  ExcIndexRange(line,0,GeometryInfo<3>::lines_per_face));
3934 
3935 
3936  // make sure we handle
3937  // non-standard faces correctly
3938  //
3939  // so set up a table that for each line (of a
3940  // quad) describes which line to take
3941  //
3942  // first index: four lines 0...3
3943  //
3944  // second index: face_orientation; 0:
3945  // opposite normal, 1: standard
3946  //
3947  // third index: face_flip; 0: standard, 1:
3948  // face rotated by 180 degrees
3949  //
3950  // forth index: face_rotation: 0: standard,
3951  // 1: face rotated by 90 degrees
3952 
3953  const unsigned int line_translation[4][2][2][2] =
3954  {
3955  { { { 2, 0 }, // line 0, face_orientation=false, face_flip=false, face_rotation=false and true
3956  { 3, 1 }
3957  }, // line 0, face_orientation=false, face_flip=true, face_rotation=false and true
3958  { { 0, 2 }, // line 0, face_orientation=true, face_flip=false, face_rotation=false and true
3959  { 1, 3 }
3960  }
3961  },// line 0, face_orientation=true, face_flip=true, face_rotation=false and true
3962 
3963  { { { 3, 1 }, // line 1 ...
3964  { 2, 0 }
3965  },
3966  { { 1, 3 },
3967  { 0, 2 }
3968  }
3969  },
3970 
3971  { { { 0, 3 }, // line 2 ...
3972  { 1, 2 }
3973  },
3974  { { 2, 1 },
3975  { 3, 0 }
3976  }
3977  },
3978 
3979  { { { 1, 2 }, // line 3 ...
3980  { 0, 3 }
3981  },
3982  { { 3, 0 },
3983  { 2, 1 }
3984  }
3985  }
3986  };
3987 
3988  return line_translation[line][face_orientation][face_flip][face_rotation];
3989 }
3990 
3991 
3992 
3993 template <int dim>
3994 inline
3995 unsigned int
3996 GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
3997  const bool,
3998  const bool,
3999  const bool)
4000 {
4001  Assert(false, ExcNotImplemented());
4002  return line;
4003 }
4004 
4005 
4006 
4007 template <>
4008 inline
4009 unsigned int
4010 GeometryInfo<1>::child_cell_on_face (const RefinementCase<1> &,
4011  const unsigned int face,
4012  const unsigned int subface,
4013  const bool, const bool, const bool,
4014  const RefinementCase<0> &)
4015 {
4016  (void)subface;
4017  Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
4018  Assert (subface<max_children_per_face,
4019  ExcIndexRange(subface, 0, max_children_per_face));
4020 
4021  return face;
4022 }
4023 
4024 
4025 
4026 template <>
4027 inline
4028 unsigned int
4029 GeometryInfo<2>::child_cell_on_face (const RefinementCase<2> &ref_case,
4030  const unsigned int face,
4031  const unsigned int subface,
4032  const bool /*face_orientation*/,
4033  const bool face_flip,
4034  const bool /*face_rotation*/,
4035  const RefinementCase<1> &)
4036 {
4037  Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
4038  Assert (subface<max_children_per_face,
4039  ExcIndexRange(subface, 0, max_children_per_face));
4040 
4041  // always return the child adjacent to the specified
4042  // subface. if the face of a cell is not refined, don't
4043  // throw an assertion but deliver the child adjacent to
4044  // the face nevertheless, i.e. deliver the child of
4045  // this cell adjacent to the subface of a possibly
4046  // refined neighbor. this simplifies setting neighbor
4047  // information in execute_refinement.
4048  const unsigned int
4049  subcells[2][RefinementCase<2>::isotropic_refinement][faces_per_cell][max_children_per_face] =
4050  {
4051  {
4052  // Normal orientation (face_flip = false)
4053  {{0,0},{1,1},{0,1},{0,1}}, // cut_x
4054  {{0,1},{0,1},{0,0},{1,1}}, // cut_y
4055  {{0,2},{1,3},{0,1},{2,3}} // cut_xy, i.e., isotropic
4056  },
4057  {
4058  // Flipped orientation (face_flip = true)
4059  {{0,0},{1,1},{1,0},{1,0}}, // cut_x
4060  {{1,0},{1,0},{0,0},{1,1}}, // cut_y
4061  {{2,0},{3,1},{1,0},{3,2}} // cut_xy, i.e., isotropic
4062  }
4063  };
4064 
4065  return subcells[face_flip][ref_case-1][face][subface];
4066 }
4067 
4068 
4069 
4070 template <>
4071 inline
4072 unsigned int
4073 GeometryInfo<3>::child_cell_on_face (const RefinementCase<3> &ref_case,
4074  const unsigned int face,
4075  const unsigned int subface,
4076  const bool face_orientation,
4077  const bool face_flip,
4078  const bool face_rotation,
4079  const RefinementCase<2> &face_ref_case)
4080 {
4081  const unsigned int dim = 3;
4082 
4083  Assert (ref_case>RefinementCase<dim-1>::no_refinement, ExcMessage("Cell has no children."));
4084  Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
4085  Assert (subface<GeometryInfo<dim-1>::n_children(face_ref_case) ||
4086  (subface==0 && face_ref_case==RefinementCase<dim-1>::no_refinement),
4087  ExcIndexRange(subface, 0, GeometryInfo<2>::n_children(face_ref_case)));
4088 
4089  // invalid number used for invalid cases,
4090  // e.g. when the children are more refined at
4091  // a given face than the face itself
4092  const unsigned int e = numbers::invalid_unsigned_int;
4093 
4094  // the whole process of finding a child cell
4095  // at a given subface considering the
4096  // possibly anisotropic refinement cases of
4097  // the cell and the face as well as
4098  // orientation, flip and rotation of the face
4099  // is quite complicated. thus, we break it
4100  // down into several steps.
4101 
4102  // first step: convert the given face refine
4103  // case to a face refine case concerning the
4104  // face in standard orientation (, flip and
4105  // rotation). This only affects cut_x and
4106  // cut_y
4107  const RefinementCase<dim-1> flip[4]=
4108  {
4109  RefinementCase<dim-1>::no_refinement,
4110  RefinementCase<dim-1>::cut_y,
4111  RefinementCase<dim-1>::cut_x,
4112  RefinementCase<dim-1>::cut_xy
4113  };
4114  // for face_orientation, 'true' is the
4115  // default value whereas for face_rotation,
4116  // 'false' is standard. If
4117  // <tt>face_rotation==face_orientation</tt>,
4118  // then one of them is non-standard and we
4119  // have to swap cut_x and cut_y, otherwise no
4120  // change is necessary.
4121  const RefinementCase<dim-1> std_face_ref = (face_orientation==face_rotation) ? flip[face_ref_case] : face_ref_case;
4122 
4123  // second step: convert the given subface
4124  // index to the one for a standard face
4125  // respecting face_orientation, face_flip and
4126  // face_rotation
4127 
4128  // first index: face_ref_case
4129  // second index: face_orientation
4130  // third index: face_flip
4131  // forth index: face_rotation
4132  // fifth index: subface index
4133  const unsigned int subface_exchange[4][2][2][2][4]=
4134  {
4135  // no_refinement (subface 0 stays 0,
4136  // all others are invalid)
4137  { { { {0,e,e,e},
4138  {0,e,e,e}
4139  },
4140  { {0,e,e,e},
4141  {0,e,e,e}
4142  }
4143  },
4144  { { {0,e,e,e},
4145  {0,e,e,e}
4146  },
4147  { {0,e,e,e},
4148  {0,e,e,e}
4149  }
4150  }
4151  },
4152  // cut_x (here, if the face is only
4153  // rotated OR only falsely oriented,
4154  // then subface 0 of the non-standard
4155  // face does NOT correspond to one of
4156  // the subfaces of a standard
4157  // face. Thus we indicate the subface
4158  // which is located at the lower left
4159  // corner (the origin of the face's
4160  // local coordinate system) with
4161  // '0'. The rest of this issue is
4162  // taken care of using the above
4163  // conversion to a 'standard face
4164  // refine case')
4165  { { { {0,1,e,e},
4166  {0,1,e,e}
4167  },
4168  { {1,0,e,e},
4169  {1,0,e,e}
4170  }
4171  },
4172  { { {0,1,e,e},
4173  {0,1,e,e}
4174  },
4175  { {1,0,e,e},
4176  {1,0,e,e}
4177  }
4178  }
4179  },
4180  // cut_y (the same applies as for
4181  // cut_x)
4182  { { { {0,1,e,e},
4183  {1,0,e,e}
4184  },
4185  { {1,0,e,e},
4186  {0,1,e,e}
4187  }
4188  },
4189  { { {0,1,e,e},
4190  {1,0,e,e}
4191  },
4192  { {1,0,e,e},
4193  {0,1,e,e}
4194  }
4195  }
4196  },
4197  // cut_xyz: this information is
4198  // identical to the information
4199  // returned by
4200  // GeometryInfo<3>::real_to_standard_face_vertex()
4201  { { { {0,2,1,3}, // face_orientation=false, face_flip=false, face_rotation=false, subfaces 0,1,2,3
4202  {2,3,0,1}
4203  }, // face_orientation=false, face_flip=false, face_rotation=true, subfaces 0,1,2,3
4204  { {3,1,2,0}, // face_orientation=false, face_flip=true, face_rotation=false, subfaces 0,1,2,3
4205  {1,0,3,2}
4206  }
4207  }, // face_orientation=false, face_flip=true, face_rotation=true, subfaces 0,1,2,3
4208  { { {0,1,2,3}, // face_orientation=true, face_flip=false, face_rotation=false, subfaces 0,1,2,3
4209  {1,3,0,2}
4210  }, // face_orientation=true, face_flip=false, face_rotation=true, subfaces 0,1,2,3
4211  { {3,2,1,0}, // face_orientation=true, face_flip=true, face_rotation=false, subfaces 0,1,2,3
4212  {2,0,3,1}
4213  }
4214  }
4215  }
4216  };// face_orientation=true, face_flip=true, face_rotation=true, subfaces 0,1,2,3
4217 
4218  const unsigned int std_subface=subface_exchange
4219  [face_ref_case]
4220  [face_orientation]
4221  [face_flip]
4222  [face_rotation]
4223  [subface];
4224  Assert (std_subface!=e, ExcInternalError());
4225 
4226  // third step: these are the children, which
4227  // can be found at the given subfaces of an
4228  // isotropically refined (standard) face
4229  //
4230  // first index: (refinement_case-1)
4231  // second index: face_index
4232  // third index: subface_index (isotropic refinement)
4233  const unsigned int
4234  iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell][max_children_per_face] =
4235  {
4236  // cut_x
4237  { {0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4238  {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4239  {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4240  {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4241  {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4242  {0, 1, 0, 1}
4243  }, // face 5, subfaces 0,1,2,3
4244  // cut_y
4245  { {0, 1, 0, 1},
4246  {0, 1, 0, 1},
4247  {0, 0, 0, 0},
4248  {1, 1, 1, 1},
4249  {0, 0, 1, 1},
4250  {0, 0, 1, 1}
4251  },
4252  // cut_xy
4253  { {0, 2, 0, 2},
4254  {1, 3, 1, 3},
4255  {0, 0, 1, 1},
4256  {2, 2, 3, 3},
4257  {0, 1, 2, 3},
4258  {0, 1, 2, 3}
4259  },
4260  // cut_z
4261  { {0, 0, 1, 1},
4262  {0, 0, 1, 1},
4263  {0, 1, 0, 1},
4264  {0, 1, 0, 1},
4265  {0, 0, 0, 0},
4266  {1, 1, 1, 1}
4267  },
4268  // cut_xz
4269  { {0, 0, 1, 1},
4270  {2, 2, 3, 3},
4271  {0, 1, 2, 3},
4272  {0, 1, 2, 3},
4273  {0, 2, 0, 2},
4274  {1, 3, 1, 3}
4275  },
4276  // cut_yz
4277  { {0, 1, 2, 3},
4278  {0, 1, 2, 3},
4279  {0, 2, 0, 2},
4280  {1, 3, 1, 3},
4281  {0, 0, 1, 1},
4282  {2, 2, 3, 3}
4283  },
4284  // cut_xyz
4285  { {0, 2, 4, 6},
4286  {1, 3, 5, 7},
4287  {0, 4, 1, 5},
4288  {2, 6, 3, 7},
4289  {0, 1, 2, 3},
4290  {4, 5, 6, 7}
4291  }
4292  };
4293 
4294  // forth step: check, whether the given face
4295  // refine case is valid for the given cell
4296  // refine case. this is the case, if the
4297  // given face refine case is at least as
4298  // refined as the face is for the given cell
4299  // refine case
4300 
4301  // note, that we are considering standard
4302  // face refinement cases here and thus must
4303  // not pass the given orientation, flip and
4304  // rotation flags
4305  if ((std_face_ref & face_refinement_case(ref_case, face))
4306  == face_refinement_case(ref_case, face))
4307  {
4308  // all is fine. for anisotropic face
4309  // refine cases, select one of the
4310  // isotropic subfaces which neighbors the
4311  // same child
4312 
4313  // first index: (standard) face refine case
4314  // second index: subface index
4315  const unsigned int equivalent_iso_subface[4][4]=
4316  {
4317  {0,e,e,e}, // no_refinement
4318  {0,3,e,e}, // cut_x
4319  {0,3,e,e}, // cut_y
4320  {0,1,2,3}
4321  }; // cut_xy
4322 
4323  const unsigned int equ_std_subface
4324  =equivalent_iso_subface[std_face_ref][std_subface];
4325  Assert (equ_std_subface!=e, ExcInternalError());
4326 
4327  return iso_children[ref_case-1][face][equ_std_subface];
4328  }
4329  else
4330  {
4331  // the face_ref_case was too coarse,
4332  // throw an error
4333  Assert(false,
4334  ExcMessage("The face RefineCase is too coarse "
4335  "for the given cell RefineCase."));
4336  }
4337  // we only get here in case of an error
4338  return e;
4339 }
4340 
4341 
4342 
4343 template <>
4344 inline
4345 unsigned int
4346 GeometryInfo<4>::child_cell_on_face (const RefinementCase<4> &,
4347  const unsigned int,
4348  const unsigned int,
4349  const bool, const bool, const bool,
4350  const RefinementCase<3> &)
4351 {
4352  Assert(false, ExcNotImplemented());
4353  return numbers::invalid_unsigned_int;
4354 }
4355 
4356 
4357 
4358 template <>
4359 inline
4360 unsigned int
4361 GeometryInfo<1>::face_to_cell_lines (const unsigned int face,
4362  const unsigned int line,
4363  const bool, const bool, const bool)
4364 {
4365  (void)face;
4366  (void)line;
4367  Assert (face+1<faces_per_cell+1, ExcIndexRange(face, 0, faces_per_cell));
4368  Assert (line+1<lines_per_face+1, ExcIndexRange(line, 0, lines_per_face));
4369 
4370  // There is only a single line, so
4371  // it must be this.
4372  return 0;
4373 }
4374 
4375 
4376 
4377 template <>
4378 inline
4379 unsigned int
4380 GeometryInfo<2>::face_to_cell_lines (const unsigned int face,
4381  const unsigned int line,
4382  const bool, const bool, const bool)
4383 {
4384  (void)line;
4385  Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
4386  Assert (line<lines_per_face, ExcIndexRange(line, 0, lines_per_face));
4387 
4388  // The face is a line itself.
4389  return face;
4390 }
4391 
4392 
4393 
4394 template <>
4395 inline
4396 unsigned int
4397 GeometryInfo<3>::face_to_cell_lines (const unsigned int face,
4398  const unsigned int line,
4399  const bool face_orientation,
4400  const bool face_flip,
4401  const bool face_rotation)
4402 {
4403  Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
4404  Assert (line<lines_per_face, ExcIndexRange(line, 0, lines_per_face));
4405 
4406  const unsigned
4407  lines[faces_per_cell][lines_per_face] = {{8,10, 0, 4}, // left face
4408  {9,11, 1, 5}, // right face
4409  {2, 6, 8, 9}, // front face
4410  {3, 7,10,11}, // back face
4411  {0, 1, 2, 3}, // bottom face
4412  {4, 5, 6, 7}
4413  };// top face
4414  return lines[face][real_to_standard_face_line(line,
4415  face_orientation,
4416  face_flip,
4417  face_rotation)];
4418 }
4419 
4420 
4421 
4422 template <int dim>
4423 inline
4424 unsigned int
4425 GeometryInfo<dim>::face_to_cell_lines (const unsigned int,
4426  const unsigned int,
4427  const bool, const bool, const bool)
4428 {
4429  Assert(false, ExcNotImplemented());
4430  return numbers::invalid_unsigned_int;
4431 }
4432 
4433 
4434 
4435 template <int dim>
4436 inline
4437 unsigned int
4438 GeometryInfo<dim>::face_to_cell_vertices (const unsigned int face,
4439  const unsigned int vertex,
4440  const bool face_orientation,
4441  const bool face_flip,
4442  const bool face_rotation)
4443 {
4444  return child_cell_on_face(RefinementCase<dim>::isotropic_refinement, face, vertex,
4445  face_orientation, face_flip, face_rotation);
4446 }
4447 
4448 
4449 
4450 template <int dim>
4451 inline
4452 Point<dim>
4453 GeometryInfo<dim>::project_to_unit_cell (const Point<dim> &q)
4454 {
4455  Point<dim> p = q;
4456  for (unsigned int i=0; i<dim; i++)
4457  if (p[i] < 0.) p[i] = 0.;
4458  else if (p[i] > 1.) p[i] = 1.;
4459 
4460  return p;
4461 }
4462 
4463 
4464 
4465 template <int dim>
4466 inline
4467 double
4468 GeometryInfo<dim>::distance_to_unit_cell (const Point<dim> &p)
4469 {
4470  double result = 0.0;
4471 
4472  for (unsigned int i=0; i<dim; i++)
4473  if ((-p[i]) > result)
4474  result = -p[i];
4475  else if ((p[i]-1.) > result)
4476  result = (p[i] - 1.);
4477 
4478  return result;
4479 }
4480 
4481 
4482 
4483 template <int dim>
4484 inline
4485 double
4486 GeometryInfo<dim>::
4487 d_linear_shape_function (const Point<dim> &xi,
4488  const unsigned int i)
4489 {
4490  Assert (i < GeometryInfo<dim>::vertices_per_cell,
4491  ExcIndexRange (i, 0, GeometryInfo<dim>::vertices_per_cell));
4492 
4493  switch (dim)
4494  {
4495  case 1:
4496  {
4497  const double x = xi[0];
4498  switch (i)
4499  {
4500  case 0:
4501  return 1-x;
4502  case 1:
4503  return x;
4504  }
4505  break;
4506  }
4507 
4508  case 2:
4509  {
4510  const double x = xi[0];
4511  const double y = xi[1];
4512  switch (i)
4513  {
4514  case 0:
4515  return (1-x)*(1-y);
4516  case 1:
4517  return x*(1-y);
4518  case 2:
4519  return (1-x)*y;
4520  case 3:
4521  return x*y;
4522  }
4523  break;
4524  }
4525 
4526  case 3:
4527  {
4528  const double x = xi[0];
4529  const double y = xi[1];
4530  const double z = xi[2];
4531  switch (i)
4532  {
4533  case 0:
4534  return (1-x)*(1-y)*(1-z);
4535  case 1:
4536  return x*(1-y)*(1-z);
4537  case 2:
4538  return (1-x)*y*(1-z);
4539  case 3:
4540  return x*y*(1-z);
4541  case 4:
4542  return (1-x)*(1-y)*z;
4543  case 5:
4544  return x*(1-y)*z;
4545  case 6:
4546  return (1-x)*y*z;
4547  case 7:
4548  return x*y*z;
4549  }
4550  break;
4551  }
4552 
4553  default:
4554  Assert (false, ExcNotImplemented());
4555  }
4556  return -1e9;
4557 }
4558 
4559 
4560 
4561 template <>
4562 Tensor<1,1>
4563 inline
4564 GeometryInfo<1>::
4565 d_linear_shape_function_gradient (const Point<1> &,
4566  const unsigned int i)
4567 {
4568  Assert (i < GeometryInfo<1>::vertices_per_cell,
4569  ExcIndexRange (i, 0, GeometryInfo<1>::vertices_per_cell));
4570 
4571  switch (i)
4572  {
4573  case 0:
4574  return Point<1>(-1.);
4575  case 1:
4576  return Point<1>(1.);
4577  }
4578 
4579  return Point<1>(-1e9);
4580 }
4581 
4582 
4583 
4584 template <>
4585 Tensor<1,2>
4586 inline
4587 GeometryInfo<2>::
4588 d_linear_shape_function_gradient (const Point<2> &xi,
4589  const unsigned int i)
4590 {
4591  Assert (i < GeometryInfo<2>::vertices_per_cell,
4592  ExcIndexRange (i, 0, GeometryInfo<2>::vertices_per_cell));
4593 
4594  const double x = xi[0];
4595  const double y = xi[1];
4596  switch (i)
4597  {
4598  case 0:
4599  return Point<2>(-(1-y),-(1-x));
4600  case 1:
4601  return Point<2>(1-y,-x);
4602  case 2:
4603  return Point<2>(-y, 1-x);
4604  case 3:
4605  return Point<2>(y,x);
4606  }
4607  return Point<2> (-1e9, -1e9);
4608 }
4609 
4610 
4611 
4612 template <>
4613 Tensor<1,3>
4614 inline
4615 GeometryInfo<3>::
4616 d_linear_shape_function_gradient (const Point<3> &xi,
4617  const unsigned int i)
4618 {
4619  Assert (i < GeometryInfo<3>::vertices_per_cell,
4620  ExcIndexRange (i, 0, GeometryInfo<3>::vertices_per_cell));
4621 
4622  const double x = xi[0];
4623  const double y = xi[1];
4624  const double z = xi[2];
4625  switch (i)
4626  {
4627  case 0:
4628  return Point<3>(-(1-y)*(1-z),
4629  -(1-x)*(1-z),
4630  -(1-x)*(1-y));
4631  case 1:
4632  return Point<3>((1-y)*(1-z),
4633  -x*(1-z),
4634  -x*(1-y));
4635  case 2:
4636  return Point<3>(-y*(1-z),
4637  (1-x)*(1-z),
4638  -(1-x)*y);
4639  case 3:
4640  return Point<3>(y*(1-z),
4641  x*(1-z),
4642  -x*y);
4643  case 4:
4644  return Point<3>(-(1-y)*z,
4645  -(1-x)*z,
4646  (1-x)*(1-y));
4647  case 5:
4648  return Point<3>((1-y)*z,
4649  -x*z,
4650  x*(1-y));
4651  case 6:
4652  return Point<3>(-y*z,
4653  (1-x)*z,
4654  (1-x)*y);
4655  case 7:
4656  return Point<3>(y*z, x*z, x*y);
4657  }
4658 
4659  return Point<3> (-1e9, -1e9, -1e9);
4660 }
4661 
4662 
4663 
4664 template <int dim>
4665 inline
4666 Tensor<1,dim>
4667 GeometryInfo<dim>::
4668 d_linear_shape_function_gradient (const Point<dim> &,
4669  const unsigned int)
4670 {
4671  Assert (false, ExcNotImplemented());
4672  return Tensor<1,dim>();
4673 }
4674 
4675 
4676 unsigned int
4677 inline
4678 GeometryInfo<0>::n_children(const RefinementCase<0> &)
4679 {
4680  return 0;
4681 }
4682 
4683 
4684 namespace internal
4685 {
4686  namespace GeometryInfoHelper
4687  {
4688  // wedge product of a single
4689  // vector in 2d: we just have to
4690  // rotate it by 90 degrees to the
4691  // right
4692  inline
4693  Tensor<1,2>
4694  wedge_product (const Tensor<1,2> (&derivative)[1])
4695  {
4696  Tensor<1,2> result;
4697  result[0] = derivative[0][1];
4698  result[1] = -derivative[0][0];
4699 
4700  return result;
4701  }
4702 
4703 
4704  // wedge product of 2 vectors in
4705  // 3d is the cross product
4706  inline
4707  Tensor<1,3>
4708  wedge_product (const Tensor<1,3> (&derivative)[2])
4709  {
4710  return cross_product_3d (derivative[0], derivative[1]);
4711  }
4712 
4713 
4714  // wedge product of dim vectors
4715  // in dim-d: that's the
4716  // determinant of the matrix
4717  template <int dim>
4718  inline
4719  Tensor<0,dim>
4720  wedge_product (const Tensor<1,dim> (&derivative)[dim])
4721  {
4722  Tensor<2,dim> jacobian;
4723  for (unsigned int i=0; i<dim; ++i)
4724  jacobian[i] = derivative[i];
4725 
4726  return determinant (jacobian);
4727  }
4728  }
4729 }
4730 
4731 
4732 template <int dim>
4733 template <int spacedim>
4734 inline
4735 void
4736 GeometryInfo<dim>::
4737 alternating_form_at_vertices
4738 #ifndef DEAL_II_CONSTEXPR_BUG
4739 (const Point<spacedim> (&vertices)[vertices_per_cell],
4740  Tensor<spacedim-dim,spacedim> (&forms)[vertices_per_cell])
4741 #else
4742 (const Point<spacedim> *vertices,
4743  Tensor<spacedim-dim,spacedim> *forms)
4744 #endif
4745 {
4746  // for each of the vertices,
4747  // compute the alternating form
4748  // of the mapped unit
4749  // vectors. consider for
4750  // example the case of a quad
4751  // in spacedim==3: to do so, we
4752  // need to see how the
4753  // infinitesimal vectors
4754  // (d\xi_1,0) and (0,d\xi_2)
4755  // are transformed into
4756  // spacedim-dimensional space
4757  // and then form their cross
4758  // product (i.e. the wedge product
4759  // of two vectors). to this end, note
4760  // that
4761  // \vec x = sum_i \vec v_i phi_i(\vec xi)
4762  // so the transformed vectors are
4763  // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
4764  // and
4765  // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
4766  // which boils down to the columns
4767  // of the 3x2 matrix \grad_\xi x(\xi)
4768  //
4769  // a similar reasoning would
4770  // hold for all dim,spacedim
4771  // pairs -- we only have to
4772  // compute the wedge product of
4773  // the columns of the
4774  // derivatives
4775  for (unsigned int i=0; i<vertices_per_cell; ++i)
4776  {
4777  Tensor<1,spacedim> derivatives[dim];
4778 
4779  for (unsigned int j=0; j<vertices_per_cell; ++j)
4780  {
4781  const Tensor<1,dim> grad_phi_j
4782  = d_linear_shape_function_gradient (unit_cell_vertex(i),
4783  j);
4784  for (unsigned int l=0; l<dim; ++l)
4785  derivatives[l] += vertices[j] * grad_phi_j[l];
4786  }
4787 
4788  forms[i] = internal::GeometryInfoHelper::wedge_product (derivatives);
4789  }
4790 }
4791 
4792 #endif // DOXYGEN
4793 DEAL_II_NAMESPACE_CLOSE
4794 
4795 #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::vertices_per_face
static constexpr unsigned int vertices_per_face
Definition: geometry_info.h:1709
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:2617
GeometryInfo::child_cell_from_point
static unsigned int child_cell_from_point(const Point< dim > &p)
GeometryInfo::ucd_to_deal
static constexpr std::array< unsigned int, vertices_per_cell > ucd_to_deal
Definition: geometry_info.h:1772
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:1699
GeometryInfo::max_children_per_cell
static constexpr unsigned int max_children_per_cell
Definition: geometry_info.h:1685
GeometryInfo::is_inside_unit_cell
static bool is_inside_unit_cell(const Point< dim > &p)
GeometryPrimitive::GeometryPrimitive
GeometryPrimitive(const Object object)
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double &tolerance=1e-10)
Definition: utilities.cc:183
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::quads_per_face
static constexpr unsigned int quads_per_face
Definition: geometry_info.h:1721
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::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::unit_cell_vertex
static Point< dim > unit_cell_vertex(const unsigned int vertex)
internal::SubfaceCase::memory_consumption
static constexpr std::size_t memory_consumption()
StandardExceptions::ExcIndexRange
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
internal::SubfaceCase::ExcInvalidSubfaceCase
static ::ExceptionBase & ExcInvalidSubfaceCase(int arg1)
GeometryInfo< 0 >::ucd_to_deal
static const std::array< unsigned int, vertices_per_cell > ucd_to_deal
Definition: geometry_info.h:1128
GeometryInfo::vertices_per_cell
static constexpr unsigned int vertices_per_cell
Definition: geometry_info.h:1704
RefinementCase::operator~
RefinementCase operator~() const
GeometryInfo::dx_to_deal
static constexpr std::array< unsigned int, vertices_per_cell > dx_to_deal
Definition: geometry_info.h:1788
GeometryInfo::d_linear_shape_function_gradient
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
GeometryInfo::distance_to_unit_cell
static double distance_to_unit_cell(const Point< dim > &p)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
StandardExceptions::ExcImpossibleInDim
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
GeometryInfo::ExcInvalidSubface
static ::ExceptionBase & ExcInvalidSubface(int arg1, int arg2, int arg3)
RefinementCase::operator|
RefinementCase operator|(const RefinementCase &r) const
GeometryInfo::unit_normal_direction
static constexpr std::array< unsigned int, faces_per_cell > unit_normal_direction
Definition: geometry_info.h:2242
DeclException1
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:346
GeometryInfo::lines_per_cell
static constexpr unsigned int lines_per_cell
Definition: geometry_info.h:1733
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:1142
GeometryInfo::faces_per_cell
static constexpr unsigned int faces_per_cell
Definition: geometry_info.h:1690
GeometryInfo::hexes_per_cell
static constexpr unsigned int hexes_per_cell
Definition: geometry_info.h:1751
GeometryInfo::ExcInvalidCoordinate
static ::ExceptionBase & ExcInvalidCoordinate(double arg1)
RefinementCase::operator&
RefinementCase operator&(const RefinementCase &r) const
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::alternating_form_at_vertices
static void alternating_form_at_vertices(const Point< spacedim >(&vertices)[vertices_per_cell], Tensor< spacedim-dim, spacedim >(&forms)[vertices_per_cell])
GeometryInfo::project_to_unit_cell
static Point< dim > project_to_unit_cell(const Point< dim > &p)
internal::SubfaceCase::SubfaceCase
SubfaceCase(const typename SubfacePossibilities< dim >::Possibilities subface_possibility)
GeometryInfo::quads_per_cell
static constexpr unsigned int quads_per_cell
Definition: geometry_info.h:1744
GeometryInfo< 0 >::dx_to_deal
static const std::array< unsigned int, vertices_per_cell > dx_to_deal
Definition: geometry_info.h:1143
GeometryInfo::unit_normal_orientation
static constexpr std::array< int, faces_per_cell > unit_normal_orientation
Definition: geometry_info.h:2261
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:1715
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:685
DeclException3
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:370
GeometryInfo::opposite_face
static constexpr std::array< unsigned int, faces_per_cell > opposite_face
Definition: geometry_info.h:2269
RefinementCase::ExcInvalidRefinementCase
static ::ExceptionBase & ExcInvalidRefinementCase(int arg1)
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)
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
GeometryInfo::vertex_to_face
static constexpr std::array< std::array< unsigned int, dim >, vertices_per_cell > vertex_to_face
Definition: geometry_info.h:1801
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
Generated by   doxygen 1.8.14