Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
geometry_info.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_geometry_info_h
17 #define dealii_geometry_info_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/point.h>
24 #include <deal.II/base/std_cxx20/iota_view.h>
25 
26 #include <array>
27 #include <cstdint>
28 
29 
30 
31 DEAL_II_NAMESPACE_OPEN
32 
33 #ifndef DOXYGEN
34 namespace internal
35 {
36  namespace GeometryInfoHelper
37  {
38  // A struct that holds the values for all the arrays we want to initialize
39  // in GeometryInfo
40  template <int dim>
41  struct Initializers;
42 
43  template <>
44  struct Initializers<1>
45  {
46  static constexpr std::array<unsigned int, 2>
47  ucd_to_deal()
48  {
49  return {{0, 1}};
50  }
51 
52  static constexpr std::array<unsigned int, 2>
53  unit_normal_direction()
54  {
55  return {{0, 0}};
56  }
57 
58  static constexpr std::array<int, 2>
59  unit_normal_orientation()
60  {
61  return {{-1, 1}};
62  }
63 
64  static constexpr std::array<Tensor<1, 1>, 2>
65  unit_normal_vector()
66  {
67  return {{Tensor<1, 1>{{-1}}, Tensor<1, 1>{{1}}}};
68  }
69 
70  static constexpr std::array<std::array<Tensor<1, 1>, 0>, 2>
71  unit_tangential_vectors()
72  {
73  return {{{{}}, {{}}}};
74  }
75 
76  static constexpr std::array<unsigned int, 2>
77  opposite_face()
78  {
79  return {{1, 0}};
80  }
81 
82  static constexpr std::array<unsigned int, 2>
83  dx_to_deal()
84  {
85  return {{0, 1}};
86  }
87 
88  static constexpr std::array<std::array<unsigned int, 1>, 2>
89  vertex_to_face()
90  {
91  return {{{{0}}, {{1}}}};
92  }
93  };
94 
95  template <>
96  struct Initializers<2>
97  {
98  static constexpr std::array<unsigned int, 4>
99  ucd_to_deal()
100  {
101  return {{0, 1, 3, 2}};
102  }
103 
104  static constexpr std::array<unsigned int, 4>
105  unit_normal_direction()
106  {
107  return {{0, 0, 1, 1}};
108  }
109 
110  static constexpr std::array<int, 4>
111  unit_normal_orientation()
112  {
113  return {{-1, 1, -1, 1}};
114  }
115 
116  static constexpr std::array<Tensor<1, 2>, 4>
117  unit_normal_vector()
118  {
119  return {{Tensor<1, 2>{{-1., 0.}},
120  Tensor<1, 2>{{1., 0.}},
121  Tensor<1, 2>{{0., -1.}},
122  Tensor<1, 2>{{0., 1.}}}};
123  }
124 
125  static constexpr std::array<std::array<Tensor<1, 2>, 1>, 4>
126  unit_tangential_vectors()
127  {
128  return {{{{Tensor<1, 2>{{0, -1}}}},
129  {{Tensor<1, 2>{{0, 1}}}},
130  {{Tensor<1, 2>{{1, 0}}}},
131  {{Tensor<1, 2>{{-1, 0}}}}}};
132  }
133 
134  static constexpr std::array<unsigned int, 4>
135  opposite_face()
136  {
137  return {{1, 0, 3, 2}};
138  }
139 
140  static constexpr std::array<unsigned int, 4>
141  dx_to_deal()
142  {
143  return {{0, 2, 1, 3}};
144  }
145 
146  static constexpr std::array<std::array<unsigned int, 2>, 4>
147  vertex_to_face()
148  {
149  return {{{{0, 2}}, {{1, 2}}, {{0, 3}}, {{1, 3}}}};
150  }
151  };
152 
153  template <>
154  struct Initializers<3>
155  {
156  static constexpr std::array<unsigned int, 8>
157  ucd_to_deal()
158  {
159  return {{0, 1, 5, 4, 2, 3, 7, 6}};
160  }
161 
162  static constexpr std::array<unsigned int, 6>
163  unit_normal_direction()
164  {
165  return {{0, 0, 1, 1, 2, 2}};
166  }
167 
168  static constexpr std::array<int, 6>
169  unit_normal_orientation()
170  {
171  return {{-1, 1, -1, 1, -1, 1}};
172  }
173 
174  static constexpr std::array<Tensor<1, 3>, 6>
175  unit_normal_vector()
176  {
177  return {{Tensor<1, 3>{{-1, 0, 0}},
178  Tensor<1, 3>{{1, 0, 0}},
179  Tensor<1, 3>{{0, -1, 0}},
180  Tensor<1, 3>{{0, 1, 0}},
181  Tensor<1, 3>{{0, 0, -1}},
182  Tensor<1, 3>{{0, 0, 1}}}};
183  }
184 
185  static constexpr std::array<std::array<Tensor<1, 3>, 2>, 6>
186  unit_tangential_vectors()
187  {
188  return {{{{Tensor<1, 3>{{0, -1, 0}}, Tensor<1, 3>{{0, 0, 1}}}},
189  {{Tensor<1, 3>{{0, 1, 0}}, Tensor<1, 3>{{0, 0, 1}}}},
190  {{Tensor<1, 3>{{0, 0, -1}}, Tensor<1, 3>{{1, 0, 0}}}},
191  {{Tensor<1, 3>{{0, 0, 1}}, Tensor<1, 3>{{1, 0, 0}}}},
192  {{Tensor<1, 3>{{-1, 0, 0}}, Tensor<1, 3>{{0, 1, 0}}}},
193  {{Tensor<1, 3>{{1, 0, 0}}, Tensor<1, 3>{{0, 1, 0}}}}}};
194  }
195 
196  static constexpr std::array<unsigned int, 6>
197  opposite_face()
198  {
199  return {{1, 0, 3, 2, 5, 4}};
200  }
201 
202  static constexpr std::array<unsigned int, 8>
203  dx_to_deal()
204  {
205  return {{0, 4, 2, 6, 1, 5, 3, 7}};
206  }
207 
208  static constexpr std::array<std::array<unsigned int, 3>, 8>
209  vertex_to_face()
210  {
211  return {{{{0, 2, 4}},
212  {{1, 2, 4}},
213  {{0, 3, 4}},
214  {{1, 3, 4}},
215  {{0, 2, 5}},
216  {{1, 2, 5}},
217  {{0, 3, 5}},
218  {{1, 3, 5}}}};
219  }
220  };
221 
222  template <>
223  struct Initializers<4>
224  {
225  static constexpr std::array<unsigned int, 16>
226  ucd_to_deal()
227  {
228  return {{numbers::invalid_unsigned_int,
229  numbers::invalid_unsigned_int,
230  numbers::invalid_unsigned_int,
231  numbers::invalid_unsigned_int,
232  numbers::invalid_unsigned_int,
233  numbers::invalid_unsigned_int,
234  numbers::invalid_unsigned_int,
235  numbers::invalid_unsigned_int,
236  numbers::invalid_unsigned_int,
237  numbers::invalid_unsigned_int,
238  numbers::invalid_unsigned_int,
239  numbers::invalid_unsigned_int,
240  numbers::invalid_unsigned_int,
241  numbers::invalid_unsigned_int,
242  numbers::invalid_unsigned_int,
243  numbers::invalid_unsigned_int}};
244  }
245 
246  static constexpr std::array<unsigned int, 8>
247  unit_normal_direction()
248  {
249  return {{0, 0, 1, 1, 2, 2, 3, 3}};
250  }
251 
252  static constexpr std::array<int, 8>
253  unit_normal_orientation()
254  {
255  return {{-1, 1, -1, 1, -1, 1, -1, 1}};
256  }
257 
258  static constexpr std::array<Tensor<1, 4>, 8>
259  unit_normal_vector()
260  {
261  return {{Tensor<1, 4>{{-1, 0, 0, 0}},
262  Tensor<1, 4>{{1, 0, 0, 0}},
263  Tensor<1, 4>{{0, -1, 0, 0}},
264  Tensor<1, 4>{{0, 1, 0, 0}},
265  Tensor<1, 4>{{0, 0, -1, 0}},
266  Tensor<1, 4>{{0, 0, 1, 0}},
267  Tensor<1, 4>{{0, 0, 0, -1}},
268  Tensor<1, 4>{{0, 0, 0, 1}}}};
269  }
270 
271  static constexpr std::array<std::array<Tensor<1, 4>, 3>, 8>
272  unit_tangential_vectors()
273  {
274  return {{{{Tensor<1, 4>{{0, -1, 0, 0}},
275  Tensor<1, 4>{{0, 0, 1, 0}},
276  Tensor<1, 4>{{0, 0, 0, 1}}}},
277  {{Tensor<1, 4>{{0, 1, 0, 0}},
278  Tensor<1, 4>{{0, 0, 1, 0}},
279  Tensor<1, 4>{{0, 0, 0, 1}}}},
280  {{Tensor<1, 4>{{0, 0, -1, 0}},
281  Tensor<1, 4>{{0, 0, 0, 1}},
282  Tensor<1, 4>{{1, 0, 0, 0}}}},
283  {{Tensor<1, 4>{{0, 0, 1, 0}},
284  Tensor<1, 4>{{0, 0, 0, 1}},
285  Tensor<1, 4>{{1, 0, 0, 0}}}},
286  {{Tensor<1, 4>{{0, 0, 0, -1}},
287  Tensor<1, 4>{{1, 0, 0, 0}},
288  Tensor<1, 4>{{0, 1, 0, 0}}}},
289  {{Tensor<1, 4>{{0, 0, 0, 1}},
290  Tensor<1, 4>{{1, 0, 0, 0}},
291  Tensor<1, 4>{{0, 1, 0, 0}}}},
292  {{Tensor<1, 4>{{-1, 0, 0, 0}},
293  Tensor<1, 4>{{0, 1, 0, 0}},
294  Tensor<1, 4>{{0, 0, 1, 0}}}},
295  {{Tensor<1, 4>{{1, 0, 0, 0}},
296  Tensor<1, 4>{{0, 1, 0, 0}},
297  Tensor<1, 4>{{0, 0, 1, 0}}}}}};
298  }
299 
300  static constexpr std::array<unsigned int, 8>
301  opposite_face()
302  {
303  return {{1, 0, 3, 2, 5, 4, 7, 6}};
304  }
305 
306  static constexpr std::array<unsigned int, 16>
307  dx_to_deal()
308  {
309  return {{numbers::invalid_unsigned_int,
310  numbers::invalid_unsigned_int,
311  numbers::invalid_unsigned_int,
312  numbers::invalid_unsigned_int,
313  numbers::invalid_unsigned_int,
314  numbers::invalid_unsigned_int,
315  numbers::invalid_unsigned_int,
316  numbers::invalid_unsigned_int,
317  numbers::invalid_unsigned_int,
318  numbers::invalid_unsigned_int,
319  numbers::invalid_unsigned_int,
320  numbers::invalid_unsigned_int,
321  numbers::invalid_unsigned_int,
322  numbers::invalid_unsigned_int,
323  numbers::invalid_unsigned_int,
324  numbers::invalid_unsigned_int}};
325  }
326 
327  static constexpr std::array<std::array<unsigned int, 4>, 16>
328  vertex_to_face()
329  {
330  return {{{{numbers::invalid_unsigned_int,
331  numbers::invalid_unsigned_int,
332  numbers::invalid_unsigned_int,
333  numbers::invalid_unsigned_int}},
334  {{numbers::invalid_unsigned_int,
335  numbers::invalid_unsigned_int,
336  numbers::invalid_unsigned_int,
337  numbers::invalid_unsigned_int}},
338  {{numbers::invalid_unsigned_int,
339  numbers::invalid_unsigned_int,
340  numbers::invalid_unsigned_int,
341  numbers::invalid_unsigned_int}},
342  {{numbers::invalid_unsigned_int,
343  numbers::invalid_unsigned_int,
344  numbers::invalid_unsigned_int,
345  numbers::invalid_unsigned_int}},
346  {{numbers::invalid_unsigned_int,
347  numbers::invalid_unsigned_int,
348  numbers::invalid_unsigned_int,
349  numbers::invalid_unsigned_int}},
350  {{numbers::invalid_unsigned_int,
351  numbers::invalid_unsigned_int,
352  numbers::invalid_unsigned_int,
353  numbers::invalid_unsigned_int}},
354  {{numbers::invalid_unsigned_int,
355  numbers::invalid_unsigned_int,
356  numbers::invalid_unsigned_int,
357  numbers::invalid_unsigned_int}},
358  {{numbers::invalid_unsigned_int,
359  numbers::invalid_unsigned_int,
360  numbers::invalid_unsigned_int,
361  numbers::invalid_unsigned_int}},
362  {{numbers::invalid_unsigned_int,
363  numbers::invalid_unsigned_int,
364  numbers::invalid_unsigned_int,
365  numbers::invalid_unsigned_int}},
366  {{numbers::invalid_unsigned_int,
367  numbers::invalid_unsigned_int,
368  numbers::invalid_unsigned_int,
369  numbers::invalid_unsigned_int}},
370  {{numbers::invalid_unsigned_int,
371  numbers::invalid_unsigned_int,
372  numbers::invalid_unsigned_int,
373  numbers::invalid_unsigned_int}},
374  {{numbers::invalid_unsigned_int,
375  numbers::invalid_unsigned_int,
376  numbers::invalid_unsigned_int,
377  numbers::invalid_unsigned_int}},
378  {{numbers::invalid_unsigned_int,
379  numbers::invalid_unsigned_int,
380  numbers::invalid_unsigned_int,
381  numbers::invalid_unsigned_int}},
382  {{numbers::invalid_unsigned_int,
383  numbers::invalid_unsigned_int,
384  numbers::invalid_unsigned_int,
385  numbers::invalid_unsigned_int}},
386  {{numbers::invalid_unsigned_int,
387  numbers::invalid_unsigned_int,
388  numbers::invalid_unsigned_int,
389  numbers::invalid_unsigned_int}},
390  {{numbers::invalid_unsigned_int,
391  numbers::invalid_unsigned_int,
392  numbers::invalid_unsigned_int,
393  numbers::invalid_unsigned_int}}}};
394  }
395  };
396  } // namespace GeometryInfoHelper
397 } // namespace internal
398 #endif // DOXYGEN
399 
400 
418 class GeometryPrimitive
419 {
420 public:
427  enum Object
428  {
432  vertex = 0,
436  line = 1,
440  quad = 2,
444  hex = 3
445  };
446 
451  GeometryPrimitive(const Object object);
452 
458  GeometryPrimitive(const unsigned int object_dimension);
459 
464  operator unsigned int() const;
465 
466 private:
470  Object object;
471 };
472 
473 
474 
488 template <int dim>
489 struct RefinementPossibilities
490 {
525  enum Possibilities
526  {
530  no_refinement = 0,
531 
543  isotropic_refinement = 0xFF
544  };
545 };
546 
547 
548 
559 template <>
560 struct RefinementPossibilities<1>
561 {
596  enum Possibilities
597  {
601  no_refinement = 0,
605  cut_x = 1,
609  isotropic_refinement = cut_x
610  };
611 };
612 
613 
614 
626 template <>
627 struct RefinementPossibilities<2>
628 {
663  enum Possibilities
664  {
668  no_refinement = 0,
672  cut_x = 1,
676  cut_y = 2,
680  cut_xy = cut_x | cut_y,
681 
685  isotropic_refinement = cut_xy
686  };
687 };
688 
689 
690 
702 template <>
703 struct RefinementPossibilities<3>
704 {
739  enum Possibilities
740  {
744  no_refinement = 0,
748  cut_x = 1,
752  cut_y = 2,
756  cut_xy = cut_x | cut_y,
760  cut_z = 4,
764  cut_xz = cut_x | cut_z,
768  cut_yz = cut_y | cut_z,
772  cut_xyz = cut_x | cut_y | cut_z,
773 
777  isotropic_refinement = cut_xyz
778  };
779 };
780 
781 
782 
794 template <int dim>
795 class RefinementCase : public RefinementPossibilities<dim>
796 {
797 public:
801  RefinementCase();
802 
807  RefinementCase(
808  const typename RefinementPossibilities<dim>::Possibilities refinement_case);
809 
815  explicit RefinementCase(const std::uint8_t refinement_case);
816 
828  operator std::uint8_t() const;
829 
834  RefinementCase
835  operator|(const RefinementCase &r) const;
836 
841  RefinementCase operator&(const RefinementCase &r) const;
842 
850  RefinementCase
851  operator~() const;
852 
853 
859  static RefinementCase
860  cut_axis(const unsigned int i);
861 
865  static std::size_t
866  memory_consumption();
867 
872  template <class Archive>
873  void
874  serialize(Archive &ar, const unsigned int version);
875 
879  DeclException1(
880  ExcInvalidRefinementCase,
881  int,
882  << "The refinement flags given (" << arg1
883  << ") contain set bits that do not "
884  << "make sense for the space dimension of the object to which they are applied.");
885 
886 private:
891  std::uint8_t value : (dim > 0 ? dim : 1);
892 };
893 
894 
895 namespace internal
896 {
917  template <int dim>
918  struct SubfacePossibilities
919  {
923  enum Possibilities
924  {
928  case_none = 0,
929 
933  case_isotropic = static_cast<std::uint8_t>(-1)
934  };
935  };
936 
937 
947  template <>
948  struct SubfacePossibilities<0>
949  {
955  enum Possibilities
956  {
960  case_none = 0,
961 
965  case_isotropic = case_none
966  };
967  };
968 
969 
970 
981  template <>
982  struct SubfacePossibilities<1>
983  {
989  enum Possibilities
990  {
994  case_none = 0,
995 
999  case_isotropic = case_none
1000  };
1001  };
1002 
1003 
1004 
1016  template <>
1017  struct SubfacePossibilities<2>
1018  {
1025  enum Possibilities
1026  {
1030  case_none = 0,
1034  case_x = 1,
1038  case_isotropic = case_x
1039  };
1040  };
1041 
1042 
1043 
1136  template <>
1137  struct SubfacePossibilities<3>
1138  {
1145  enum Possibilities
1146  {
1147  case_none = 0,
1148  case_x = 1,
1149  case_x1y = 2,
1150  case_x2y = 3,
1151  case_x1y2y = 4,
1152  case_y = 5,
1153  case_y1x = 6,
1154  case_y2x = 7,
1155  case_y1x2x = 8,
1156  case_xy = 9,
1157 
1158  case_isotropic = case_xy
1159  };
1160  };
1161 
1162 
1163 
1171  template <int dim>
1172  class SubfaceCase : public SubfacePossibilities<dim>
1173  {
1174  public:
1180  SubfaceCase(const typename SubfacePossibilities<dim>::Possibilities
1181  subface_possibility);
1182 
1194  operator std::uint8_t() const;
1195 
1199  static constexpr std::size_t
1200  memory_consumption();
1201 
1205  DeclException1(
1206  ExcInvalidSubfaceCase,
1207  int,
1208  << "The subface case given (" << arg1 << ") does not make sense "
1209  << "for the space dimension of the object to which they are applied.");
1210 
1211  private:
1216  std::uint8_t value : (dim == 3 ? 4 : 1);
1217  };
1218 
1219 } // namespace internal
1220 
1221 
1222 
1223 template <int dim>
1224 struct GeometryInfo;
1225 
1226 
1227 
1245 template <>
1246 struct GeometryInfo<0>
1247 {
1255  static constexpr unsigned int max_children_per_cell = 1;
1256 
1260  static constexpr unsigned int faces_per_cell = 0;
1261 
1278  static std::array<unsigned int, 0>
1279  face_indices();
1280 
1288  static constexpr unsigned int max_children_per_face = 0;
1289 
1295  static unsigned int
1296  n_children(const RefinementCase<0> &refinement_case);
1297 
1301  static constexpr unsigned int vertices_per_cell = 1;
1302 
1321  static std::array<unsigned int, vertices_per_cell>
1322  vertex_indices();
1323 
1330  static constexpr unsigned int vertices_per_face = 0;
1331 
1335  static constexpr unsigned int lines_per_face = 0;
1336 
1340  static constexpr unsigned int quads_per_face = 0;
1341 
1345  static constexpr unsigned int lines_per_cell = 0;
1346 
1350  static constexpr unsigned int quads_per_cell = 0;
1351 
1355  static constexpr unsigned int hexes_per_cell = 0;
1356 
1374  static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1375 
1389  static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1390 };
1391 
1392 
1393 
1918 template <int dim>
1919 struct GeometryInfo
1920 {
1928  static constexpr unsigned int max_children_per_cell = 1 << dim;
1929 
1933  static constexpr unsigned int faces_per_cell = 2 * dim;
1934 
1951  static std_cxx20::ranges::iota_view<unsigned int, unsigned int>
1952  face_indices();
1953 
1961  static constexpr unsigned int max_children_per_face =
1962  GeometryInfo<dim - 1>::max_children_per_cell;
1963 
1967  static constexpr unsigned int vertices_per_cell = 1 << dim;
1968 
1984  static std_cxx20::ranges::iota_view<unsigned int, unsigned int>
1985  vertex_indices();
1986 
1990  static constexpr unsigned int vertices_per_face =
1991  GeometryInfo<dim - 1>::vertices_per_cell;
1992 
1996  static constexpr unsigned int lines_per_face =
1997  GeometryInfo<dim - 1>::lines_per_cell;
1998 
2002  static constexpr unsigned int quads_per_face =
2003  GeometryInfo<dim - 1>::quads_per_cell;
2004 
2014  static constexpr unsigned int lines_per_cell =
2015  (2 * GeometryInfo<dim - 1>::lines_per_cell +
2016  GeometryInfo<dim - 1>::vertices_per_cell);
2017 
2025  static constexpr unsigned int quads_per_cell =
2026  (2 * GeometryInfo<dim - 1>::quads_per_cell +
2027  GeometryInfo<dim - 1>::lines_per_cell);
2028 
2032  static constexpr unsigned int hexes_per_cell =
2033  (2 * GeometryInfo<dim - 1>::hexes_per_cell +
2034  GeometryInfo<dim - 1>::quads_per_cell);
2035 
2053  static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
2054  internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal();
2055 
2069  static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
2070  internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal();
2071 
2082  static constexpr std::array<std::array<unsigned int, dim>, vertices_per_cell>
2083  vertex_to_face =
2084  internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face();
2085 
2090  static unsigned int
2091  n_children(const RefinementCase<dim> &refinement_case);
2092 
2097  static unsigned int
2098  n_subfaces(const internal::SubfaceCase<dim> &subface_case);
2099 
2109  static double
2110  subface_ratio(const internal::SubfaceCase<dim> &subface_case,
2111  const unsigned int subface_no);
2112 
2118  static RefinementCase<dim - 1>
2119  face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2120  const unsigned int face_no,
2121  const bool face_orientation = true,
2122  const bool face_flip = false,
2123  const bool face_rotation = false);
2124 
2130  static RefinementCase<dim>
2131  min_cell_refinement_case_for_face_refinement(
2132  const RefinementCase<dim - 1> &face_refinement_case,
2133  const unsigned int face_no,
2134  const bool face_orientation = true,
2135  const bool face_flip = false,
2136  const bool face_rotation = false);
2137 
2142  static RefinementCase<1>
2143  line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2144  const unsigned int line_no);
2145 
2150  static RefinementCase<dim>
2151  min_cell_refinement_case_for_line_refinement(const unsigned int line_no);
2152 
2199  static unsigned int
2200  child_cell_on_face(const RefinementCase<dim> & ref_case,
2201  const unsigned int face,
2202  const unsigned int subface,
2203  const bool face_orientation = true,
2204  const bool face_flip = false,
2205  const bool face_rotation = false,
2206  const RefinementCase<dim - 1> &face_refinement_case =
2207  RefinementCase<dim - 1>::isotropic_refinement);
2208 
2222  static unsigned int
2223  line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
2224 
2245  static unsigned int
2246  face_to_cell_vertices(const unsigned int face,
2247  const unsigned int vertex,
2248  const bool face_orientation = true,
2249  const bool face_flip = false,
2250  const bool face_rotation = false);
2251 
2263  static unsigned int
2264  face_to_cell_lines(const unsigned int face,
2265  const unsigned int line,
2266  const bool face_orientation = true,
2267  const bool face_flip = false,
2268  const bool face_rotation = false);
2269 
2279  static unsigned int
2280  standard_to_real_face_vertex(const unsigned int vertex,
2281  const bool face_orientation = true,
2282  const bool face_flip = false,
2283  const bool face_rotation = false);
2284 
2294  static unsigned int
2295  real_to_standard_face_vertex(const unsigned int vertex,
2296  const bool face_orientation = true,
2297  const bool face_flip = false,
2298  const bool face_rotation = false);
2299 
2309  static unsigned int
2310  standard_to_real_face_line(const unsigned int line,
2311  const bool face_orientation = true,
2312  const bool face_flip = false,
2313  const bool face_rotation = false);
2314 
2324  static unsigned int
2325  real_to_standard_face_line(const unsigned int line,
2326  const bool face_orientation = true,
2327  const bool face_flip = false,
2328  const bool face_rotation = false);
2329 
2335  static Point<dim>
2336  unit_cell_vertex(const unsigned int vertex);
2337 
2347  static unsigned int
2348  child_cell_from_point(const Point<dim> &p);
2349 
2357  static Point<dim>
2358  cell_to_child_coordinates(const Point<dim> & p,
2359  const unsigned int child_index,
2360  const RefinementCase<dim> refine_case =
2361  RefinementCase<dim>::isotropic_refinement);
2362 
2368  static Point<dim>
2369  child_to_cell_coordinates(const Point<dim> & p,
2370  const unsigned int child_index,
2371  const RefinementCase<dim> refine_case =
2372  RefinementCase<dim>::isotropic_refinement);
2373 
2378  static bool
2379  is_inside_unit_cell(const Point<dim> &p);
2380 
2392  static bool
2393  is_inside_unit_cell(const Point<dim> &p, const double eps);
2394 
2399  static Point<dim>
2400  project_to_unit_cell(const Point<dim> &p);
2401 
2407  static double
2408  distance_to_unit_cell(const Point<dim> &p);
2409 
2414  static double
2415  d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2416 
2421  static Tensor<1, dim>
2422  d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2423 
2475  template <int spacedim>
2476  static void
2477  alternating_form_at_vertices
2478 #ifndef DEAL_II_CONSTEXPR_BUG
2479  (const Point<spacedim> (&vertices)[vertices_per_cell],
2480  Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2481 #else
2482  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2483 #endif
2484 
2494  static constexpr std::array<unsigned int, faces_per_cell>
2495  unit_normal_direction =
2496  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction();
2497 
2514  static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2515  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation();
2516 
2527  static constexpr std::array<Tensor<1, dim>, faces_per_cell>
2528  unit_normal_vector =
2529  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_vector();
2530 
2544  static constexpr std::array<std::array<Tensor<1, dim>, dim - 1>,
2545  faces_per_cell>
2546  unit_tangential_vectors = internal::GeometryInfoHelper::Initializers<
2547  dim>::unit_tangential_vectors();
2548 
2554  static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2555  internal::GeometryInfoHelper::Initializers<dim>::opposite_face();
2556 
2557 
2561  DeclException1(ExcInvalidCoordinate,
2562  double,
2563  << "The coordinates must satisfy 0 <= x_i <= 1, "
2564  << "but here we have x_i=" << arg1);
2565 
2569  DeclException3(ExcInvalidSubface,
2570  int,
2571  int,
2572  int,
2573  << "RefinementCase<dim> " << arg1 << ": face " << arg2
2574  << " has no subface " << arg3);
2575 };
2576 
2577 
2578 
2579 #ifndef DOXYGEN
2580 
2581 
2582 /* -------------- declaration of explicit specializations ------------- */
2583 
2584 
2585 template <>
2586 Tensor<1, 1>
2587 GeometryInfo<1>::d_linear_shape_function_gradient(const Point<1> & xi,
2588  const unsigned int i);
2589 template <>
2590 Tensor<1, 2>
2591 GeometryInfo<2>::d_linear_shape_function_gradient(const Point<2> & xi,
2592  const unsigned int i);
2593 template <>
2594 Tensor<1, 3>
2595 GeometryInfo<3>::d_linear_shape_function_gradient(const Point<3> & xi,
2596  const unsigned int i);
2597 
2598 
2599 
2600 /* -------------- inline functions ------------- */
2601 
2602 
2603 inline GeometryPrimitive::GeometryPrimitive(const Object object)
2604  : object(object)
2605 {}
2606 
2607 
2608 
2609 inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2610  : object(static_cast<Object>(object_dimension))
2611 {}
2612 
2613 
2614 inline GeometryPrimitive::operator unsigned int() const
2615 {
2616  return static_cast<unsigned int>(object);
2617 }
2618 
2619 
2620 
2621 namespace internal
2622 {
2623  template <int dim>
2624  inline SubfaceCase<dim>::SubfaceCase(
2625  const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2626  : value(subface_possibility)
2627  {}
2628 
2629 
2630  template <int dim>
2631  inline SubfaceCase<dim>::operator std::uint8_t() const
2632  {
2633  return value;
2634  }
2635 
2636 
2637 } // namespace internal
2638 
2639 
2640 template <int dim>
2641 inline RefinementCase<dim>
2642 RefinementCase<dim>::cut_axis(const unsigned int)
2643 {
2644  Assert(false, ExcInternalError());
2645  return static_cast<std::uint8_t>(-1);
2646 }
2647 
2648 
2649 template <>
2650 inline RefinementCase<1>
2651 RefinementCase<1>::cut_axis(const unsigned int i)
2652 {
2653  AssertIndexRange(i, 1);
2654 
2655  const RefinementCase options[1] = {RefinementPossibilities<1>::cut_x};
2656  return options[i];
2657 }
2658 
2659 
2660 
2661 template <>
2662 inline RefinementCase<2>
2663 RefinementCase<2>::cut_axis(const unsigned int i)
2664 {
2665  AssertIndexRange(i, 2);
2666 
2667  const RefinementCase options[2] = {RefinementPossibilities<2>::cut_x,
2668  RefinementPossibilities<2>::cut_y};
2669  return options[i];
2670 }
2671 
2672 
2673 
2674 template <>
2675 inline RefinementCase<3>
2676 RefinementCase<3>::cut_axis(const unsigned int i)
2677 {
2678  AssertIndexRange(i, 3);
2679 
2680  const RefinementCase options[3] = {RefinementPossibilities<3>::cut_x,
2681  RefinementPossibilities<3>::cut_y,
2682  RefinementPossibilities<3>::cut_z};
2683  return options[i];
2684 }
2685 
2686 
2687 
2688 template <int dim>
2689 inline RefinementCase<dim>::RefinementCase()
2690  : value(RefinementPossibilities<dim>::no_refinement)
2691 {}
2692 
2693 
2694 
2695 template <int dim>
2696 inline RefinementCase<dim>::RefinementCase(
2697  const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2698  : value(refinement_case)
2699 {
2700  // check that only those bits of
2701  // the given argument are set that
2702  // make sense for a given space
2703  // dimension
2704  Assert((refinement_case &
2705  RefinementPossibilities<dim>::isotropic_refinement) ==
2706  refinement_case,
2707  ExcInvalidRefinementCase(refinement_case));
2708 }
2709 
2710 
2711 
2712 template <int dim>
2713 inline RefinementCase<dim>::RefinementCase(const std::uint8_t refinement_case)
2714  : value(refinement_case)
2715 {
2716  // check that only those bits of
2717  // the given argument are set that
2718  // make sense for a given space
2719  // dimension
2720  Assert((refinement_case &
2721  RefinementPossibilities<dim>::isotropic_refinement) ==
2722  refinement_case,
2723  ExcInvalidRefinementCase(refinement_case));
2724 }
2725 
2726 
2727 
2728 template <int dim>
2729 inline RefinementCase<dim>::operator std::uint8_t() const
2730 {
2731  return value;
2732 }
2733 
2734 
2735 
2736 template <int dim>
2737 inline RefinementCase<dim>
2738 RefinementCase<dim>::operator|(const RefinementCase<dim> &r) const
2739 {
2740  return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2741 }
2742 
2743 
2744 
2745 template <int dim>
2746 inline RefinementCase<dim> RefinementCase<dim>::
2747  operator&(const RefinementCase<dim> &r) const
2748 {
2749  return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2750 }
2751 
2752 
2753 
2754 template <int dim>
2755 inline RefinementCase<dim>
2756 RefinementCase<dim>::operator~() const
2757 {
2758  return RefinementCase<dim>(static_cast<std::uint8_t>(
2759  (~value) & RefinementPossibilities<dim>::isotropic_refinement));
2760 }
2761 
2762 
2763 
2764 template <int dim>
2765 inline std::size_t
2766 RefinementCase<dim>::memory_consumption()
2767 {
2768  return sizeof(RefinementCase<dim>);
2769 }
2770 
2771 
2772 
2773 template <int dim>
2774 template <class Archive>
2775 inline void
2776 RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2777 {
2778  // serialization can't deal with bitfields, so copy from/to a full sized
2779  // std::uint8_t
2780  std::uint8_t uchar_value = value;
2781  ar & uchar_value;
2782  value = uchar_value;
2783 }
2784 
2785 
2786 
2787 template <>
2788 inline Point<1>
2789 GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2790 {
2791  AssertIndexRange(vertex, vertices_per_cell);
2792 
2793  return Point<1>(static_cast<double>(vertex));
2794 }
2795 
2796 
2797 
2798 template <>
2799 inline Point<2>
2800 GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2801 {
2802  AssertIndexRange(vertex, vertices_per_cell);
2803 
2804  return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2805 }
2806 
2807 
2808 
2809 template <>
2810 inline Point<3>
2811 GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2812 {
2813  AssertIndexRange(vertex, vertices_per_cell);
2814 
2815  return {static_cast<double>(vertex % 2),
2816  static_cast<double>(vertex / 2 % 2),
2817  static_cast<double>(vertex / 4)};
2818 }
2819 
2820 
2821 
2822 inline std::array<unsigned int, 0>
2823 GeometryInfo<0>::face_indices()
2824 {
2825  return {{}};
2826 }
2827 
2828 
2829 
2830 inline std::array<unsigned int, 1>
2831 GeometryInfo<0>::vertex_indices()
2832 {
2833  return {{0}};
2834 }
2835 
2836 
2837 
2838 template <int dim>
2839 inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
2840 GeometryInfo<dim>::face_indices()
2841 {
2842  return {0U, faces_per_cell};
2843 }
2844 
2845 
2846 
2847 template <int dim>
2848 inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
2849 GeometryInfo<dim>::vertex_indices()
2850 {
2851  return {0U, vertices_per_cell};
2852 }
2853 
2854 
2855 
2856 template <int dim>
2857 inline Point<dim>
2858 GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2859 {
2860  Assert(false, ExcNotImplemented());
2861 
2862  return Point<dim>();
2863 }
2864 
2865 
2866 
2867 template <>
2868 inline unsigned int
2869 GeometryInfo<1>::child_cell_from_point(const Point<1> &p)
2870 {
2871  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2872 
2873  return (p[0] <= 0.5 ? 0 : 1);
2874 }
2875 
2876 
2877 
2878 template <>
2879 inline unsigned int
2880 GeometryInfo<2>::child_cell_from_point(const Point<2> &p)
2881 {
2882  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2883  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2884 
2885  return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
2886 }
2887 
2888 
2889 
2890 template <>
2891 inline unsigned int
2892 GeometryInfo<3>::child_cell_from_point(const Point<3> &p)
2893 {
2894  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2895  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2896  Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
2897 
2898  return (p[0] <= 0.5 ?
2899  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
2900  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
2901 }
2902 
2903 
2904 template <int dim>
2905 inline unsigned int
2906 GeometryInfo<dim>::child_cell_from_point(const Point<dim> &)
2907 {
2908  Assert(false, ExcNotImplemented());
2909 
2910  return 0;
2911 }
2912 
2913 
2914 
2915 template <>
2916 inline Point<1>
2917 GeometryInfo<1>::cell_to_child_coordinates(const Point<1> & p,
2918  const unsigned int child_index,
2919  const RefinementCase<1> refine_case)
2920 
2921 {
2922  AssertIndexRange(child_index, 2);
2923  Assert(refine_case == RefinementCase<1>::cut_x, ExcInternalError());
2924  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
2925 
2926  return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
2927 }
2928 
2929 
2930 
2931 template <>
2932 inline Point<2>
2933 GeometryInfo<2>::cell_to_child_coordinates(const Point<2> & p,
2934  const unsigned int child_index,
2935  const RefinementCase<2> refine_case)
2936 
2937 {
2938  AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
2939 
2940  Point<2> point = p;
2941  switch (refine_case)
2942  {
2943  case RefinementCase<2>::cut_x:
2944  point[0] *= 2.0;
2945  if (child_index == 1)
2946  point[0] -= 1.0;
2947  break;
2948  case RefinementCase<2>::cut_y:
2949  point[1] *= 2.0;
2950  if (child_index == 1)
2951  point[1] -= 1.0;
2952  break;
2953  case RefinementCase<2>::cut_xy:
2954  point *= 2.0;
2955  point -= unit_cell_vertex(child_index);
2956  break;
2957  default:
2958  Assert(false, ExcInternalError());
2959  }
2960 
2961  return point;
2962 }
2963 
2964 
2965 
2966 template <>
2967 inline Point<3>
2968 GeometryInfo<3>::cell_to_child_coordinates(const Point<3> & p,
2969  const unsigned int child_index,
2970  const RefinementCase<3> refine_case)
2971 
2972 {
2973  AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
2974 
2975  Point<3> point = p;
2976  // there might be a cleverer way to do
2977  // this, but since this function is called
2978  // in very few places for initialization
2979  // purposes only, I don't care at the
2980  // moment
2981  switch (refine_case)
2982  {
2983  case RefinementCase<3>::cut_x:
2984  point[0] *= 2.0;
2985  if (child_index == 1)
2986  point[0] -= 1.0;
2987  break;
2988  case RefinementCase<3>::cut_y:
2989  point[1] *= 2.0;
2990  if (child_index == 1)
2991  point[1] -= 1.0;
2992  break;
2993  case RefinementCase<3>::cut_z:
2994  point[2] *= 2.0;
2995  if (child_index == 1)
2996  point[2] -= 1.0;
2997  break;
2998  case RefinementCase<3>::cut_xy:
2999  point[0] *= 2.0;
3000  point[1] *= 2.0;
3001  if (child_index % 2 == 1)
3002  point[0] -= 1.0;
3003  if (child_index / 2 == 1)
3004  point[1] -= 1.0;
3005  break;
3006  case RefinementCase<3>::cut_xz:
3007  // careful, this is slightly
3008  // different from xy and yz due to
3009  // different internal numbering of
3010  // children!
3011  point[0] *= 2.0;
3012  point[2] *= 2.0;
3013  if (child_index / 2 == 1)
3014  point[0] -= 1.0;
3015  if (child_index % 2 == 1)
3016  point[2] -= 1.0;
3017  break;
3018  case RefinementCase<3>::cut_yz:
3019  point[1] *= 2.0;
3020  point[2] *= 2.0;
3021  if (child_index % 2 == 1)
3022  point[1] -= 1.0;
3023  if (child_index / 2 == 1)
3024  point[2] -= 1.0;
3025  break;
3026  case RefinementCase<3>::cut_xyz:
3027  point *= 2.0;
3028  point -= unit_cell_vertex(child_index);
3029  break;
3030  default:
3031  Assert(false, ExcInternalError());
3032  }
3033 
3034  return point;
3035 }
3036 
3037 
3038 
3039 template <int dim>
3040 inline Point<dim>
3041 GeometryInfo<dim>::cell_to_child_coordinates(
3042  const Point<dim> & /*p*/,
3043  const unsigned int /*child_index*/,
3044  const RefinementCase<dim> /*refine_case*/)
3045 
3046 {
3047  Assert(false, ExcNotImplemented());
3048  return Point<dim>();
3049 }
3050 
3051 
3052 
3053 template <>
3054 inline Point<1>
3055 GeometryInfo<1>::child_to_cell_coordinates(const Point<1> & p,
3056  const unsigned int child_index,
3057  const RefinementCase<1> refine_case)
3058 
3059 {
3060  AssertIndexRange(child_index, 2);
3061  Assert(refine_case == RefinementCase<1>::cut_x, ExcInternalError());
3062  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3063 
3064  return (p + unit_cell_vertex(child_index)) * 0.5;
3065 }
3066 
3067 
3068 
3069 template <>
3070 inline Point<3>
3071 GeometryInfo<3>::child_to_cell_coordinates(const Point<3> & p,
3072  const unsigned int child_index,
3073  const RefinementCase<3> refine_case)
3074 
3075 {
3076  AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3077 
3078  Point<3> point = p;
3079  // there might be a cleverer way to do
3080  // this, but since this function is called
3081  // in very few places for initialization
3082  // purposes only, I don't care at the
3083  // moment
3084  switch (refine_case)
3085  {
3086  case RefinementCase<3>::cut_x:
3087  if (child_index == 1)
3088  point[0] += 1.0;
3089  point[0] *= 0.5;
3090  break;
3091  case RefinementCase<3>::cut_y:
3092  if (child_index == 1)
3093  point[1] += 1.0;
3094  point[1] *= 0.5;
3095  break;
3096  case RefinementCase<3>::cut_z:
3097  if (child_index == 1)
3098  point[2] += 1.0;
3099  point[2] *= 0.5;
3100  break;
3101  case RefinementCase<3>::cut_xy:
3102  if (child_index % 2 == 1)
3103  point[0] += 1.0;
3104  if (child_index / 2 == 1)
3105  point[1] += 1.0;
3106  point[0] *= 0.5;
3107  point[1] *= 0.5;
3108  break;
3109  case RefinementCase<3>::cut_xz:
3110  // careful, this is slightly
3111  // different from xy and yz due to
3112  // different internal numbering of
3113  // children!
3114  if (child_index / 2 == 1)
3115  point[0] += 1.0;
3116  if (child_index % 2 == 1)
3117  point[2] += 1.0;
3118  point[0] *= 0.5;
3119  point[2] *= 0.5;
3120  break;
3121  case RefinementCase<3>::cut_yz:
3122  if (child_index % 2 == 1)
3123  point[1] += 1.0;
3124  if (child_index / 2 == 1)
3125  point[2] += 1.0;
3126  point[1] *= 0.5;
3127  point[2] *= 0.5;
3128  break;
3129  case RefinementCase<3>::cut_xyz:
3130  point += unit_cell_vertex(child_index);
3131  point *= 0.5;
3132  break;
3133  default:
3134  Assert(false, ExcInternalError());
3135  }
3136 
3137  return point;
3138 }
3139 
3140 
3141 
3142 template <>
3143 inline Point<2>
3144 GeometryInfo<2>::child_to_cell_coordinates(const Point<2> & p,
3145  const unsigned int child_index,
3146  const RefinementCase<2> refine_case)
3147 {
3148  AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3149 
3150  Point<2> point = p;
3151  switch (refine_case)
3152  {
3153  case RefinementCase<2>::cut_x:
3154  if (child_index == 1)
3155  point[0] += 1.0;
3156  point[0] *= 0.5;
3157  break;
3158  case RefinementCase<2>::cut_y:
3159  if (child_index == 1)
3160  point[1] += 1.0;
3161  point[1] *= 0.5;
3162  break;
3163  case RefinementCase<2>::cut_xy:
3164  point += unit_cell_vertex(child_index);
3165  point *= 0.5;
3166  break;
3167  default:
3168  Assert(false, ExcInternalError());
3169  }
3170 
3171  return point;
3172 }
3173 
3174 
3175 
3176 template <int dim>
3177 inline Point<dim>
3178 GeometryInfo<dim>::child_to_cell_coordinates(
3179  const Point<dim> & /*p*/,
3180  const unsigned int /*child_index*/,
3181  const RefinementCase<dim> /*refine_case*/)
3182 {
3183  Assert(false, ExcNotImplemented());
3184  return Point<dim>();
3185 }
3186 
3187 
3188 
3189 template <int dim>
3190 inline bool
3191 GeometryInfo<dim>::is_inside_unit_cell(const Point<dim> &)
3192 {
3193  Assert(false, ExcNotImplemented());
3194  return false;
3195 }
3196 
3197 template <>
3198 inline bool
3199 GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p)
3200 {
3201  return (p[0] >= 0.) && (p[0] <= 1.);
3202 }
3203 
3204 
3205 
3206 template <>
3207 inline bool
3208 GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p)
3209 {
3210  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
3211 }
3212 
3213 
3214 
3215 template <>
3216 inline bool
3217 GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p)
3218 {
3219  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
3220  (p[2] >= 0.) && (p[2] <= 1.);
3221 }
3222 
3223 
3224 
3225 template <int dim>
3226 inline bool
3227 GeometryInfo<dim>::is_inside_unit_cell(const Point<dim> &, const double)
3228 {
3229  Assert(false, ExcNotImplemented());
3230  return false;
3231 }
3232 
3233 template <>
3234 inline bool
3235 GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
3236 {
3237  return (p[0] >= -eps) && (p[0] <= 1. + eps);
3238 }
3239 
3240 
3241 
3242 template <>
3243 inline bool
3244 GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
3245 {
3246  const double l = -eps, u = 1 + eps;
3247  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
3248 }
3249 
3250 
3251 
3252 template <>
3253 inline bool
3254 GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
3255 {
3256  const double l = -eps, u = 1.0 + eps;
3257  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
3258  (p[2] >= l) && (p[2] <= u);
3259 }
3260 
3261 
3262 
3263 template <>
3264 inline unsigned int
3265 GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
3266  const unsigned int vertex)
3267 {
3268  (void)line;
3269  AssertIndexRange(line, lines_per_cell);
3270  AssertIndexRange(vertex, 2);
3271 
3272  return vertex;
3273 }
3274 
3275 
3276 template <>
3277 inline unsigned int
3278 GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
3279  const unsigned int vertex)
3280 {
3281  constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
3282  return cell_vertices[line][vertex];
3283 }
3284 
3285 
3286 
3287 template <>
3288 inline unsigned int
3289 GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
3290  const unsigned int vertex)
3291 {
3292  AssertIndexRange(line, lines_per_cell);
3293  AssertIndexRange(vertex, 2);
3294 
3295  constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3296  {1, 3},
3297  {0, 1},
3298  {2, 3},
3299  {4, 6}, // top face
3300  {5, 7},
3301  {4, 5},
3302  {6, 7},
3303  {0,
3304  4}, // connects of bottom
3305  {1, 5}, // top face
3306  {2, 6},
3307  {3, 7}};
3308  return vertices[line][vertex];
3309 }
3310 
3311 
3312 template <>
3313 inline unsigned int
3314 GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3315 {
3316  Assert(false, ExcNotImplemented());
3317  return numbers::invalid_unsigned_int;
3318 }
3319 
3320 template <>
3321 inline unsigned int
3322 GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3323  const bool face_orientation,
3324  const bool face_flip,
3325  const bool face_rotation)
3326 {
3327  AssertIndexRange(vertex, GeometryInfo<3>::vertices_per_face);
3328 
3329  // set up a table to make sure that
3330  // we handle non-standard faces correctly
3331  //
3332  // so set up a table that for each vertex (of
3333  // a quad in standard position) describes
3334  // which vertex to take
3335  //
3336  // first index: four vertices 0...3
3337  //
3338  // second index: face_orientation; 0:
3339  // opposite normal, 1: standard
3340  //
3341  // third index: face_flip; 0: standard, 1:
3342  // face rotated by 180 degrees
3343  //
3344  // forth index: face_rotation: 0: standard,
3345  // 1: face rotated by 90 degrees
3346 
3347  constexpr unsigned int vertex_translation[4][2][2][2] = {
3348  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3349  // face_rotation=false and true
3350  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3351  // face_rotation=false and true
3352  {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3353  // face_rotation=false and true
3354  {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3355  // face_rotation=false and true
3356 
3357  {{{2, 3}, // vertex 1 ...
3358  {1, 0}},
3359  {{1, 0}, {2, 3}}},
3360 
3361  {{{1, 0}, // vertex 2 ...
3362  {2, 3}},
3363  {{2, 3}, {1, 0}}},
3364 
3365  {{{3, 1}, // vertex 3 ...
3366  {0, 2}},
3367  {{3, 1}, {0, 2}}}};
3368 
3369  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3370 }
3371 
3372 
3373 
3374 template <int dim>
3375 inline unsigned int
3376 GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3377  const bool,
3378  const bool,
3379  const bool)
3380 {
3381  Assert(dim > 1, ExcImpossibleInDim(dim));
3382  AssertIndexRange(vertex, GeometryInfo<dim>::vertices_per_face);
3383  return vertex;
3384 }
3385 
3386 template <int dim>
3387 inline unsigned int
3388 GeometryInfo<dim>::n_children(const RefinementCase<dim> &ref_case)
3389 {
3390  constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3391  0, 2, 2, 4, 2, 4, 4, 8};
3392 
3393  return n_children[ref_case];
3394 }
3395 
3396 
3397 
3398 template <int dim>
3399 inline unsigned int
3400 GeometryInfo<dim>::n_subfaces(const internal::SubfaceCase<dim> &)
3401 {
3402  Assert(false, ExcNotImplemented());
3403  return 0;
3404 }
3405 
3406 template <>
3407 inline unsigned int
3408 GeometryInfo<1>::n_subfaces(const internal::SubfaceCase<1> &)
3409 {
3410  Assert(false, ExcImpossibleInDim(1));
3411  return 0;
3412 }
3413 
3414 template <>
3415 inline unsigned int
3416 GeometryInfo<2>::n_subfaces(const internal::SubfaceCase<2> &subface_case)
3417 {
3418  return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3419 }
3420 
3421 
3422 
3423 template <>
3424 inline unsigned int
3425 GeometryInfo<3>::n_subfaces(const internal::SubfaceCase<3> &subface_case)
3426 {
3427  const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3428  0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3429  return nsubs[subface_case];
3430 }
3431 
3432 
3433 
3434 template <int dim>
3435 inline double
3436 GeometryInfo<dim>::subface_ratio(const internal::SubfaceCase<dim> &,
3437  const unsigned int)
3438 {
3439  Assert(false, ExcNotImplemented());
3440  return 0.;
3441 }
3442 
3443 template <>
3444 inline double
3445 GeometryInfo<1>::subface_ratio(const internal::SubfaceCase<1> &,
3446  const unsigned int)
3447 {
3448  return 1;
3449 }
3450 
3451 
3452 template <>
3453 inline double
3454 GeometryInfo<2>::subface_ratio(const internal::SubfaceCase<2> &subface_case,
3455  const unsigned int)
3456 {
3457  double ratio = 1;
3458  switch (subface_case)
3459  {
3460  case internal::SubfaceCase<2>::case_none:
3461  // Here, an
3462  // Assert(false,ExcInternalError())
3463  // would be the right
3464  // choice, but
3465  // unfortunately the
3466  // current function is
3467  // also called for faces
3468  // without children (see
3469  // tests/fe/mapping.cc).
3470  // Assert(false, ExcMessage("Face has no subfaces."));
3471  // Furthermore, assign
3472  // following value as
3473  // otherwise the
3474  // bits/volume_x tests
3475  // break
3476  ratio = 1. / GeometryInfo<2>::max_children_per_face;
3477  break;
3478  case internal::SubfaceCase<2>::case_x:
3479  ratio = 0.5;
3480  break;
3481  default:
3482  // there should be no
3483  // cases left
3484  Assert(false, ExcInternalError());
3485  break;
3486  }
3487 
3488  return ratio;
3489 }
3490 
3491 
3492 template <>
3493 inline double
3494 GeometryInfo<3>::subface_ratio(const internal::SubfaceCase<3> &subface_case,
3495  const unsigned int subface_no)
3496 {
3497  double ratio = 1;
3498  switch (subface_case)
3499  {
3500  case internal::SubfaceCase<3>::case_none:
3501  // Here, an
3502  // Assert(false,ExcInternalError())
3503  // would be the right
3504  // choice, but
3505  // unfortunately the
3506  // current function is
3507  // also called for faces
3508  // without children (see
3509  // tests/bits/mesh_3d_16.cc). Add
3510  // following switch to
3511  // avoid diffs in
3512  // tests/bits/mesh_3d_16
3513  ratio = 1. / GeometryInfo<3>::max_children_per_face;
3514  break;
3515  case internal::SubfaceCase<3>::case_x:
3516  case internal::SubfaceCase<3>::case_y:
3517  ratio = 0.5;
3518  break;
3519  case internal::SubfaceCase<3>::case_xy:
3520  case internal::SubfaceCase<3>::case_x1y2y:
3521  case internal::SubfaceCase<3>::case_y1x2x:
3522  ratio = 0.25;
3523  break;
3524  case internal::SubfaceCase<3>::case_x1y:
3525  case internal::SubfaceCase<3>::case_y1x:
3526  if (subface_no < 2)
3527  ratio = 0.25;
3528  else
3529  ratio = 0.5;
3530  break;
3531  case internal::SubfaceCase<3>::case_x2y:
3532  case internal::SubfaceCase<3>::case_y2x:
3533  if (subface_no == 0)
3534  ratio = 0.5;
3535  else
3536  ratio = 0.25;
3537  break;
3538  default:
3539  // there should be no
3540  // cases left
3541  Assert(false, ExcInternalError());
3542  break;
3543  }
3544 
3545  return ratio;
3546 }
3547 
3548 
3549 
3550 template <int dim>
3551 RefinementCase<dim - 1> inline GeometryInfo<dim>::face_refinement_case(
3552  const RefinementCase<dim> &,
3553  const unsigned int,
3554  const bool,
3555  const bool,
3556  const bool)
3557 {
3558  Assert(false, ExcNotImplemented());
3559  return RefinementCase<dim - 1>::no_refinement;
3560 }
3561 
3562 template <>
3563 RefinementCase<0> inline GeometryInfo<1>::face_refinement_case(
3564  const RefinementCase<1> &,
3565  const unsigned int,
3566  const bool,
3567  const bool,
3568  const bool)
3569 {
3570  Assert(false, ExcImpossibleInDim(1));
3571 
3572  return RefinementCase<0>::no_refinement;
3573 }
3574 
3575 
3576 template <>
3577 inline RefinementCase<1>
3578 GeometryInfo<2>::face_refinement_case(
3579  const RefinementCase<2> &cell_refinement_case,
3580  const unsigned int face_no,
3581  const bool,
3582  const bool,
3583  const bool)
3584 {
3585  const unsigned int dim = 2;
3586  AssertIndexRange(cell_refinement_case,
3587  RefinementCase<dim>::isotropic_refinement + 1);
3588  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
3589 
3590  const RefinementCase<dim - 1>
3591  ref_cases[RefinementCase<dim>::isotropic_refinement +
3592  1][GeometryInfo<dim>::faces_per_cell / 2] = {
3593  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3594  RefinementCase<dim - 1>::no_refinement},
3595 
3596  {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3597 
3598  {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3599 
3600  {RefinementCase<dim - 1>::cut_x, // cut_xy
3601  RefinementCase<dim - 1>::cut_x}};
3602 
3603  return ref_cases[cell_refinement_case][face_no / 2];
3604 }
3605 
3606 
3607 template <>
3608 inline RefinementCase<2>
3609 GeometryInfo<3>::face_refinement_case(
3610  const RefinementCase<3> &cell_refinement_case,
3611  const unsigned int face_no,
3612  const bool face_orientation,
3613  const bool /*face_flip*/,
3614  const bool face_rotation)
3615 {
3616  const unsigned int dim = 3;
3617  AssertIndexRange(cell_refinement_case,
3618  RefinementCase<dim>::isotropic_refinement + 1);
3619  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
3620 
3621  const RefinementCase<dim - 1>
3622  ref_cases[RefinementCase<dim>::isotropic_refinement + 1]
3623  [GeometryInfo<dim>::faces_per_cell / 2] = {
3624  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3625  RefinementCase<dim - 1>::no_refinement,
3626  RefinementCase<dim - 1>::no_refinement},
3627 
3628  {RefinementCase<dim - 1>::no_refinement, // cut_x
3629  RefinementCase<dim - 1>::cut_y,
3630  RefinementCase<dim - 1>::cut_x},
3631 
3632  {RefinementCase<dim - 1>::cut_x, // cut_y
3633  RefinementCase<dim - 1>::no_refinement,
3634  RefinementCase<dim - 1>::cut_y},
3635 
3636  {RefinementCase<dim - 1>::cut_x, // cut_xy
3637  RefinementCase<dim - 1>::cut_y,
3638  RefinementCase<dim - 1>::cut_xy},
3639 
3640  {RefinementCase<dim - 1>::cut_y, // cut_z
3641  RefinementCase<dim - 1>::cut_x,
3642  RefinementCase<dim - 1>::no_refinement},
3643 
3644  {RefinementCase<dim - 1>::cut_y, // cut_xz
3645  RefinementCase<dim - 1>::cut_xy,
3646  RefinementCase<dim - 1>::cut_x},
3647 
3648  {RefinementCase<dim - 1>::cut_xy, // cut_yz
3649  RefinementCase<dim - 1>::cut_x,
3650  RefinementCase<dim - 1>::cut_y},
3651 
3652  {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3653  RefinementCase<dim - 1>::cut_xy,
3654  RefinementCase<dim - 1>::cut_xy},
3655  };
3656 
3657  const RefinementCase<dim - 1> ref_case =
3658  ref_cases[cell_refinement_case][face_no / 2];
3659 
3660  const RefinementCase<dim - 1> flip[4] = {
3661  RefinementCase<dim - 1>::no_refinement,
3662  RefinementCase<dim - 1>::cut_y,
3663  RefinementCase<dim - 1>::cut_x,
3664  RefinementCase<dim - 1>::cut_xy};
3665 
3666  // correct the ref_case for face_orientation
3667  // and face_rotation. for face_orientation,
3668  // 'true' is the default value whereas for
3669  // face_rotation, 'false' is standard. If
3670  // <tt>face_rotation==face_orientation</tt>,
3671  // then one of them is non-standard and we
3672  // have to swap cut_x and cut_y, otherwise no
3673  // change is necessary. face_flip has no
3674  // influence. however, in order to keep the
3675  // interface consistent with other functions,
3676  // we still include it as an argument to this
3677  // function
3678  return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3679 }
3680 
3681 
3682 
3683 template <int dim>
3684 inline RefinementCase<1>
3685 GeometryInfo<dim>::line_refinement_case(const RefinementCase<dim> &,
3686  const unsigned int)
3687 {
3688  Assert(false, ExcNotImplemented());
3689  return RefinementCase<1>::no_refinement;
3690 }
3691 
3692 template <>
3693 inline RefinementCase<1>
3694 GeometryInfo<1>::line_refinement_case(
3695  const RefinementCase<1> &cell_refinement_case,
3696  const unsigned int line_no)
3697 {
3698  (void)line_no;
3699  const unsigned int dim = 1;
3700  (void)dim;
3701  AssertIndexRange(cell_refinement_case,
3702  RefinementCase<dim>::isotropic_refinement + 1);
3703  AssertIndexRange(line_no, GeometryInfo<dim>::lines_per_cell);
3704 
3705  return cell_refinement_case;
3706 }
3707 
3708 
3709 template <>
3710 inline RefinementCase<1>
3711 GeometryInfo<2>::line_refinement_case(
3712  const RefinementCase<2> &cell_refinement_case,
3713  const unsigned int line_no)
3714 {
3715  // Assertions are in face_refinement_case()
3716  return face_refinement_case(cell_refinement_case, line_no);
3717 }
3718 
3719 
3720 template <>
3721 inline RefinementCase<1>
3722 GeometryInfo<3>::line_refinement_case(
3723  const RefinementCase<3> &cell_refinement_case,
3724  const unsigned int line_no)
3725 {
3726  const unsigned int dim = 3;
3727  AssertIndexRange(cell_refinement_case,
3728  RefinementCase<dim>::isotropic_refinement + 1);
3729  AssertIndexRange(line_no, GeometryInfo<dim>::lines_per_cell);
3730 
3731  // array indicating, which simple refine
3732  // case cuts a line in direction x, y or
3733  // z. For example, cut_y and everything
3734  // containing cut_y (cut_xy, cut_yz,
3735  // cut_xyz) cuts lines, which are in y
3736  // direction.
3737  const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3738  RefinementCase<dim>::cut_y,
3739  RefinementCase<dim>::cut_z};
3740 
3741  // order the direction of lines
3742  // 0->x, 1->y, 2->z
3743  const unsigned int direction[lines_per_cell] = {
3744  1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3745 
3746  return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3747  RefinementCase<1>::cut_x :
3748  RefinementCase<1>::no_refinement);
3749 }
3750 
3751 
3752 
3753 template <int dim>
3754 inline RefinementCase<dim>
3755 GeometryInfo<dim>::min_cell_refinement_case_for_face_refinement(
3756  const RefinementCase<dim - 1> &,
3757  const unsigned int,
3758  const bool,
3759  const bool,
3760  const bool)
3761 {
3762  Assert(false, ExcNotImplemented());
3763 
3764  return RefinementCase<dim>::no_refinement;
3765 }
3766 
3767 template <>
3768 inline RefinementCase<1>
3769 GeometryInfo<1>::min_cell_refinement_case_for_face_refinement(
3770  const RefinementCase<0> &,
3771  const unsigned int,
3772  const bool,
3773  const bool,
3774  const bool)
3775 {
3776  const unsigned int dim = 1;
3777  Assert(false, ExcImpossibleInDim(dim));
3778 
3779  return RefinementCase<dim>::no_refinement;
3780 }
3781 
3782 
3783 template <>
3784 inline RefinementCase<2>
3785 GeometryInfo<2>::min_cell_refinement_case_for_face_refinement(
3786  const RefinementCase<1> &face_refinement_case,
3787  const unsigned int face_no,
3788  const bool,
3789  const bool,
3790  const bool)
3791 {
3792  const unsigned int dim = 2;
3793  AssertIndexRange(face_refinement_case,
3794  RefinementCase<dim - 1>::isotropic_refinement + 1);
3795  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
3796 
3797  if (face_refinement_case == RefinementCase<dim>::cut_x)
3798  return (face_no / 2) ? RefinementCase<dim>::cut_x :
3799  RefinementCase<dim>::cut_y;
3800  else
3801  return RefinementCase<dim>::no_refinement;
3802 }
3803 
3804 
3805 template <>
3806 inline RefinementCase<3>
3807 GeometryInfo<3>::min_cell_refinement_case_for_face_refinement(
3808  const RefinementCase<2> &face_refinement_case,
3809  const unsigned int face_no,
3810  const bool face_orientation,
3811  const bool /*face_flip*/,
3812  const bool face_rotation)
3813 {
3814  const unsigned int dim = 3;
3815  AssertIndexRange(face_refinement_case,
3816  RefinementCase<dim - 1>::isotropic_refinement + 1);
3817  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
3818 
3819  const RefinementCase<2> flip[4] = {RefinementCase<2>::no_refinement,
3820  RefinementCase<2>::cut_y,
3821  RefinementCase<2>::cut_x,
3822  RefinementCase<2>::cut_xy};
3823 
3824  // correct the face_refinement_case for
3825  // face_orientation and face_rotation. for
3826  // face_orientation, 'true' is the default
3827  // value whereas for face_rotation, 'false'
3828  // is standard. If
3829  // <tt>face_rotation==face_orientation</tt>,
3830  // then one of them is non-standard and we
3831  // have to swap cut_x and cut_y, otherwise no
3832  // change is necessary. face_flip has no
3833  // influence. however, in order to keep the
3834  // interface consistent with other functions,
3835  // we still include it as an argument to this
3836  // function
3837  const RefinementCase<dim - 1> std_face_ref =
3838  (face_orientation == face_rotation) ? flip[face_refinement_case] :
3839  face_refinement_case;
3840 
3841  const RefinementCase<dim> face_to_cell[3][4] = {
3842  {RefinementCase<dim>::no_refinement, // faces 0 and 1
3843  RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3844  RefinementCase<dim>::cut_z,
3845  RefinementCase<dim>::cut_yz},
3846 
3847  {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3848  // "exchanged on faces 2 and 3")
3849  RefinementCase<dim>::cut_z,
3850  RefinementCase<dim>::cut_x,
3851  RefinementCase<dim>::cut_xz},
3852 
3853  {RefinementCase<dim>::no_refinement, // faces 4 and 5
3854  RefinementCase<dim>::cut_x,
3855  RefinementCase<dim>::cut_y,
3856  RefinementCase<dim>::cut_xy}};
3857 
3858  return face_to_cell[face_no / 2][std_face_ref];
3859 }
3860 
3861 
3862 
3863 template <int dim>
3864 inline RefinementCase<dim>
3865 GeometryInfo<dim>::min_cell_refinement_case_for_line_refinement(
3866  const unsigned int)
3867 {
3868  Assert(false, ExcNotImplemented());
3869 
3870  return RefinementCase<dim>::no_refinement;
3871 }
3872 
3873 template <>
3874 inline RefinementCase<1>
3875 GeometryInfo<1>::min_cell_refinement_case_for_line_refinement(
3876  const unsigned int line_no)
3877 {
3878  (void)line_no;
3879  AssertIndexRange(line_no, 1);
3880 
3881  return RefinementCase<1>::cut_x;
3882 }
3883 
3884 
3885 template <>
3886 inline RefinementCase<2>
3887 GeometryInfo<2>::min_cell_refinement_case_for_line_refinement(
3888  const unsigned int line_no)
3889 {
3890  const unsigned int dim = 2;
3891  (void)dim;
3892  AssertIndexRange(line_no, GeometryInfo<dim>::lines_per_cell);
3893 
3894  return (line_no / 2) ? RefinementCase<2>::cut_x : RefinementCase<2>::cut_y;
3895 }
3896 
3897 
3898 template <>
3899 inline RefinementCase<3>
3900 GeometryInfo<3>::min_cell_refinement_case_for_line_refinement(
3901  const unsigned int line_no)
3902 {
3903  const unsigned int dim = 3;
3904  AssertIndexRange(line_no, GeometryInfo<dim>::lines_per_cell);
3905 
3906  const RefinementCase<dim> ref_cases[6] = {
3907  RefinementCase<dim>::cut_y, // lines 0 and 1
3908  RefinementCase<dim>::cut_x, // lines 2 and 3
3909  RefinementCase<dim>::cut_y, // lines 4 and 5
3910  RefinementCase<dim>::cut_x, // lines 6 and 7
3911  RefinementCase<dim>::cut_z, // lines 8 and 9
3912  RefinementCase<dim>::cut_z}; // lines 10 and 11
3913 
3914  return ref_cases[line_no / 2];
3915 }
3916 
3917 
3918 
3919 template <>
3920 inline unsigned int
3921 GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
3922  const bool face_orientation,
3923  const bool face_flip,
3924  const bool face_rotation)
3925 {
3926  AssertIndexRange(vertex, GeometryInfo<3>::vertices_per_face);
3927 
3928  // set up a table to make sure that
3929  // we handle non-standard faces correctly
3930  //
3931  // so set up a table that for each vertex (of
3932  // a quad in standard position) describes
3933  // which vertex to take
3934  //
3935  // first index: four vertices 0...3
3936  //
3937  // second index: face_orientation; 0:
3938  // opposite normal, 1: standard
3939  //
3940  // third index: face_flip; 0: standard, 1:
3941  // face rotated by 180 degrees
3942  //
3943  // forth index: face_rotation: 0: standard,
3944  // 1: face rotated by 90 degrees
3945 
3946  const unsigned int vertex_translation[4][2][2][2] = {
3947  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3948  // face_rotation=false and true
3949  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3950  // face_rotation=false and true
3951  {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
3952  // face_rotation=false and true
3953  {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
3954  // face_rotation=false and true
3955 
3956  {{{2, 3}, // vertex 1 ...
3957  {1, 0}},
3958  {{1, 3}, {2, 0}}},
3959 
3960  {{{1, 0}, // vertex 2 ...
3961  {2, 3}},
3962  {{2, 0}, {1, 3}}},
3963 
3964  {{{3, 1}, // vertex 3 ...
3965  {0, 2}},
3966  {{3, 2}, {0, 1}}}};
3967 
3968  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3969 }
3970 
3971 
3972 
3973 template <int dim>
3974 inline unsigned int
3975 GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
3976  const bool,
3977  const bool,
3978  const bool)
3979 {
3980  Assert(dim > 1, ExcImpossibleInDim(dim));
3981  AssertIndexRange(vertex, GeometryInfo<dim>::vertices_per_face);
3982  return vertex;
3983 }
3984 
3985 
3986 
3987 template <>
3988 inline unsigned int
3989 GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
3990  const bool face_orientation,
3991  const bool face_flip,
3992  const bool face_rotation)
3993 {
3994  AssertIndexRange(line, GeometryInfo<3>::lines_per_face);
3995 
3996 
3997  // make sure we handle
3998  // non-standard faces correctly
3999  //
4000  // so set up a table that for each line (of a
4001  // quad) describes which line to take
4002  //
4003  // first index: four lines 0...3
4004  //
4005  // second index: face_orientation; 0:
4006  // opposite normal, 1: standard
4007  //
4008  // third index: face_flip; 0: standard, 1:
4009  // face rotated by 180 degrees
4010  //
4011  // forth index: face_rotation: 0: standard,
4012  // 1: face rotated by 90 degrees
4013 
4014  const unsigned int line_translation[4][2][2][2] = {
4015  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4016  // face_rotation=false and true
4017  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4018  // face_rotation=false and true
4019  {{0, 3}, // line 0, face_orientation=true, face_flip=false,
4020  // face_rotation=false and true
4021  {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
4022  // face_rotation=false and true
4023 
4024  {{{3, 1}, // line 1 ...
4025  {2, 0}},
4026  {{1, 2}, {0, 3}}},
4027 
4028  {{{0, 3}, // line 2 ...
4029  {1, 2}},
4030  {{2, 0}, {3, 1}}},
4031 
4032  {{{1, 2}, // line 3 ...
4033  {0, 3}},
4034  {{3, 1}, {2, 0}}}};
4035 
4036  return line_translation[line][face_orientation][face_flip][face_rotation];
4037 }
4038 
4039 
4040 
4041 template <int dim>
4042 inline unsigned int
4043 GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
4044  const bool,
4045  const bool,
4046  const bool)
4047 {
4048  Assert(false, ExcNotImplemented());
4049  return line;
4050 }
4051 
4052 
4053 
4054 template <>
4055 inline unsigned int
4056 GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
4057  const bool face_orientation,
4058  const bool face_flip,
4059  const bool face_rotation)
4060 {
4061  AssertIndexRange(line, GeometryInfo<3>::lines_per_face);
4062 
4063 
4064  // make sure we handle
4065  // non-standard faces correctly
4066  //
4067  // so set up a table that for each line (of a
4068  // quad) describes which line to take
4069  //
4070  // first index: four lines 0...3
4071  //
4072  // second index: face_orientation; 0:
4073  // opposite normal, 1: standard
4074  //
4075  // third index: face_flip; 0: standard, 1:
4076  // face rotated by 180 degrees
4077  //
4078  // forth index: face_rotation: 0: standard,
4079  // 1: face rotated by 90 degrees
4080 
4081  const unsigned int line_translation[4][2][2][2] = {
4082  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4083  // face_rotation=false and true
4084  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4085  // face_rotation=false and true
4086  {{0, 2}, // line 0, face_orientation=true, face_flip=false,
4087  // face_rotation=false and true
4088  {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
4089  // face_rotation=false and true
4090 
4091  {{{3, 1}, // line 1 ...
4092  {2, 0}},
4093  {{1, 3}, {0, 2}}},
4094 
4095  {{{0, 3}, // line 2 ...
4096  {1, 2}},
4097  {{2, 1}, {3, 0}}},
4098 
4099  {{{1, 2}, // line 3 ...
4100  {0, 3}},
4101  {{3, 0}, {2, 1}}}};
4102 
4103  return line_translation[line][face_orientation][face_flip][face_rotation];
4104 }
4105 
4106 
4107 
4108 template <int dim>
4109 inline unsigned int
4110 GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
4111  const bool,
4112  const bool,
4113  const bool)
4114 {
4115  Assert(false, ExcNotImplemented());
4116  return line;
4117 }
4118 
4119 
4120 
4121 template <>
4122 inline unsigned int
4123 GeometryInfo<1>::child_cell_on_face(const RefinementCase<1> &,
4124  const unsigned int face,
4125  const unsigned int subface,
4126  const bool,
4127  const bool,
4128  const bool,
4129  const RefinementCase<0> &)
4130 {
4131  (void)subface;
4132  AssertIndexRange(face, faces_per_cell);
4133  AssertIndexRange(subface, max_children_per_face);
4134 
4135  return face;
4136 }
4137 
4138 
4139 
4140 template <>
4141 inline unsigned int
4142 GeometryInfo<2>::child_cell_on_face(const RefinementCase<2> &ref_case,
4143  const unsigned int face,
4144  const unsigned int subface,
4145  const bool /*face_orientation*/,
4146  const bool face_flip,
4147  const bool /*face_rotation*/,
4148  const RefinementCase<1> &)
4149 {
4150  AssertIndexRange(face, faces_per_cell);
4151  AssertIndexRange(subface, max_children_per_face);
4152 
4153  // always return the child adjacent to the specified
4154  // subface. if the face of a cell is not refined, don't
4155  // throw an assertion but deliver the child adjacent to
4156  // the face nevertheless, i.e. deliver the child of
4157  // this cell adjacent to the subface of a possibly
4158  // refined neighbor. this simplifies setting neighbor
4159  // information in execute_refinement.
4160  const unsigned int
4161  subcells[2][RefinementCase<2>::isotropic_refinement][faces_per_cell]
4162  [max_children_per_face] = {
4163  {
4164  // Normal orientation (face_flip = false)
4165  {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
4166  {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
4167  {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
4168  },
4169  {
4170  // Flipped orientation (face_flip = true)
4171  {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
4172  {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
4173  {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
4174  }};
4175 
4176  return subcells[face_flip][ref_case - 1][face][subface];
4177 }
4178 
4179 
4180 
4181 template <>
4182 inline unsigned int
4183 GeometryInfo<3>::child_cell_on_face(const RefinementCase<3> &ref_case,
4184  const unsigned int face,
4185  const unsigned int subface,
4186  const bool face_orientation,
4187  const bool face_flip,
4188  const bool face_rotation,
4189  const RefinementCase<2> &face_ref_case)
4190 {
4191  const unsigned int dim = 3;
4192 
4193  Assert(ref_case > RefinementCase<dim - 1>::no_refinement,
4194  ExcMessage("Cell has no children."));
4195  AssertIndexRange(face, faces_per_cell);
4196  if (!(subface == 0 &&
4197  face_ref_case == RefinementCase<dim - 1>::no_refinement))
4198  {
4199  AssertIndexRange(subface,
4200  GeometryInfo<dim - 1>::n_children(face_ref_case));
4201  }
4202 
4203  // invalid number used for invalid cases,
4204  // e.g. when the children are more refined at
4205  // a given face than the face itself
4206  const unsigned int e = numbers::invalid_unsigned_int;
4207 
4208  // the whole process of finding a child cell
4209  // at a given subface considering the
4210  // possibly anisotropic refinement cases of
4211  // the cell and the face as well as
4212  // orientation, flip and rotation of the face
4213  // is quite complicated. thus, we break it
4214  // down into several steps.
4215 
4216  // first step: convert the given face refine
4217  // case to a face refine case concerning the
4218  // face in standard orientation (, flip and
4219  // rotation). This only affects cut_x and
4220  // cut_y
4221  const RefinementCase<dim - 1> flip[4] = {
4222  RefinementCase<dim - 1>::no_refinement,
4223  RefinementCase<dim - 1>::cut_y,
4224  RefinementCase<dim - 1>::cut_x,
4225  RefinementCase<dim - 1>::cut_xy};
4226  // for face_orientation, 'true' is the
4227  // default value whereas for face_rotation,
4228  // 'false' is standard. If
4229  // <tt>face_rotation==face_orientation</tt>,
4230  // then one of them is non-standard and we
4231  // have to swap cut_x and cut_y, otherwise no
4232  // change is necessary.
4233  const RefinementCase<dim - 1> std_face_ref =
4234  (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4235 
4236  // second step: convert the given subface
4237  // index to the one for a standard face
4238  // respecting face_orientation, face_flip and
4239  // face_rotation
4240 
4241  // first index: face_ref_case
4242  // second index: face_orientation
4243  // third index: face_flip
4244  // forth index: face_rotation
4245  // fifth index: subface index
4246  const unsigned int subface_exchange[4][2][2][2][4] = {
4247  // no_refinement (subface 0 stays 0,
4248  // all others are invalid)
4249  {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4250  {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4251  // cut_x (here, if the face is only
4252  // rotated OR only falsely oriented,
4253  // then subface 0 of the non-standard
4254  // face does NOT correspond to one of
4255  // the subfaces of a standard
4256  // face. Thus we indicate the subface
4257  // which is located at the lower left
4258  // corner (the origin of the face's
4259  // local coordinate system) with
4260  // '0'. The rest of this issue is
4261  // taken care of using the above
4262  // conversion to a 'standard face
4263  // refine case')
4264  {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4265  {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4266  // cut_y (the same applies as for
4267  // cut_x)
4268  {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4269  {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4270  // cut_xyz: this information is
4271  // identical to the information
4272  // returned by
4273  // GeometryInfo<3>::real_to_standard_face_vertex()
4274  {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4275  // face_rotation=false, subfaces 0,1,2,3
4276  {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4277  // face_rotation=true, subfaces 0,1,2,3
4278  {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4279  // face_rotation=false, subfaces 0,1,2,3
4280  {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4281  // face_rotation=true, subfaces 0,1,2,3
4282  {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4283  // face_rotation=false, subfaces 0,1,2,3
4284  {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4285  // face_rotation=true, subfaces 0,1,2,3
4286  {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4287  // face_rotation=false, subfaces 0,1,2,3
4288  {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4289  // face_rotation=true, subfaces 0,1,2,3
4290 
4291  const unsigned int std_subface =
4292  subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4293  [subface];
4294  Assert(std_subface != e, ExcInternalError());
4295 
4296  // third step: these are the children, which
4297  // can be found at the given subfaces of an
4298  // isotropically refined (standard) face
4299  //
4300  // first index: (refinement_case-1)
4301  // second index: face_index
4302  // third index: subface_index (isotropic refinement)
4303  const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4304  [max_children_per_face] = {
4305  // cut_x
4306  {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4307  {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4308  {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4309  {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4310  {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4311  {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4312  // cut_y
4313  {{0, 1, 0, 1},
4314  {0, 1, 0, 1},
4315  {0, 0, 0, 0},
4316  {1, 1, 1, 1},
4317  {0, 0, 1, 1},
4318  {0, 0, 1, 1}},
4319  // cut_xy
4320  {{0, 2, 0, 2},
4321  {1, 3, 1, 3},
4322  {0, 0, 1, 1},
4323  {2, 2, 3, 3},
4324  {0, 1, 2, 3},
4325  {0, 1, 2, 3}},
4326  // cut_z
4327  {{0, 0, 1, 1},
4328  {0, 0, 1, 1},
4329  {0, 1, 0, 1},
4330  {0, 1, 0, 1},
4331  {0, 0, 0, 0},
4332  {1, 1, 1, 1}},
4333  // cut_xz
4334  {{0, 0, 1, 1},
4335  {2, 2, 3, 3},
4336  {0, 1, 2, 3},
4337  {0, 1, 2, 3},
4338  {0, 2, 0, 2},
4339  {1, 3, 1, 3}},
4340  // cut_yz
4341  {{0, 1, 2, 3},
4342  {0, 1, 2, 3},
4343  {0, 2, 0, 2},
4344  {1, 3, 1, 3},
4345  {0, 0, 1, 1},
4346  {2, 2, 3, 3}},
4347  // cut_xyz
4348  {{0, 2, 4, 6},
4349  {1, 3, 5, 7},
4350  {0, 4, 1, 5},
4351  {2, 6, 3, 7},
4352  {0, 1, 2, 3},
4353  {4, 5, 6, 7}}};
4354 
4355  // forth step: check, whether the given face
4356  // refine case is valid for the given cell
4357  // refine case. this is the case, if the
4358  // given face refine case is at least as
4359  // refined as the face is for the given cell
4360  // refine case
4361 
4362  // note, that we are considering standard
4363  // face refinement cases here and thus must
4364  // not pass the given orientation, flip and
4365  // rotation flags
4366  if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4367  face_refinement_case(ref_case, face))
4368  {
4369  // all is fine. for anisotropic face
4370  // refine cases, select one of the
4371  // isotropic subfaces which neighbors the
4372  // same child
4373 
4374  // first index: (standard) face refine case
4375  // second index: subface index
4376  const unsigned int equivalent_iso_subface[4][4] = {
4377  {0, e, e, e}, // no_refinement
4378  {0, 3, e, e}, // cut_x
4379  {0, 3, e, e}, // cut_y
4380  {0, 1, 2, 3}}; // cut_xy
4381 
4382  const unsigned int equ_std_subface =
4383  equivalent_iso_subface[std_face_ref][std_subface];
4384  Assert(equ_std_subface != e, ExcInternalError());
4385 
4386  return iso_children[ref_case - 1][face][equ_std_subface];
4387  }
4388  else
4389  {
4390  // the face_ref_case was too coarse,
4391  // throw an error
4392  Assert(false,
4393  ExcMessage("The face RefineCase is too coarse "
4394  "for the given cell RefineCase."));
4395  }
4396  // we only get here in case of an error
4397  return e;
4398 }
4399 
4400 
4401 
4402 template <>
4403 inline unsigned int
4404 GeometryInfo<4>::child_cell_on_face(const RefinementCase<4> &,
4405  const unsigned int,
4406  const unsigned int,
4407  const bool,
4408  const bool,
4409  const bool,
4410  const RefinementCase<3> &)
4411 {
4412  Assert(false, ExcNotImplemented());
4413  return numbers::invalid_unsigned_int;
4414 }
4415 
4416 
4417 
4418 template <>
4419 inline unsigned int
4420 GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4421  const unsigned int line,
4422  const bool,
4423  const bool,
4424  const bool)
4425 {
4426  (void)line;
4427  AssertIndexRange(face, faces_per_cell);
4428  AssertIndexRange(line, lines_per_face);
4429 
4430  // The face is a line itself.
4431  return face;
4432 }
4433 
4434 
4435 
4436 template <>
4437 inline unsigned int
4438 GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4439  const unsigned int line,
4440  const bool face_orientation,
4441  const bool face_flip,
4442  const bool face_rotation)
4443 {
4444  AssertIndexRange(face, faces_per_cell);
4445  AssertIndexRange(line, lines_per_face);
4446 
4447  const unsigned lines[faces_per_cell][lines_per_face] = {
4448  {8, 10, 0, 4}, // left face
4449  {9, 11, 1, 5}, // right face
4450  {2, 6, 8, 9}, // front face
4451  {3, 7, 10, 11}, // back face
4452  {0, 1, 2, 3}, // bottom face
4453  {4, 5, 6, 7}}; // top face
4454  return lines[face][real_to_standard_face_line(
4455  line, face_orientation, face_flip, face_rotation)];
4456 }
4457 
4458 
4459 
4460 template <int dim>
4461 inline unsigned int
4462 GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4463  const unsigned int,
4464  const bool,
4465  const bool,
4466  const bool)
4467 {
4468  Assert(false, ExcNotImplemented());
4469  return numbers::invalid_unsigned_int;
4470 }
4471 
4472 
4473 
4474 template <int dim>
4475 inline unsigned int
4476 GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4477  const unsigned int vertex,
4478  const bool face_orientation,
4479  const bool face_flip,
4480  const bool face_rotation)
4481 {
4482  return child_cell_on_face(RefinementCase<dim>::isotropic_refinement,
4483  face,
4484  vertex,
4485  face_orientation,
4486  face_flip,
4487  face_rotation);
4488 }
4489 
4490 
4491 
4492 template <int dim>
4493 inline Point<dim>
4494 GeometryInfo<dim>::project_to_unit_cell(const Point<dim> &q)
4495 {
4496  Point<dim> p = q;
4497  for (unsigned int i = 0; i < dim; i++)
4498  if (p[i] < 0.)
4499  p[i] = 0.;
4500  else if (p[i] > 1.)
4501  p[i] = 1.;
4502 
4503  return p;
4504 }
4505 
4506 
4507 
4508 template <int dim>
4509 inline double
4510 GeometryInfo<dim>::distance_to_unit_cell(const Point<dim> &p)
4511 {
4512  double result = 0.0;
4513 
4514  for (unsigned int i = 0; i < dim; i++)
4515  if ((-p[i]) > result)
4516  result = -p[i];
4517  else if ((p[i] - 1.) > result)
4518  result = (p[i] - 1.);
4519 
4520  return result;
4521 }
4522 
4523 
4524 
4525 template <int dim>
4526 inline double
4527 GeometryInfo<dim>::d_linear_shape_function(const Point<dim> & xi,
4528  const unsigned int i)
4529 {
4530  AssertIndexRange(i, GeometryInfo<dim>::vertices_per_cell);
4531 
4532  switch (dim)
4533  {
4534  case 1:
4535  {
4536  const double x = xi[0];
4537  switch (i)
4538  {
4539  case 0:
4540  return 1 - x;
4541  case 1:
4542  return x;
4543  }
4544  break;
4545  }
4546 
4547  case 2:
4548  {
4549  const double x = xi[0];
4550  const double y = xi[1];
4551  switch (i)
4552  {
4553  case 0:
4554  return (1 - x) * (1 - y);
4555  case 1:
4556  return x * (1 - y);
4557  case 2:
4558  return (1 - x) * y;
4559  case 3:
4560  return x * y;
4561  }
4562  break;
4563  }
4564 
4565  case 3:
4566  {
4567  const double x = xi[0];
4568  const double y = xi[1];
4569  const double z = xi[2];
4570  switch (i)
4571  {
4572  case 0:
4573  return (1 - x) * (1 - y) * (1 - z);
4574  case 1:
4575  return x * (1 - y) * (1 - z);
4576  case 2:
4577  return (1 - x) * y * (1 - z);
4578  case 3:
4579  return x * y * (1 - z);
4580  case 4:
4581  return (1 - x) * (1 - y) * z;
4582  case 5:
4583  return x * (1 - y) * z;
4584  case 6:
4585  return (1 - x) * y * z;
4586  case 7:
4587  return x * y * z;
4588  }
4589  break;
4590  }
4591 
4592  default:
4593  Assert(false, ExcNotImplemented());
4594  }
4595  return -1e9;
4596 }
4597 
4598 
4599 
4600 template <>
4601 Tensor<1, 1> inline GeometryInfo<1>::d_linear_shape_function_gradient(
4602  const Point<1> &,
4603  const unsigned int i)
4604 {
4605  AssertIndexRange(i, GeometryInfo<1>::vertices_per_cell);
4606 
4607  switch (i)
4608  {
4609  case 0:
4610  return Point<1>(-1.);
4611  case 1:
4612  return Point<1>(1.);
4613  }
4614 
4615  return Point<1>(-1e9);
4616 }
4617 
4618 
4619 
4620 template <>
4621 Tensor<1, 2> inline GeometryInfo<2>::d_linear_shape_function_gradient(
4622  const Point<2> & xi,
4623  const unsigned int i)
4624 {
4625  AssertIndexRange(i, GeometryInfo<2>::vertices_per_cell);
4626 
4627  const double x = xi[0];
4628  const double y = xi[1];
4629  switch (i)
4630  {
4631  case 0:
4632  return Point<2>(-(1 - y), -(1 - x));
4633  case 1:
4634  return Point<2>(1 - y, -x);
4635  case 2:
4636  return Point<2>(-y, 1 - x);
4637  case 3:
4638  return Point<2>(y, x);
4639  }
4640  return Point<2>(-1e9, -1e9);
4641 }
4642 
4643 
4644 
4645 template <>
4646 Tensor<1, 3> inline GeometryInfo<3>::d_linear_shape_function_gradient(
4647  const Point<3> & xi,
4648  const unsigned int i)
4649 {
4650  AssertIndexRange(i, GeometryInfo<3>::vertices_per_cell);
4651 
4652  const double x = xi[0];
4653  const double y = xi[1];
4654  const double z = xi[2];
4655  switch (i)
4656  {
4657  case 0:
4658  return Point<3>(-(1 - y) * (1 - z),
4659  -(1 - x) * (1 - z),
4660  -(1 - x) * (1 - y));
4661  case 1:
4662  return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4663  case 2:
4664  return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4665  case 3:
4666  return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4667  case 4:
4668  return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4669  case 5:
4670  return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4671  case 6:
4672  return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4673  case 7:
4674  return Point<3>(y * z, x * z, x * y);
4675  }
4676 
4677  return Point<3>(-1e9, -1e9, -1e9);
4678 }
4679 
4680 
4681 
4682 template <int dim>
4683 inline Tensor<1, dim>
4684 GeometryInfo<dim>::d_linear_shape_function_gradient(const Point<dim> &,
4685  const unsigned int)
4686 {
4687  Assert(false, ExcNotImplemented());
4688  return Tensor<1, dim>();
4689 }
4690 
4691 
4692 unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4693 {
4694  return 0;
4695 }
4696 
4697 
4698 namespace internal
4699 {
4700  namespace GeometryInfoHelper
4701  {
4702  // wedge product of a single
4703  // vector in 2d: we just have to
4704  // rotate it by 90 degrees to the
4705  // right
4706  inline Tensor<1, 2>
4707  wedge_product(const Tensor<1, 2> (&derivative)[1])
4708  {
4709  Tensor<1, 2> result;
4710  result[0] = derivative[0][1];
4711  result[1] = -derivative[0][0];
4712 
4713  return result;
4714  }
4715 
4716 
4717  // wedge product of 2 vectors in
4718  // 3d is the cross product
4719  inline Tensor<1, 3>
4720  wedge_product(const Tensor<1, 3> (&derivative)[2])
4721  {
4722  return cross_product_3d(derivative[0], derivative[1]);
4723  }
4724 
4725 
4726  // wedge product of dim vectors
4727  // in dim-d: that's the
4728  // determinant of the matrix
4729  template <int dim>
4730  inline Tensor<0, dim>
4731  wedge_product(const Tensor<1, dim> (&derivative)[dim])
4732  {
4733  Tensor<2, dim> jacobian;
4734  for (unsigned int i = 0; i < dim; ++i)
4735  jacobian[i] = derivative[i];
4736 
4737  return determinant(jacobian);
4738  }
4739  } // namespace GeometryInfoHelper
4740 } // namespace internal
4741 
4742 
4743 template <int dim>
4744 template <int spacedim>
4745 inline void
4746 GeometryInfo<dim>::alternating_form_at_vertices
4747 # ifndef DEAL_II_CONSTEXPR_BUG
4748  (const Point<spacedim> (&vertices)[vertices_per_cell],
4749  Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell])
4750 # else
4751  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms)
4752 # endif
4753 {
4754  // for each of the vertices,
4755  // compute the alternating form
4756  // of the mapped unit
4757  // vectors. consider for
4758  // example the case of a quad
4759  // in spacedim==3: to do so, we
4760  // need to see how the
4761  // infinitesimal vectors
4762  // (d\xi_1,0) and (0,d\xi_2)
4763  // are transformed into
4764  // spacedim-dimensional space
4765  // and then form their cross
4766  // product (i.e. the wedge product
4767  // of two vectors). to this end, note
4768  // that
4769  // \vec x = sum_i \vec v_i phi_i(\vec xi)
4770  // so the transformed vectors are
4771  // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
4772  // and
4773  // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
4774  // which boils down to the columns
4775  // of the 3x2 matrix \grad_\xi x(\xi)
4776  //
4777  // a similar reasoning would
4778  // hold for all dim,spacedim
4779  // pairs -- we only have to
4780  // compute the wedge product of
4781  // the columns of the
4782  // derivatives
4783  for (unsigned int i = 0; i < vertices_per_cell; ++i)
4784  {
4785  Tensor<1, spacedim> derivatives[dim];
4786 
4787  for (unsigned int j = 0; j < vertices_per_cell; ++j)
4788  {
4789  const Tensor<1, dim> grad_phi_j =
4790  d_linear_shape_function_gradient(unit_cell_vertex(i), j);
4791  for (unsigned int l = 0; l < dim; ++l)
4792  derivatives[l] += vertices[j] * grad_phi_j[l];
4793  }
4794 
4795  forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
4796  }
4797 }
4798 
4799 #endif // DOXYGEN
4800 DEAL_II_NAMESPACE_CLOSE
4801 
4802 #endif
GeometryInfo::ExcInvalidCoordinate
static ::ExceptionBase & ExcInvalidCoordinate(double arg1)
iota_view.h
GeometryInfo< 0 >::ucd_to_deal
static const std::array< unsigned int, vertices_per_cell > ucd_to_deal
Definition: geometry_info.h:1374
GeometryInfo::min_cell_refinement_case_for_line_refinement
static RefinementCase< dim > min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
GeometryInfo::dx_to_deal
static constexpr std::array< unsigned int, vertices_per_cell > dx_to_deal
Definition: geometry_info.h:2069
GeometryInfo::n_children
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
RefinementPossibilities::no_refinement
@ no_refinement
Definition: geometry_info.h:530
GeometryInfo::max_children_per_cell
static constexpr unsigned int max_children_per_cell
Definition: geometry_info.h:1928
RefinementCase::value
std::uint8_t value
Definition: geometry_info.h:891
GeometryInfo::face_indices
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
GeometryInfo::unit_normal_direction
static constexpr std::array< unsigned int, faces_per_cell > unit_normal_direction
Definition: geometry_info.h:2495
RefinementCase::RefinementCase
RefinementCase()
GeometryInfo::unit_cell_vertex
static Point< dim > unit_cell_vertex(const unsigned int vertex)
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
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])
determinant
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:2645
GeometryInfo::vertices_per_cell
static constexpr unsigned int vertices_per_cell
Definition: geometry_info.h:1967
GeometryInfo::d_linear_shape_function
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
LAPACKSupport::U
static const char U
Definition: lapack_support.h:167
GeometryInfo::project_to_unit_cell
static Point< dim > project_to_unit_cell(const Point< dim > &p)
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
Definition: geometry_info.h:1224
GeometryInfo::lines_per_face
static constexpr unsigned int lines_per_face
Definition: geometry_info.h:1996
GeometryInfo::quads_per_face
static constexpr unsigned int quads_per_face
Definition: geometry_info.h:2002
GeometryInfo< 0 >::dx_to_deal
static const std::array< unsigned int, vertices_per_cell > dx_to_deal
Definition: geometry_info.h:1389
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
internal::SubfacePossibilities::case_isotropic
@ case_isotropic
Definition: geometry_info.h:933
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
GeometryInfo::ExcInvalidSubface
static ::ExceptionBase & ExcInvalidSubface(int arg1, int arg2, int arg3)
internal::SubfacePossibilities
Definition: geometry_info.h:918
internal::SubfaceCase::ExcInvalidSubfaceCase
static ::ExceptionBase & ExcInvalidSubfaceCase(int arg1)
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
RefinementPossibilities< 3 >::Possibilities
Possibilities
Definition: geometry_info.h:739
GeometryInfo::is_inside_unit_cell
static bool is_inside_unit_cell(const Point< dim > &p)
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)
GeometryInfo::vertices_per_face
static constexpr unsigned int vertices_per_face
Definition: geometry_info.h:1990
GeometryPrimitive::vertex
@ vertex
Definition: geometry_info.h:432
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)
internal::SubfacePossibilities< 2 >::Possibilities
Possibilities
Definition: geometry_info.h:1025
RefinementCase::operator~
RefinementCase operator~() const
internal::SubfaceCase::memory_consumption
static constexpr std::size_t memory_consumption()
GeometryPrimitive::Object
Object
Definition: geometry_info.h:427
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::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::child_cell_from_point
static unsigned int child_cell_from_point(const Point< dim > &p)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
internal::SubfaceCase
Definition: geometry_info.h:1172
RefinementCase::operator&
RefinementCase operator&(const RefinementCase &r) const
GeometryInfo::vertex_indices
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
Tensor
Definition: tensor.h:450
internal::SubfacePossibilities::case_none
@ case_none
Definition: geometry_info.h:928
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
RefinementCase
Definition: geometry_info.h:795
GeometryInfo::n_subfaces
static unsigned int n_subfaces(const internal::SubfaceCase< dim > &subface_case)
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)
RefinementPossibilities::isotropic_refinement
@ isotropic_refinement
Definition: geometry_info.h:543
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
GeometryPrimitive
Definition: geometry_info.h:418
RefinementPossibilities< 2 >::Possibilities
Possibilities
Definition: geometry_info.h:663
DeclException3
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:564
RefinementPossibilities
Definition: geometry_info.h:489
DeclException1
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:518
internal::SubfacePossibilities< 1 >::Possibilities
Possibilities
Definition: geometry_info.h:989
GeometryInfo::line_to_cell_vertices
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
GeometryInfo::unit_tangential_vectors
static constexpr std::array< std::array< Tensor< 1, dim >, dim - 1 >, faces_per_cell > unit_tangential_vectors
Definition: geometry_info.h:2546
RefinementCase::serialize
void serialize(Archive &ar, const unsigned int version)
GeometryPrimitive::GeometryPrimitive
GeometryPrimitive(const Object object)
exceptions.h
RefinementCase::operator|
RefinementCase operator|(const RefinementCase &r) const
internal::SubfaceCase::SubfaceCase
SubfaceCase(const typename SubfacePossibilities< dim >::Possibilities subface_possibility)
GeometryInfo::line_refinement_case
static RefinementCase< 1 > line_refinement_case(const RefinementCase< dim > &cell_refinement_case, const unsigned int line_no)
value
static const bool value
Definition: dof_tools_constraints.cc:433
vertices
Point< 3 > vertices[4]
Definition: data_out_base.cc:174
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
GeometryInfo::max_children_per_face
static constexpr unsigned int max_children_per_face
Definition: geometry_info.h:1961
GeometryInfo::unit_normal_vector
static constexpr std::array< Tensor< 1, dim >, faces_per_cell > unit_normal_vector
Definition: geometry_info.h:2528
cross_product_3d
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2407
internal::SubfacePossibilities::Possibilities
Possibilities
Definition: geometry_info.h:923
DataOutBase::eps
@ eps
Definition: data_out_base.h:1582
GeometryInfo::lines_per_cell
static constexpr unsigned int lines_per_cell
Definition: geometry_info.h:2014
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
GeometryInfo::d_linear_shape_function_gradient
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
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)
std_cxx20::ranges::iota_view
boost::integer_range< IncrementableType > iota_view
Definition: iota_view.h:48
GeometryPrimitive::hex
@ hex
Definition: geometry_info.h:444
GeometryInfo::faces_per_cell
static constexpr unsigned int faces_per_cell
Definition: geometry_info.h:1933
GeometryInfo::subface_ratio
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
internal::SubfacePossibilities< 0 >::Possibilities
Possibilities
Definition: geometry_info.h:955
StandardExceptions::ExcImpossibleInDim
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
GeometryInfo::ucd_to_deal
static constexpr std::array< unsigned int, vertices_per_cell > ucd_to_deal
Definition: geometry_info.h:2053
GeometryPrimitive::object
Object object
Definition: geometry_info.h:470
RefinementCase::memory_consumption
static std::size_t memory_consumption()
RefinementPossibilities::Possibilities
Possibilities
Definition: geometry_info.h:525
GeometryInfo::quads_per_cell
static constexpr unsigned int quads_per_cell
Definition: geometry_info.h:2025
GeometryInfo::opposite_face
static constexpr std::array< unsigned int, faces_per_cell > opposite_face
Definition: geometry_info.h:2554
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:191
internal::SubfaceCase::value
std::uint8_t value
Definition: geometry_info.h:1216
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)
Point< dim >
GeometryPrimitive::line
@ line
Definition: geometry_info.h:436
config.h
RefinementPossibilities< 1 >::Possibilities
Possibilities
Definition: geometry_info.h:596
GeometryInfo::vertex_to_face
static constexpr std::array< std::array< unsigned int, dim >, vertices_per_cell > vertex_to_face
Definition: geometry_info.h:2083
internal
Definition: aligned_vector.h:369
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)
GeometryPrimitive::quad
@ quad
Definition: geometry_info.h:440
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)
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
GeometryInfo::distance_to_unit_cell
static double distance_to_unit_cell(const Point< dim > &p)
GeometryInfo::cell_to_child_coordinates
static Point< dim > cell_to_child_coordinates(const Point< dim > &p, const unsigned int child_index, const RefinementCase< dim > refine_case=RefinementCase< dim >::isotropic_refinement)
GeometryInfo::unit_normal_orientation
static constexpr std::array< int, faces_per_cell > unit_normal_orientation
Definition: geometry_info.h:2514
RefinementCase::ExcInvalidRefinementCase
static ::ExceptionBase & ExcInvalidRefinementCase(int arg1)
GeometryInfo::hexes_per_cell
static constexpr unsigned int hexes_per_cell
Definition: geometry_info.h:2032
int
point.h
internal::SubfacePossibilities< 3 >::Possibilities
Possibilities
Definition: geometry_info.h:1145
Generated by   doxygen 1.8.17