Reference documentation for deal.II version 9.4.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 - 2022 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/ndarray.h>
24 #include <deal.II/base/point.h>
25 #include <deal.II/base/std_cxx20/iota_view.h>
26 
27 #include <array>
28 #include <cstdint>
29 
30 
31 
32 DEAL_II_NAMESPACE_OPEN
33 
34 #ifndef DOXYGEN
35 namespace internal
36 {
37  namespace GeometryInfoHelper
38  {
39  // A struct that holds the values for all the arrays we want to initialize
40  // in GeometryInfo
41  template <int dim>
42  struct Initializers;
43 
44  template <>
45  struct Initializers<1>
46  {
47  static constexpr std::array<unsigned int, 2>
48  ucd_to_deal()
49  {
50  return {{0, 1}};
51  }
52 
53  static constexpr std::array<unsigned int, 2>
54  unit_normal_direction()
55  {
56  return {{0, 0}};
57  }
58 
59  static constexpr std::array<int, 2>
60  unit_normal_orientation()
61  {
62  return {{-1, 1}};
63  }
64 
65  static constexpr std::array<Tensor<1, 1>, 2>
66  unit_normal_vector()
67  {
68  return {{Tensor<1, 1>{{-1}}, Tensor<1, 1>{{1}}}};
69  }
70 
71  static constexpr ::ndarray<Tensor<1, 1>, 2, 0>
72  unit_tangential_vectors()
73  {
74  return {{{{}}, {{}}}};
75  }
76 
77  static constexpr std::array<unsigned int, 2>
78  opposite_face()
79  {
80  return {{1, 0}};
81  }
82 
83  static constexpr std::array<unsigned int, 2>
84  dx_to_deal()
85  {
86  return {{0, 1}};
87  }
88 
89  static constexpr ::ndarray<unsigned int, 2, 1>
90  vertex_to_face()
91  {
92  return {{{{0}}, {{1}}}};
93  }
94  };
95 
96  template <>
97  struct Initializers<2>
98  {
99  static constexpr std::array<unsigned int, 4>
100  ucd_to_deal()
101  {
102  return {{0, 1, 3, 2}};
103  }
104 
105  static constexpr std::array<unsigned int, 4>
106  unit_normal_direction()
107  {
108  return {{0, 0, 1, 1}};
109  }
110 
111  static constexpr std::array<int, 4>
112  unit_normal_orientation()
113  {
114  return {{-1, 1, -1, 1}};
115  }
116 
117  static constexpr std::array<Tensor<1, 2>, 4>
118  unit_normal_vector()
119  {
120  return {{Tensor<1, 2>{{-1., 0.}},
121  Tensor<1, 2>{{1., 0.}},
122  Tensor<1, 2>{{0., -1.}},
123  Tensor<1, 2>{{0., 1.}}}};
124  }
125 
126  static constexpr ::ndarray<Tensor<1, 2>, 4, 1>
127  unit_tangential_vectors()
128  {
129  return {{{{Tensor<1, 2>{{0, -1}}}},
130  {{Tensor<1, 2>{{0, 1}}}},
131  {{Tensor<1, 2>{{1, 0}}}},
132  {{Tensor<1, 2>{{-1, 0}}}}}};
133  }
134 
135  static constexpr std::array<unsigned int, 4>
136  opposite_face()
137  {
138  return {{1, 0, 3, 2}};
139  }
140 
141  static constexpr std::array<unsigned int, 4>
142  dx_to_deal()
143  {
144  return {{0, 2, 1, 3}};
145  }
146 
147  static constexpr ::ndarray<unsigned int, 4, 2>
148  vertex_to_face()
149  {
150  return {{{{0, 2}}, {{1, 2}}, {{0, 3}}, {{1, 3}}}};
151  }
152  };
153 
154  template <>
155  struct Initializers<3>
156  {
157  static constexpr std::array<unsigned int, 8>
158  ucd_to_deal()
159  {
160  return {{0, 1, 5, 4, 2, 3, 7, 6}};
161  }
162 
163  static constexpr std::array<unsigned int, 6>
164  unit_normal_direction()
165  {
166  return {{0, 0, 1, 1, 2, 2}};
167  }
168 
169  static constexpr std::array<int, 6>
170  unit_normal_orientation()
171  {
172  return {{-1, 1, -1, 1, -1, 1}};
173  }
174 
175  static constexpr std::array<Tensor<1, 3>, 6>
176  unit_normal_vector()
177  {
178  return {{Tensor<1, 3>{{-1, 0, 0}},
179  Tensor<1, 3>{{1, 0, 0}},
180  Tensor<1, 3>{{0, -1, 0}},
181  Tensor<1, 3>{{0, 1, 0}},
182  Tensor<1, 3>{{0, 0, -1}},
183  Tensor<1, 3>{{0, 0, 1}}}};
184  }
185 
186  static constexpr ::ndarray<Tensor<1, 3>, 6, 2>
187  unit_tangential_vectors()
188  {
189  return {{{{Tensor<1, 3>{{0, -1, 0}}, Tensor<1, 3>{{0, 0, 1}}}},
190  {{Tensor<1, 3>{{0, 1, 0}}, Tensor<1, 3>{{0, 0, 1}}}},
191  {{Tensor<1, 3>{{0, 0, -1}}, Tensor<1, 3>{{1, 0, 0}}}},
192  {{Tensor<1, 3>{{0, 0, 1}}, Tensor<1, 3>{{1, 0, 0}}}},
193  {{Tensor<1, 3>{{-1, 0, 0}}, Tensor<1, 3>{{0, 1, 0}}}},
194  {{Tensor<1, 3>{{1, 0, 0}}, Tensor<1, 3>{{0, 1, 0}}}}}};
195  }
196 
197  static constexpr std::array<unsigned int, 6>
198  opposite_face()
199  {
200  return {{1, 0, 3, 2, 5, 4}};
201  }
202 
203  static constexpr std::array<unsigned int, 8>
204  dx_to_deal()
205  {
206  return {{0, 4, 2, 6, 1, 5, 3, 7}};
207  }
208 
209  static constexpr ::ndarray<unsigned int, 8, 3>
210  vertex_to_face()
211  {
212  return {{{{0, 2, 4}},
213  {{1, 2, 4}},
214  {{0, 3, 4}},
215  {{1, 3, 4}},
216  {{0, 2, 5}},
217  {{1, 2, 5}},
218  {{0, 3, 5}},
219  {{1, 3, 5}}}};
220  }
221  };
222 
223  template <>
224  struct Initializers<4>
225  {
226  static constexpr std::array<unsigned int, 16>
227  ucd_to_deal()
228  {
229  return {{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  numbers::invalid_unsigned_int}};
245  }
246 
247  static constexpr std::array<unsigned int, 8>
248  unit_normal_direction()
249  {
250  return {{0, 0, 1, 1, 2, 2, 3, 3}};
251  }
252 
253  static constexpr std::array<int, 8>
254  unit_normal_orientation()
255  {
256  return {{-1, 1, -1, 1, -1, 1, -1, 1}};
257  }
258 
259  static constexpr std::array<Tensor<1, 4>, 8>
260  unit_normal_vector()
261  {
262  return {{Tensor<1, 4>{{-1, 0, 0, 0}},
263  Tensor<1, 4>{{1, 0, 0, 0}},
264  Tensor<1, 4>{{0, -1, 0, 0}},
265  Tensor<1, 4>{{0, 1, 0, 0}},
266  Tensor<1, 4>{{0, 0, -1, 0}},
267  Tensor<1, 4>{{0, 0, 1, 0}},
268  Tensor<1, 4>{{0, 0, 0, -1}},
269  Tensor<1, 4>{{0, 0, 0, 1}}}};
270  }
271 
272  static constexpr ::ndarray<Tensor<1, 4>, 8, 3>
273  unit_tangential_vectors()
274  {
275  return {{{{Tensor<1, 4>{{0, -1, 0, 0}},
276  Tensor<1, 4>{{0, 0, 1, 0}},
277  Tensor<1, 4>{{0, 0, 0, 1}}}},
278  {{Tensor<1, 4>{{0, 1, 0, 0}},
279  Tensor<1, 4>{{0, 0, 1, 0}},
280  Tensor<1, 4>{{0, 0, 0, 1}}}},
281  {{Tensor<1, 4>{{0, 0, -1, 0}},
282  Tensor<1, 4>{{0, 0, 0, 1}},
283  Tensor<1, 4>{{1, 0, 0, 0}}}},
284  {{Tensor<1, 4>{{0, 0, 1, 0}},
285  Tensor<1, 4>{{0, 0, 0, 1}},
286  Tensor<1, 4>{{1, 0, 0, 0}}}},
287  {{Tensor<1, 4>{{0, 0, 0, -1}},
288  Tensor<1, 4>{{1, 0, 0, 0}},
289  Tensor<1, 4>{{0, 1, 0, 0}}}},
290  {{Tensor<1, 4>{{0, 0, 0, 1}},
291  Tensor<1, 4>{{1, 0, 0, 0}},
292  Tensor<1, 4>{{0, 1, 0, 0}}}},
293  {{Tensor<1, 4>{{-1, 0, 0, 0}},
294  Tensor<1, 4>{{0, 1, 0, 0}},
295  Tensor<1, 4>{{0, 0, 1, 0}}}},
296  {{Tensor<1, 4>{{1, 0, 0, 0}},
297  Tensor<1, 4>{{0, 1, 0, 0}},
298  Tensor<1, 4>{{0, 0, 1, 0}}}}}};
299  }
300 
301  static constexpr std::array<unsigned int, 8>
302  opposite_face()
303  {
304  return {{1, 0, 3, 2, 5, 4, 7, 6}};
305  }
306 
307  static constexpr std::array<unsigned int, 16>
308  dx_to_deal()
309  {
310  return {{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  numbers::invalid_unsigned_int}};
326  }
327 
328  static constexpr ::ndarray<unsigned int, 16, 4>
329  vertex_to_face()
330  {
331  return {{{{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  numbers::invalid_unsigned_int}}}};
395  }
396  };
397  } // namespace GeometryInfoHelper
398 } // namespace internal
399 #endif // DOXYGEN
400 
401 
417 class GeometryPrimitive
418 {
419 public:
426  enum Object
427  {
431  vertex = 0,
435  line = 1,
439  quad = 2,
443  hex = 3
444  };
445 
450  GeometryPrimitive(const Object object);
451 
457  GeometryPrimitive(const unsigned int object_dimension);
458 
463  operator unsigned int() const;
464 
465 private:
469  Object object;
470 };
471 
472 
473 
486 template <int dim>
487 struct RefinementPossibilities
488 {
523  enum Possibilities
524  {
528  no_refinement = 0,
529 
541  isotropic_refinement = 0xFF
542  };
543 };
544 
545 
546 
556 template <>
557 struct RefinementPossibilities<1>
558 {
593  enum Possibilities
594  {
598  no_refinement = 0,
602  cut_x = 1,
606  isotropic_refinement = cut_x
607  };
608 };
609 
610 
611 
622 template <>
623 struct RefinementPossibilities<2>
624 {
659  enum Possibilities
660  {
664  no_refinement = 0,
668  cut_x = 1,
672  cut_y = 2,
676  cut_xy = cut_x | cut_y,
677 
681  isotropic_refinement = cut_xy
682  };
683 };
684 
685 
686 
697 template <>
698 struct RefinementPossibilities<3>
699 {
734  enum Possibilities
735  {
739  no_refinement = 0,
743  cut_x = 1,
747  cut_y = 2,
751  cut_xy = cut_x | cut_y,
755  cut_z = 4,
759  cut_xz = cut_x | cut_z,
763  cut_yz = cut_y | cut_z,
767  cut_xyz = cut_x | cut_y | cut_z,
768 
772  isotropic_refinement = cut_xyz
773  };
774 };
775 
776 
777 
788 template <int dim>
789 class RefinementCase : public RefinementPossibilities<dim>
790 {
791 public:
795  RefinementCase();
796 
801  RefinementCase(
802  const typename RefinementPossibilities<dim>::Possibilities refinement_case);
803 
809  explicit RefinementCase(const std::uint8_t refinement_case);
810 
822  operator std::uint8_t() const;
823 
828  RefinementCase
829  operator|(const RefinementCase &r) const;
830 
835  RefinementCase
836  operator&(const RefinementCase &r) const;
837 
845  RefinementCase
846  operator~() const;
847 
848 
854  static RefinementCase
855  cut_axis(const unsigned int i);
856 
860  static std::size_t
861  memory_consumption();
862 
868  template <class Archive>
869  void
870  serialize(Archive &ar, const unsigned int version);
871 
875  DeclException1(
876  ExcInvalidRefinementCase,
877  int,
878  << "The refinement flags given (" << arg1
879  << ") contain set bits that do not "
880  << "make sense for the space dimension of the object to which they are applied.");
881 
882 private:
887  std::uint8_t value : (dim > 0 ? dim : 1);
888 };
889 
890 
891 namespace internal
892 {
912  template <int dim>
913  struct SubfacePossibilities
914  {
918  enum Possibilities
919  {
923  case_none = 0,
924 
928  case_isotropic = static_cast<std::uint8_t>(-1)
929  };
930  };
931 
932 
941  template <>
942  struct SubfacePossibilities<0>
943  {
949  enum Possibilities
950  {
954  case_none = 0,
955 
959  case_isotropic = case_none
960  };
961  };
962 
963 
964 
974  template <>
975  struct SubfacePossibilities<1>
976  {
982  enum Possibilities
983  {
987  case_none = 0,
988 
992  case_isotropic = case_none
993  };
994  };
995 
996 
997 
1008  template <>
1009  struct SubfacePossibilities<2>
1010  {
1017  enum Possibilities
1018  {
1022  case_none = 0,
1026  case_x = 1,
1030  case_isotropic = case_x
1031  };
1032  };
1033 
1034 
1035 
1127  template <>
1128  struct SubfacePossibilities<3>
1129  {
1136  enum Possibilities
1137  {
1138  case_none = 0,
1139  case_x = 1,
1140  case_x1y = 2,
1141  case_x2y = 3,
1142  case_x1y2y = 4,
1143  case_y = 5,
1144  case_y1x = 6,
1145  case_y2x = 7,
1146  case_y1x2x = 8,
1147  case_xy = 9,
1148 
1149  case_isotropic = case_xy
1150  };
1151  };
1152 
1153 
1154 
1161  template <int dim>
1162  class SubfaceCase : public SubfacePossibilities<dim>
1163  {
1164  public:
1170  SubfaceCase(const typename SubfacePossibilities<dim>::Possibilities
1171  subface_possibility);
1172 
1184  operator std::uint8_t() const;
1185 
1189  static constexpr std::size_t
1190  memory_consumption();
1191 
1195  DeclException1(
1196  ExcInvalidSubfaceCase,
1197  int,
1198  << "The subface case given (" << arg1 << ") does not make sense "
1199  << "for the space dimension of the object to which they are applied.");
1200 
1201  private:
1206  std::uint8_t value : (dim == 3 ? 4 : 1);
1207  };
1208 
1209 } // namespace internal
1210 
1211 
1212 
1213 template <int dim>
1214 struct GeometryInfo;
1215 
1216 
1217 
1237 template <>
1238 struct GeometryInfo<0>
1239 {
1247  static constexpr unsigned int max_children_per_cell = 1;
1248 
1252  static constexpr unsigned int faces_per_cell = 0;
1253 
1270  static std::array<unsigned int, 0>
1271  face_indices();
1272 
1280  static constexpr unsigned int max_children_per_face = 0;
1281 
1287  static unsigned int
1288  n_children(const RefinementCase<0> &refinement_case);
1289 
1293  static constexpr unsigned int vertices_per_cell = 1;
1294 
1313  static std::array<unsigned int, vertices_per_cell>
1314  vertex_indices();
1315 
1339  static unsigned int
1340  face_to_cell_vertices(const unsigned int face,
1341  const unsigned int vertex,
1342  const bool face_orientation = true,
1343  const bool face_flip = false,
1344  const bool face_rotation = false);
1345 
1360  static unsigned int
1361  face_to_cell_lines(const unsigned int face,
1362  const unsigned int line,
1363  const bool face_orientation = true,
1364  const bool face_flip = false,
1365  const bool face_rotation = false);
1366 
1373  static constexpr unsigned int vertices_per_face = 0;
1374 
1378  static constexpr unsigned int lines_per_face = 0;
1379 
1383  static constexpr unsigned int quads_per_face = 0;
1384 
1388  static constexpr unsigned int lines_per_cell = 0;
1389 
1393  static constexpr unsigned int quads_per_cell = 0;
1394 
1398  static constexpr unsigned int hexes_per_cell = 0;
1399 
1417  static const std::array<unsigned int, vertices_per_cell> ucd_to_deal;
1418 
1432  static const std::array<unsigned int, vertices_per_cell> dx_to_deal;
1433 };
1434 
1435 
1436 
1967 template <int dim>
1968 struct GeometryInfo
1969 {
1977  static constexpr unsigned int max_children_per_cell = 1 << dim;
1978 
1982  static constexpr unsigned int faces_per_cell = 2 * dim;
1983 
2000  static std_cxx20::ranges::iota_view<unsigned int, unsigned int>
2001  face_indices();
2002 
2010  static constexpr unsigned int max_children_per_face =
2011  GeometryInfo<dim - 1>::max_children_per_cell;
2012 
2016  static constexpr unsigned int vertices_per_cell = 1 << dim;
2017 
2033  static std_cxx20::ranges::iota_view<unsigned int, unsigned int>
2034  vertex_indices();
2035 
2039  static constexpr unsigned int vertices_per_face =
2040  GeometryInfo<dim - 1>::vertices_per_cell;
2041 
2045  static constexpr unsigned int lines_per_face =
2046  GeometryInfo<dim - 1>::lines_per_cell;
2047 
2051  static constexpr unsigned int quads_per_face =
2052  GeometryInfo<dim - 1>::quads_per_cell;
2053 
2063  static constexpr unsigned int lines_per_cell =
2064  (2 * GeometryInfo<dim - 1>::lines_per_cell +
2065  GeometryInfo<dim - 1>::vertices_per_cell);
2066 
2074  static constexpr unsigned int quads_per_cell =
2075  (2 * GeometryInfo<dim - 1>::quads_per_cell +
2076  GeometryInfo<dim - 1>::lines_per_cell);
2077 
2081  static constexpr unsigned int hexes_per_cell =
2082  (2 * GeometryInfo<dim - 1>::hexes_per_cell +
2083  GeometryInfo<dim - 1>::quads_per_cell);
2084 
2102  static constexpr std::array<unsigned int, vertices_per_cell> ucd_to_deal =
2103  internal::GeometryInfoHelper::Initializers<dim>::ucd_to_deal();
2104 
2118  static constexpr std::array<unsigned int, vertices_per_cell> dx_to_deal =
2119  internal::GeometryInfoHelper::Initializers<dim>::dx_to_deal();
2120 
2131  static constexpr ndarray<unsigned int, vertices_per_cell, dim>
2132  vertex_to_face =
2133  internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face();
2134 
2139  static unsigned int
2140  n_children(const RefinementCase<dim> &refinement_case);
2141 
2146  static unsigned int
2147  n_subfaces(const internal::SubfaceCase<dim> &subface_case);
2148 
2158  static double
2159  subface_ratio(const internal::SubfaceCase<dim> &subface_case,
2160  const unsigned int subface_no);
2161 
2167  static RefinementCase<dim - 1>
2168  face_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2169  const unsigned int face_no,
2170  const bool face_orientation = true,
2171  const bool face_flip = false,
2172  const bool face_rotation = false);
2173 
2179  static RefinementCase<dim>
2180  min_cell_refinement_case_for_face_refinement(
2181  const RefinementCase<dim - 1> &face_refinement_case,
2182  const unsigned int face_no,
2183  const bool face_orientation = true,
2184  const bool face_flip = false,
2185  const bool face_rotation = false);
2186 
2191  static RefinementCase<1>
2192  line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
2193  const unsigned int line_no);
2194 
2199  static RefinementCase<dim>
2200  min_cell_refinement_case_for_line_refinement(const unsigned int line_no);
2201 
2248  static unsigned int
2249  child_cell_on_face(const RefinementCase<dim> & ref_case,
2250  const unsigned int face,
2251  const unsigned int subface,
2252  const bool face_orientation = true,
2253  const bool face_flip = false,
2254  const bool face_rotation = false,
2255  const RefinementCase<dim - 1> &face_refinement_case =
2256  RefinementCase<dim - 1>::isotropic_refinement);
2257 
2271  static unsigned int
2272  line_to_cell_vertices(const unsigned int line, const unsigned int vertex);
2273 
2294  static unsigned int
2295  face_to_cell_vertices(const unsigned int face,
2296  const unsigned int vertex,
2297  const bool face_orientation = true,
2298  const bool face_flip = false,
2299  const bool face_rotation = false);
2300 
2312  static unsigned int
2313  face_to_cell_lines(const unsigned int face,
2314  const unsigned int line,
2315  const bool face_orientation = true,
2316  const bool face_flip = false,
2317  const bool face_rotation = false);
2318 
2328  static unsigned int
2329  standard_to_real_face_vertex(const unsigned int vertex,
2330  const bool face_orientation = true,
2331  const bool face_flip = false,
2332  const bool face_rotation = false);
2333 
2343  static unsigned int
2344  real_to_standard_face_vertex(const unsigned int vertex,
2345  const bool face_orientation = true,
2346  const bool face_flip = false,
2347  const bool face_rotation = false);
2348 
2358  static unsigned int
2359  standard_to_real_face_line(const unsigned int line,
2360  const bool face_orientation = true,
2361  const bool face_flip = false,
2362  const bool face_rotation = false);
2363 
2369  static unsigned int
2370  standard_to_real_line_vertex(const unsigned int vertex,
2371  const bool line_orientation = true);
2372 
2380  static std::array<unsigned int, 2>
2381  standard_quad_vertex_to_line_vertex_index(const unsigned int vertex);
2382 
2390  static std::array<unsigned int, 2>
2391  standard_hex_vertex_to_quad_vertex_index(const unsigned int vertex);
2392 
2400  static std::array<unsigned int, 2>
2401  standard_hex_line_to_quad_line_index(const unsigned int line);
2402 
2412  static unsigned int
2413  real_to_standard_face_line(const unsigned int line,
2414  const bool face_orientation = true,
2415  const bool face_flip = false,
2416  const bool face_rotation = false);
2417 
2423  static Point<dim>
2424  unit_cell_vertex(const unsigned int vertex);
2425 
2435  static unsigned int
2436  child_cell_from_point(const Point<dim> &p);
2437 
2445  static Point<dim>
2446  cell_to_child_coordinates(const Point<dim> & p,
2447  const unsigned int child_index,
2448  const RefinementCase<dim> refine_case =
2449  RefinementCase<dim>::isotropic_refinement);
2450 
2456  static Point<dim>
2457  child_to_cell_coordinates(const Point<dim> & p,
2458  const unsigned int child_index,
2459  const RefinementCase<dim> refine_case =
2460  RefinementCase<dim>::isotropic_refinement);
2461 
2466  static bool
2467  is_inside_unit_cell(const Point<dim> &p);
2468 
2480  static bool
2481  is_inside_unit_cell(const Point<dim> &p, const double eps);
2482 
2487  template <typename Number = double>
2488  static Point<dim, Number>
2489  project_to_unit_cell(const Point<dim, Number> &p);
2490 
2496  static double
2497  distance_to_unit_cell(const Point<dim> &p);
2498 
2503  static double
2504  d_linear_shape_function(const Point<dim> &xi, const unsigned int i);
2505 
2510  static Tensor<1, dim>
2511  d_linear_shape_function_gradient(const Point<dim> &xi, const unsigned int i);
2512 
2564  template <int spacedim>
2565  static void
2566  alternating_form_at_vertices
2567 #ifndef DEAL_II_CXX14_CONSTEXPR_BUG
2568  (const Point<spacedim> (&vertices)[vertices_per_cell],
2569  Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell]);
2570 #else
2571  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms);
2572 #endif
2573 
2583  static constexpr std::array<unsigned int, faces_per_cell>
2584  unit_normal_direction =
2585  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_direction();
2586 
2603  static constexpr std::array<int, faces_per_cell> unit_normal_orientation =
2604  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_orientation();
2605 
2616  static constexpr std::array<Tensor<1, dim>, faces_per_cell>
2617  unit_normal_vector =
2618  internal::GeometryInfoHelper::Initializers<dim>::unit_normal_vector();
2619 
2633  static constexpr ndarray<Tensor<1, dim>, faces_per_cell, dim - 1>
2634  unit_tangential_vectors = internal::GeometryInfoHelper::Initializers<
2635  dim>::unit_tangential_vectors();
2636 
2642  static constexpr std::array<unsigned int, faces_per_cell> opposite_face =
2643  internal::GeometryInfoHelper::Initializers<dim>::opposite_face();
2644 
2645 
2649  DeclException1(ExcInvalidCoordinate,
2650  double,
2651  << "The coordinates must satisfy 0 <= x_i <= 1, "
2652  << "but here we have x_i=" << arg1);
2653 
2657  DeclException3(ExcInvalidSubface,
2658  int,
2659  int,
2660  int,
2661  << "RefinementCase<dim> " << arg1 << ": face " << arg2
2662  << " has no subface " << arg3);
2663 };
2664 
2665 
2666 
2667 #ifndef DOXYGEN
2668 
2669 
2670 /* -------------- declaration of explicit specializations ------------- */
2671 
2672 
2673 template <>
2674 Tensor<1, 1>
2675 GeometryInfo<1>::d_linear_shape_function_gradient(const Point<1> & xi,
2676  const unsigned int i);
2677 template <>
2678 Tensor<1, 2>
2679 GeometryInfo<2>::d_linear_shape_function_gradient(const Point<2> & xi,
2680  const unsigned int i);
2681 template <>
2682 Tensor<1, 3>
2683 GeometryInfo<3>::d_linear_shape_function_gradient(const Point<3> & xi,
2684  const unsigned int i);
2685 
2686 
2687 
2688 /* -------------- inline functions ------------- */
2689 
2690 
2691 inline GeometryPrimitive::GeometryPrimitive(const Object object)
2692  : object(object)
2693 {}
2694 
2695 
2696 
2697 inline GeometryPrimitive::GeometryPrimitive(const unsigned int object_dimension)
2698  : object(static_cast<Object>(object_dimension))
2699 {}
2700 
2701 
2702 inline GeometryPrimitive::operator unsigned int() const
2703 {
2704  return static_cast<unsigned int>(object);
2705 }
2706 
2707 
2708 
2709 namespace internal
2710 {
2711  template <int dim>
2712  inline SubfaceCase<dim>::SubfaceCase(
2713  const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
2714  : value(subface_possibility)
2715  {}
2716 
2717 
2718  template <int dim>
2719  inline SubfaceCase<dim>::operator std::uint8_t() const
2720  {
2721  return value;
2722  }
2723 
2724 
2725 } // namespace internal
2726 
2727 
2728 template <int dim>
2729 inline RefinementCase<dim>
2730 RefinementCase<dim>::cut_axis(const unsigned int)
2731 {
2732  Assert(false, ExcInternalError());
2733  return static_cast<std::uint8_t>(-1);
2734 }
2735 
2736 
2737 template <>
2738 inline RefinementCase<1>
2739 RefinementCase<1>::cut_axis(const unsigned int i)
2740 {
2741  AssertIndexRange(i, 1);
2742 
2743  const RefinementCase options[1] = {RefinementPossibilities<1>::cut_x};
2744  return options[i];
2745 }
2746 
2747 
2748 
2749 template <>
2750 inline RefinementCase<2>
2751 RefinementCase<2>::cut_axis(const unsigned int i)
2752 {
2753  AssertIndexRange(i, 2);
2754 
2755  const RefinementCase options[2] = {RefinementPossibilities<2>::cut_x,
2756  RefinementPossibilities<2>::cut_y};
2757  return options[i];
2758 }
2759 
2760 
2761 
2762 template <>
2763 inline RefinementCase<3>
2764 RefinementCase<3>::cut_axis(const unsigned int i)
2765 {
2766  AssertIndexRange(i, 3);
2767 
2768  const RefinementCase options[3] = {RefinementPossibilities<3>::cut_x,
2769  RefinementPossibilities<3>::cut_y,
2770  RefinementPossibilities<3>::cut_z};
2771  return options[i];
2772 }
2773 
2774 
2775 
2776 template <int dim>
2777 inline RefinementCase<dim>::RefinementCase()
2778  : value(RefinementPossibilities<dim>::no_refinement)
2779 {}
2780 
2781 
2782 
2783 template <int dim>
2784 inline RefinementCase<dim>::RefinementCase(
2785  const typename RefinementPossibilities<dim>::Possibilities refinement_case)
2786  : value(refinement_case)
2787 {
2788  // check that only those bits of
2789  // the given argument are set that
2790  // make sense for a given space
2791  // dimension
2792  Assert((refinement_case &
2793  RefinementPossibilities<dim>::isotropic_refinement) ==
2794  refinement_case,
2795  ExcInvalidRefinementCase(refinement_case));
2796 }
2797 
2798 
2799 
2800 template <int dim>
2801 inline RefinementCase<dim>::RefinementCase(const std::uint8_t refinement_case)
2802  : value(refinement_case)
2803 {
2804  // check that only those bits of
2805  // the given argument are set that
2806  // make sense for a given space
2807  // dimension
2808  Assert((refinement_case &
2809  RefinementPossibilities<dim>::isotropic_refinement) ==
2810  refinement_case,
2811  ExcInvalidRefinementCase(refinement_case));
2812 }
2813 
2814 
2815 
2816 template <int dim>
2817 inline RefinementCase<dim>::operator std::uint8_t() const
2818 {
2819  return value;
2820 }
2821 
2822 
2823 
2824 template <int dim>
2825 inline RefinementCase<dim>
2826 RefinementCase<dim>::operator|(const RefinementCase<dim> &r) const
2827 {
2828  return RefinementCase<dim>(static_cast<std::uint8_t>(value | r.value));
2829 }
2830 
2831 
2832 
2833 template <int dim>
2834 inline RefinementCase<dim>
2835 RefinementCase<dim>::operator&(const RefinementCase<dim> &r) const
2836 {
2837  return RefinementCase<dim>(static_cast<std::uint8_t>(value & r.value));
2838 }
2839 
2840 
2841 
2842 template <int dim>
2843 inline RefinementCase<dim>
2844 RefinementCase<dim>::operator~() const
2845 {
2846  return RefinementCase<dim>(static_cast<std::uint8_t>(
2847  (~value) & RefinementPossibilities<dim>::isotropic_refinement));
2848 }
2849 
2850 
2851 
2852 template <int dim>
2853 inline std::size_t
2854 RefinementCase<dim>::memory_consumption()
2855 {
2856  return sizeof(RefinementCase<dim>);
2857 }
2858 
2859 
2860 
2861 template <int dim>
2862 template <class Archive>
2863 inline void
2864 RefinementCase<dim>::serialize(Archive &ar, const unsigned int)
2865 {
2866  // serialization can't deal with bitfields, so copy from/to a full sized
2867  // std::uint8_t
2868  std::uint8_t uchar_value = value;
2869  ar & uchar_value;
2870  value = uchar_value;
2871 }
2872 
2873 
2874 
2875 template <>
2876 inline Point<1>
2877 GeometryInfo<1>::unit_cell_vertex(const unsigned int vertex)
2878 {
2879  AssertIndexRange(vertex, vertices_per_cell);
2880 
2881  return Point<1>(static_cast<double>(vertex));
2882 }
2883 
2884 
2885 
2886 template <>
2887 inline Point<2>
2888 GeometryInfo<2>::unit_cell_vertex(const unsigned int vertex)
2889 {
2890  AssertIndexRange(vertex, vertices_per_cell);
2891 
2892  return {static_cast<double>(vertex % 2), static_cast<double>(vertex / 2)};
2893 }
2894 
2895 
2896 
2897 template <>
2898 inline Point<3>
2899 GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
2900 {
2901  AssertIndexRange(vertex, vertices_per_cell);
2902 
2903  return {static_cast<double>(vertex % 2),
2904  static_cast<double>(vertex / 2 % 2),
2905  static_cast<double>(vertex / 4)};
2906 }
2907 
2908 
2909 
2910 inline std::array<unsigned int, 0>
2911 GeometryInfo<0>::face_indices()
2912 {
2913  return {{}};
2914 }
2915 
2916 
2917 
2918 inline std::array<unsigned int, 1>
2919 GeometryInfo<0>::vertex_indices()
2920 {
2921  return {{0}};
2922 }
2923 
2924 
2925 
2926 template <int dim>
2927 inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
2928 GeometryInfo<dim>::face_indices()
2929 {
2930  return {0U, faces_per_cell};
2931 }
2932 
2933 
2934 
2935 template <int dim>
2936 inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
2937 GeometryInfo<dim>::vertex_indices()
2938 {
2939  return {0U, vertices_per_cell};
2940 }
2941 
2942 
2943 
2944 template <int dim>
2945 inline Point<dim>
2946 GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
2947 {
2948  Assert(false, ExcNotImplemented());
2949 
2950  return {};
2951 }
2952 
2953 
2954 
2955 template <>
2956 inline unsigned int
2957 GeometryInfo<1>::child_cell_from_point(const Point<1> &p)
2958 {
2959  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2960 
2961  return (p[0] <= 0.5 ? 0 : 1);
2962 }
2963 
2964 
2965 
2966 template <>
2967 inline unsigned int
2968 GeometryInfo<2>::child_cell_from_point(const Point<2> &p)
2969 {
2970  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2971  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2972 
2973  return (p[0] <= 0.5 ? (p[1] <= 0.5 ? 0 : 2) : (p[1] <= 0.5 ? 1 : 3));
2974 }
2975 
2976 
2977 
2978 template <>
2979 inline unsigned int
2980 GeometryInfo<3>::child_cell_from_point(const Point<3> &p)
2981 {
2982  Assert((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
2983  Assert((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
2984  Assert((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
2985 
2986  return (p[0] <= 0.5 ?
2987  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 0 : 4) : (p[2] <= 0.5 ? 2 : 6)) :
2988  (p[1] <= 0.5 ? (p[2] <= 0.5 ? 1 : 5) : (p[2] <= 0.5 ? 3 : 7)));
2989 }
2990 
2991 
2992 template <int dim>
2993 inline unsigned int
2994 GeometryInfo<dim>::child_cell_from_point(const Point<dim> &)
2995 {
2996  Assert(false, ExcNotImplemented());
2997 
2998  return 0;
2999 }
3000 
3001 
3002 
3003 template <>
3004 inline Point<1>
3005 GeometryInfo<1>::cell_to_child_coordinates(const Point<1> & p,
3006  const unsigned int child_index,
3007  const RefinementCase<1> refine_case)
3008 
3009 {
3010  AssertIndexRange(child_index, 2);
3011  Assert(refine_case == RefinementCase<1>::cut_x, ExcInternalError());
3012  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3013 
3014  return Point<1>(p * 2.0 - unit_cell_vertex(child_index));
3015 }
3016 
3017 
3018 
3019 template <>
3020 inline Point<2>
3021 GeometryInfo<2>::cell_to_child_coordinates(const Point<2> & p,
3022  const unsigned int child_index,
3023  const RefinementCase<2> refine_case)
3024 
3025 {
3026  AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3027 
3028  Point<2> point = p;
3029  switch (refine_case)
3030  {
3031  case RefinementCase<2>::cut_x:
3032  point[0] *= 2.0;
3033  if (child_index == 1)
3034  point[0] -= 1.0;
3035  break;
3036  case RefinementCase<2>::cut_y:
3037  point[1] *= 2.0;
3038  if (child_index == 1)
3039  point[1] -= 1.0;
3040  break;
3041  case RefinementCase<2>::cut_xy:
3042  point *= 2.0;
3043  point -= unit_cell_vertex(child_index);
3044  break;
3045  default:
3046  Assert(false, ExcInternalError());
3047  }
3048 
3049  return point;
3050 }
3051 
3052 
3053 
3054 template <>
3055 inline Point<3>
3056 GeometryInfo<3>::cell_to_child_coordinates(const Point<3> & p,
3057  const unsigned int child_index,
3058  const RefinementCase<3> refine_case)
3059 
3060 {
3061  AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3062 
3063  Point<3> point = p;
3064  // there might be a cleverer way to do
3065  // this, but since this function is called
3066  // in very few places for initialization
3067  // purposes only, I don't care at the
3068  // moment
3069  switch (refine_case)
3070  {
3071  case RefinementCase<3>::cut_x:
3072  point[0] *= 2.0;
3073  if (child_index == 1)
3074  point[0] -= 1.0;
3075  break;
3076  case RefinementCase<3>::cut_y:
3077  point[1] *= 2.0;
3078  if (child_index == 1)
3079  point[1] -= 1.0;
3080  break;
3081  case RefinementCase<3>::cut_z:
3082  point[2] *= 2.0;
3083  if (child_index == 1)
3084  point[2] -= 1.0;
3085  break;
3086  case RefinementCase<3>::cut_xy:
3087  point[0] *= 2.0;
3088  point[1] *= 2.0;
3089  if (child_index % 2 == 1)
3090  point[0] -= 1.0;
3091  if (child_index / 2 == 1)
3092  point[1] -= 1.0;
3093  break;
3094  case RefinementCase<3>::cut_xz:
3095  // careful, this is slightly
3096  // different from xy and yz due to
3097  // different internal numbering of
3098  // children!
3099  point[0] *= 2.0;
3100  point[2] *= 2.0;
3101  if (child_index / 2 == 1)
3102  point[0] -= 1.0;
3103  if (child_index % 2 == 1)
3104  point[2] -= 1.0;
3105  break;
3106  case RefinementCase<3>::cut_yz:
3107  point[1] *= 2.0;
3108  point[2] *= 2.0;
3109  if (child_index % 2 == 1)
3110  point[1] -= 1.0;
3111  if (child_index / 2 == 1)
3112  point[2] -= 1.0;
3113  break;
3114  case RefinementCase<3>::cut_xyz:
3115  point *= 2.0;
3116  point -= unit_cell_vertex(child_index);
3117  break;
3118  default:
3119  Assert(false, ExcInternalError());
3120  }
3121 
3122  return point;
3123 }
3124 
3125 
3126 
3127 template <int dim>
3128 inline Point<dim>
3129 GeometryInfo<dim>::cell_to_child_coordinates(
3130  const Point<dim> & /*p*/,
3131  const unsigned int /*child_index*/,
3132  const RefinementCase<dim> /*refine_case*/)
3133 
3134 {
3135  Assert(false, ExcNotImplemented());
3136  return {};
3137 }
3138 
3139 
3140 
3141 template <>
3142 inline Point<1>
3143 GeometryInfo<1>::child_to_cell_coordinates(const Point<1> & p,
3144  const unsigned int child_index,
3145  const RefinementCase<1> refine_case)
3146 
3147 {
3148  AssertIndexRange(child_index, 2);
3149  Assert(refine_case == RefinementCase<1>::cut_x, ExcInternalError());
3150  (void)refine_case; // removes -Wunused-parameter warning in optimized mode
3151 
3152  return (p + unit_cell_vertex(child_index)) * 0.5;
3153 }
3154 
3155 
3156 
3157 template <>
3158 inline Point<3>
3159 GeometryInfo<3>::child_to_cell_coordinates(const Point<3> & p,
3160  const unsigned int child_index,
3161  const RefinementCase<3> refine_case)
3162 
3163 {
3164  AssertIndexRange(child_index, GeometryInfo<3>::n_children(refine_case));
3165 
3166  Point<3> point = p;
3167  // there might be a cleverer way to do
3168  // this, but since this function is called
3169  // in very few places for initialization
3170  // purposes only, I don't care at the
3171  // moment
3172  switch (refine_case)
3173  {
3174  case RefinementCase<3>::cut_x:
3175  if (child_index == 1)
3176  point[0] += 1.0;
3177  point[0] *= 0.5;
3178  break;
3179  case RefinementCase<3>::cut_y:
3180  if (child_index == 1)
3181  point[1] += 1.0;
3182  point[1] *= 0.5;
3183  break;
3184  case RefinementCase<3>::cut_z:
3185  if (child_index == 1)
3186  point[2] += 1.0;
3187  point[2] *= 0.5;
3188  break;
3189  case RefinementCase<3>::cut_xy:
3190  if (child_index % 2 == 1)
3191  point[0] += 1.0;
3192  if (child_index / 2 == 1)
3193  point[1] += 1.0;
3194  point[0] *= 0.5;
3195  point[1] *= 0.5;
3196  break;
3197  case RefinementCase<3>::cut_xz:
3198  // careful, this is slightly
3199  // different from xy and yz due to
3200  // different internal numbering of
3201  // children!
3202  if (child_index / 2 == 1)
3203  point[0] += 1.0;
3204  if (child_index % 2 == 1)
3205  point[2] += 1.0;
3206  point[0] *= 0.5;
3207  point[2] *= 0.5;
3208  break;
3209  case RefinementCase<3>::cut_yz:
3210  if (child_index % 2 == 1)
3211  point[1] += 1.0;
3212  if (child_index / 2 == 1)
3213  point[2] += 1.0;
3214  point[1] *= 0.5;
3215  point[2] *= 0.5;
3216  break;
3217  case RefinementCase<3>::cut_xyz:
3218  point += unit_cell_vertex(child_index);
3219  point *= 0.5;
3220  break;
3221  default:
3222  Assert(false, ExcInternalError());
3223  }
3224 
3225  return point;
3226 }
3227 
3228 
3229 
3230 template <>
3231 inline Point<2>
3232 GeometryInfo<2>::child_to_cell_coordinates(const Point<2> & p,
3233  const unsigned int child_index,
3234  const RefinementCase<2> refine_case)
3235 {
3236  AssertIndexRange(child_index, GeometryInfo<2>::n_children(refine_case));
3237 
3238  Point<2> point = p;
3239  switch (refine_case)
3240  {
3241  case RefinementCase<2>::cut_x:
3242  if (child_index == 1)
3243  point[0] += 1.0;
3244  point[0] *= 0.5;
3245  break;
3246  case RefinementCase<2>::cut_y:
3247  if (child_index == 1)
3248  point[1] += 1.0;
3249  point[1] *= 0.5;
3250  break;
3251  case RefinementCase<2>::cut_xy:
3252  point += unit_cell_vertex(child_index);
3253  point *= 0.5;
3254  break;
3255  default:
3256  Assert(false, ExcInternalError());
3257  }
3258 
3259  return point;
3260 }
3261 
3262 
3263 
3264 template <int dim>
3265 inline Point<dim>
3266 GeometryInfo<dim>::child_to_cell_coordinates(
3267  const Point<dim> & /*p*/,
3268  const unsigned int /*child_index*/,
3269  const RefinementCase<dim> /*refine_case*/)
3270 {
3271  Assert(false, ExcNotImplemented());
3272  return {};
3273 }
3274 
3275 
3276 
3277 template <int dim>
3278 inline bool
3279 GeometryInfo<dim>::is_inside_unit_cell(const Point<dim> &)
3280 {
3281  Assert(false, ExcNotImplemented());
3282  return false;
3283 }
3284 
3285 template <>
3286 inline bool
3287 GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p)
3288 {
3289  return (p[0] >= 0.) && (p[0] <= 1.);
3290 }
3291 
3292 
3293 
3294 template <>
3295 inline bool
3296 GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p)
3297 {
3298  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.);
3299 }
3300 
3301 
3302 
3303 template <>
3304 inline bool
3305 GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p)
3306 {
3307  return (p[0] >= 0.) && (p[0] <= 1.) && (p[1] >= 0.) && (p[1] <= 1.) &&
3308  (p[2] >= 0.) && (p[2] <= 1.);
3309 }
3310 
3311 
3312 
3313 template <int dim>
3314 inline bool
3315 GeometryInfo<dim>::is_inside_unit_cell(const Point<dim> &, const double)
3316 {
3317  Assert(false, ExcNotImplemented());
3318  return false;
3319 }
3320 
3321 template <>
3322 inline bool
3323 GeometryInfo<1>::is_inside_unit_cell(const Point<1> &p, const double eps)
3324 {
3325  return (p[0] >= -eps) && (p[0] <= 1. + eps);
3326 }
3327 
3328 
3329 
3330 template <>
3331 inline bool
3332 GeometryInfo<2>::is_inside_unit_cell(const Point<2> &p, const double eps)
3333 {
3334  const double l = -eps, u = 1 + eps;
3335  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u);
3336 }
3337 
3338 
3339 
3340 template <>
3341 inline bool
3342 GeometryInfo<3>::is_inside_unit_cell(const Point<3> &p, const double eps)
3343 {
3344  const double l = -eps, u = 1.0 + eps;
3345  return (p[0] >= l) && (p[0] <= u) && (p[1] >= l) && (p[1] <= u) &&
3346  (p[2] >= l) && (p[2] <= u);
3347 }
3348 
3349 
3350 
3351 template <>
3352 inline unsigned int
3353 GeometryInfo<1>::line_to_cell_vertices(const unsigned int line,
3354  const unsigned int vertex)
3355 {
3356  (void)line;
3357  AssertIndexRange(line, lines_per_cell);
3358  AssertIndexRange(vertex, 2);
3359 
3360  return vertex;
3361 }
3362 
3363 
3364 template <>
3365 inline unsigned int
3366 GeometryInfo<2>::line_to_cell_vertices(const unsigned int line,
3367  const unsigned int vertex)
3368 {
3369  constexpr unsigned int cell_vertices[4][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}};
3370  return cell_vertices[line][vertex];
3371 }
3372 
3373 
3374 
3375 template <>
3376 inline unsigned int
3377 GeometryInfo<3>::line_to_cell_vertices(const unsigned int line,
3378  const unsigned int vertex)
3379 {
3380  AssertIndexRange(line, lines_per_cell);
3381  AssertIndexRange(vertex, 2);
3382 
3383  constexpr unsigned vertices[lines_per_cell][2] = {{0, 2}, // bottom face
3384  {1, 3},
3385  {0, 1},
3386  {2, 3},
3387  {4, 6}, // top face
3388  {5, 7},
3389  {4, 5},
3390  {6, 7},
3391  {0,
3392  4}, // connects of bottom
3393  {1, 5}, // top face
3394  {2, 6},
3395  {3, 7}};
3396  return vertices[line][vertex];
3397 }
3398 
3399 
3400 template <>
3401 inline unsigned int
3402 GeometryInfo<4>::line_to_cell_vertices(const unsigned int, const unsigned int)
3403 {
3404  Assert(false, ExcNotImplemented());
3405  return numbers::invalid_unsigned_int;
3406 }
3407 
3408 template <>
3409 inline unsigned int
3410 GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
3411  const bool face_orientation,
3412  const bool face_flip,
3413  const bool face_rotation)
3414 {
3415  AssertIndexRange(vertex, GeometryInfo<3>::vertices_per_face);
3416 
3417  // set up a table to make sure that
3418  // we handle non-standard faces correctly
3419  //
3420  // so set up a table that for each vertex (of
3421  // a quad in standard position) describes
3422  // which vertex to take
3423  //
3424  // first index: four vertices 0...3
3425  //
3426  // second index: face_orientation; 0:
3427  // opposite normal, 1: standard
3428  //
3429  // third index: face_flip; 0: standard, 1:
3430  // face rotated by 180 degrees
3431  //
3432  // forth index: face_rotation: 0: standard,
3433  // 1: face rotated by 90 degrees
3434 
3435  constexpr unsigned int vertex_translation[4][2][2][2] = {
3436  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
3437  // face_rotation=false and true
3438  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
3439  // face_rotation=false and true
3440  {{0, 2}, // vertex 0, face_orientation=true, face_flip=false,
3441  // face_rotation=false and true
3442  {3, 1}}}, // vertex 0, face_orientation=true, face_flip=true,
3443  // face_rotation=false and true
3444 
3445  {{{2, 3}, // vertex 1 ...
3446  {1, 0}},
3447  {{1, 0}, {2, 3}}},
3448 
3449  {{{1, 0}, // vertex 2 ...
3450  {2, 3}},
3451  {{2, 3}, {1, 0}}},
3452 
3453  {{{3, 1}, // vertex 3 ...
3454  {0, 2}},
3455  {{3, 1}, {0, 2}}}};
3456 
3457  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
3458 }
3459 
3460 
3461 
3462 template <int dim>
3463 inline unsigned int
3464 GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
3465  const bool,
3466  const bool,
3467  const bool)
3468 {
3469  Assert(dim > 1, ExcImpossibleInDim(dim));
3470  AssertIndexRange(vertex, GeometryInfo<dim>::vertices_per_face);
3471  return vertex;
3472 }
3473 
3474 template <int dim>
3475 inline unsigned int
3476 GeometryInfo<dim>::n_children(const RefinementCase<dim> &ref_case)
3477 {
3478  constexpr unsigned int n_children[RefinementCase<3>::cut_xyz + 1] = {
3479  0, 2, 2, 4, 2, 4, 4, 8};
3480 
3481  return n_children[ref_case];
3482 }
3483 
3484 
3485 
3486 template <int dim>
3487 inline unsigned int
3488 GeometryInfo<dim>::n_subfaces(const internal::SubfaceCase<dim> &)
3489 {
3490  Assert(false, ExcNotImplemented());
3491  return 0;
3492 }
3493 
3494 template <>
3495 inline unsigned int
3496 GeometryInfo<1>::n_subfaces(const internal::SubfaceCase<1> &)
3497 {
3498  Assert(false, ExcImpossibleInDim(1));
3499  return 0;
3500 }
3501 
3502 template <>
3503 inline unsigned int
3504 GeometryInfo<2>::n_subfaces(const internal::SubfaceCase<2> &subface_case)
3505 {
3506  return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
3507 }
3508 
3509 
3510 
3511 template <>
3512 inline unsigned int
3513 GeometryInfo<3>::n_subfaces(const internal::SubfaceCase<3> &subface_case)
3514 {
3515  const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic + 1] = {
3516  0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
3517  return nsubs[subface_case];
3518 }
3519 
3520 
3521 
3522 template <int dim>
3523 inline double
3524 GeometryInfo<dim>::subface_ratio(const internal::SubfaceCase<dim> &,
3525  const unsigned int)
3526 {
3527  Assert(false, ExcNotImplemented());
3528  return 0.;
3529 }
3530 
3531 template <>
3532 inline double
3533 GeometryInfo<1>::subface_ratio(const internal::SubfaceCase<1> &,
3534  const unsigned int)
3535 {
3536  return 1;
3537 }
3538 
3539 
3540 template <>
3541 inline double
3542 GeometryInfo<2>::subface_ratio(const internal::SubfaceCase<2> &subface_case,
3543  const unsigned int)
3544 {
3545  double ratio = 1;
3546  switch (subface_case)
3547  {
3548  case internal::SubfaceCase<2>::case_none:
3549  // Here, an
3550  // Assert(false,ExcInternalError())
3551  // would be the right
3552  // choice, but
3553  // unfortunately the
3554  // current function is
3555  // also called for faces
3556  // without children (see
3557  // tests/fe/mapping.cc).
3558  // Assert(false, ExcMessage("Face has no subfaces."));
3559  // Furthermore, assign
3560  // following value as
3561  // otherwise the
3562  // bits/volume_x tests
3563  // break
3564  ratio = 1. / GeometryInfo<2>::max_children_per_face;
3565  break;
3566  case internal::SubfaceCase<2>::case_x:
3567  ratio = 0.5;
3568  break;
3569  default:
3570  // there should be no
3571  // cases left
3572  Assert(false, ExcInternalError());
3573  break;
3574  }
3575 
3576  return ratio;
3577 }
3578 
3579 
3580 template <>
3581 inline double
3582 GeometryInfo<3>::subface_ratio(const internal::SubfaceCase<3> &subface_case,
3583  const unsigned int subface_no)
3584 {
3585  double ratio = 1;
3586  switch (subface_case)
3587  {
3588  case internal::SubfaceCase<3>::case_none:
3589  // Here, an
3590  // Assert(false,ExcInternalError())
3591  // would be the right
3592  // choice, but
3593  // unfortunately the
3594  // current function is
3595  // also called for faces
3596  // without children (see
3597  // tests/bits/mesh_3d_16.cc). Add
3598  // following switch to
3599  // avoid diffs in
3600  // tests/bits/mesh_3d_16
3601  ratio = 1. / GeometryInfo<3>::max_children_per_face;
3602  break;
3603  case internal::SubfaceCase<3>::case_x:
3604  case internal::SubfaceCase<3>::case_y:
3605  ratio = 0.5;
3606  break;
3607  case internal::SubfaceCase<3>::case_xy:
3608  case internal::SubfaceCase<3>::case_x1y2y:
3609  case internal::SubfaceCase<3>::case_y1x2x:
3610  ratio = 0.25;
3611  break;
3612  case internal::SubfaceCase<3>::case_x1y:
3613  case internal::SubfaceCase<3>::case_y1x:
3614  if (subface_no < 2)
3615  ratio = 0.25;
3616  else
3617  ratio = 0.5;
3618  break;
3619  case internal::SubfaceCase<3>::case_x2y:
3620  case internal::SubfaceCase<3>::case_y2x:
3621  if (subface_no == 0)
3622  ratio = 0.5;
3623  else
3624  ratio = 0.25;
3625  break;
3626  default:
3627  // there should be no
3628  // cases left
3629  Assert(false, ExcInternalError());
3630  break;
3631  }
3632 
3633  return ratio;
3634 }
3635 
3636 
3637 
3638 template <int dim>
3639 RefinementCase<dim - 1> inline GeometryInfo<dim>::face_refinement_case(
3640  const RefinementCase<dim> &,
3641  const unsigned int,
3642  const bool,
3643  const bool,
3644  const bool)
3645 {
3646  Assert(false, ExcNotImplemented());
3647  return RefinementCase<dim - 1>::no_refinement;
3648 }
3649 
3650 template <>
3651 RefinementCase<0> inline GeometryInfo<1>::face_refinement_case(
3652  const RefinementCase<1> &,
3653  const unsigned int,
3654  const bool,
3655  const bool,
3656  const bool)
3657 {
3658  Assert(false, ExcImpossibleInDim(1));
3659 
3660  return RefinementCase<0>::no_refinement;
3661 }
3662 
3663 
3664 template <>
3665 inline RefinementCase<1>
3666 GeometryInfo<2>::face_refinement_case(
3667  const RefinementCase<2> &cell_refinement_case,
3668  const unsigned int face_no,
3669  const bool,
3670  const bool,
3671  const bool)
3672 {
3673  const unsigned int dim = 2;
3674  AssertIndexRange(cell_refinement_case,
3675  RefinementCase<dim>::isotropic_refinement + 1);
3676  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
3677 
3678  const RefinementCase<dim - 1>
3679  ref_cases[RefinementCase<dim>::isotropic_refinement +
3680  1][GeometryInfo<dim>::faces_per_cell / 2] = {
3681  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3682  RefinementCase<dim - 1>::no_refinement},
3683 
3684  {RefinementCase<dim - 1>::no_refinement, RefinementCase<dim - 1>::cut_x},
3685 
3686  {RefinementCase<dim - 1>::cut_x, RefinementCase<dim - 1>::no_refinement},
3687 
3688  {RefinementCase<dim - 1>::cut_x, // cut_xy
3689  RefinementCase<dim - 1>::cut_x}};
3690 
3691  return ref_cases[cell_refinement_case][face_no / 2];
3692 }
3693 
3694 
3695 template <>
3696 inline RefinementCase<2>
3697 GeometryInfo<3>::face_refinement_case(
3698  const RefinementCase<3> &cell_refinement_case,
3699  const unsigned int face_no,
3700  const bool face_orientation,
3701  const bool /*face_flip*/,
3702  const bool face_rotation)
3703 {
3704  const unsigned int dim = 3;
3705  AssertIndexRange(cell_refinement_case,
3706  RefinementCase<dim>::isotropic_refinement + 1);
3707  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
3708 
3709  const RefinementCase<dim - 1>
3710  ref_cases[RefinementCase<dim>::isotropic_refinement + 1]
3711  [GeometryInfo<dim>::faces_per_cell / 2] = {
3712  {RefinementCase<dim - 1>::no_refinement, // no_refinement
3713  RefinementCase<dim - 1>::no_refinement,
3714  RefinementCase<dim - 1>::no_refinement},
3715 
3716  {RefinementCase<dim - 1>::no_refinement, // cut_x
3717  RefinementCase<dim - 1>::cut_y,
3718  RefinementCase<dim - 1>::cut_x},
3719 
3720  {RefinementCase<dim - 1>::cut_x, // cut_y
3721  RefinementCase<dim - 1>::no_refinement,
3722  RefinementCase<dim - 1>::cut_y},
3723 
3724  {RefinementCase<dim - 1>::cut_x, // cut_xy
3725  RefinementCase<dim - 1>::cut_y,
3726  RefinementCase<dim - 1>::cut_xy},
3727 
3728  {RefinementCase<dim - 1>::cut_y, // cut_z
3729  RefinementCase<dim - 1>::cut_x,
3730  RefinementCase<dim - 1>::no_refinement},
3731 
3732  {RefinementCase<dim - 1>::cut_y, // cut_xz
3733  RefinementCase<dim - 1>::cut_xy,
3734  RefinementCase<dim - 1>::cut_x},
3735 
3736  {RefinementCase<dim - 1>::cut_xy, // cut_yz
3737  RefinementCase<dim - 1>::cut_x,
3738  RefinementCase<dim - 1>::cut_y},
3739 
3740  {RefinementCase<dim - 1>::cut_xy, // cut_xyz
3741  RefinementCase<dim - 1>::cut_xy,
3742  RefinementCase<dim - 1>::cut_xy},
3743  };
3744 
3745  const RefinementCase<dim - 1> ref_case =
3746  ref_cases[cell_refinement_case][face_no / 2];
3747 
3748  const RefinementCase<dim - 1> flip[4] = {
3749  RefinementCase<dim - 1>::no_refinement,
3750  RefinementCase<dim - 1>::cut_y,
3751  RefinementCase<dim - 1>::cut_x,
3752  RefinementCase<dim - 1>::cut_xy};
3753 
3754  // correct the ref_case for face_orientation
3755  // and face_rotation. for face_orientation,
3756  // 'true' is the default value whereas for
3757  // face_rotation, 'false' is standard. If
3758  // <tt>face_rotation==face_orientation</tt>,
3759  // then one of them is non-standard and we
3760  // have to swap cut_x and cut_y, otherwise no
3761  // change is necessary. face_flip has no
3762  // influence. however, in order to keep the
3763  // interface consistent with other functions,
3764  // we still include it as an argument to this
3765  // function
3766  return (face_orientation == face_rotation) ? flip[ref_case] : ref_case;
3767 }
3768 
3769 
3770 
3771 template <int dim>
3772 inline RefinementCase<1>
3773 GeometryInfo<dim>::line_refinement_case(const RefinementCase<dim> &,
3774  const unsigned int)
3775 {
3776  Assert(false, ExcNotImplemented());
3777  return RefinementCase<1>::no_refinement;
3778 }
3779 
3780 template <>
3781 inline RefinementCase<1>
3782 GeometryInfo<1>::line_refinement_case(
3783  const RefinementCase<1> &cell_refinement_case,
3784  const unsigned int line_no)
3785 {
3786  (void)line_no;
3787  const unsigned int dim = 1;
3788  (void)dim;
3789  AssertIndexRange(cell_refinement_case,
3790  RefinementCase<dim>::isotropic_refinement + 1);
3791  AssertIndexRange(line_no, GeometryInfo<dim>::lines_per_cell);
3792 
3793  return cell_refinement_case;
3794 }
3795 
3796 
3797 template <>
3798 inline RefinementCase<1>
3799 GeometryInfo<2>::line_refinement_case(
3800  const RefinementCase<2> &cell_refinement_case,
3801  const unsigned int line_no)
3802 {
3803  // Assertions are in face_refinement_case()
3804  return face_refinement_case(cell_refinement_case, line_no);
3805 }
3806 
3807 
3808 template <>
3809 inline RefinementCase<1>
3810 GeometryInfo<3>::line_refinement_case(
3811  const RefinementCase<3> &cell_refinement_case,
3812  const unsigned int line_no)
3813 {
3814  const unsigned int dim = 3;
3815  AssertIndexRange(cell_refinement_case,
3816  RefinementCase<dim>::isotropic_refinement + 1);
3817  AssertIndexRange(line_no, GeometryInfo<dim>::lines_per_cell);
3818 
3819  // array indicating, which simple refine
3820  // case cuts a line in direction x, y or
3821  // z. For example, cut_y and everything
3822  // containing cut_y (cut_xy, cut_yz,
3823  // cut_xyz) cuts lines, which are in y
3824  // direction.
3825  const RefinementCase<dim> cut_one[dim] = {RefinementCase<dim>::cut_x,
3826  RefinementCase<dim>::cut_y,
3827  RefinementCase<dim>::cut_z};
3828 
3829  // order the direction of lines
3830  // 0->x, 1->y, 2->z
3831  const unsigned int direction[lines_per_cell] = {
3832  1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3833 
3834  return ((cell_refinement_case & cut_one[direction[line_no]]) ?
3835  RefinementCase<1>::cut_x :
3836  RefinementCase<1>::no_refinement);
3837 }
3838 
3839 
3840 
3841 template <int dim>
3842 inline RefinementCase<dim>
3843 GeometryInfo<dim>::min_cell_refinement_case_for_face_refinement(
3844  const RefinementCase<dim - 1> &,
3845  const unsigned int,
3846  const bool,
3847  const bool,
3848  const bool)
3849 {
3850  Assert(false, ExcNotImplemented());
3851 
3852  return RefinementCase<dim>::no_refinement;
3853 }
3854 
3855 template <>
3856 inline RefinementCase<1>
3857 GeometryInfo<1>::min_cell_refinement_case_for_face_refinement(
3858  const RefinementCase<0> &,
3859  const unsigned int,
3860  const bool,
3861  const bool,
3862  const bool)
3863 {
3864  const unsigned int dim = 1;
3865  Assert(false, ExcImpossibleInDim(dim));
3866 
3867  return RefinementCase<dim>::no_refinement;
3868 }
3869 
3870 
3871 template <>
3872 inline RefinementCase<2>
3873 GeometryInfo<2>::min_cell_refinement_case_for_face_refinement(
3874  const RefinementCase<1> &face_refinement_case,
3875  const unsigned int face_no,
3876  const bool,
3877  const bool,
3878  const bool)
3879 {
3880  const unsigned int dim = 2;
3881  AssertIndexRange(face_refinement_case,
3882  RefinementCase<dim - 1>::isotropic_refinement + 1);
3883  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
3884 
3885  if (face_refinement_case == RefinementCase<dim>::cut_x)
3886  return (face_no / 2) != 0u ? RefinementCase<dim>::cut_x :
3887  RefinementCase<dim>::cut_y;
3888  else
3889  return RefinementCase<dim>::no_refinement;
3890 }
3891 
3892 
3893 template <>
3894 inline RefinementCase<3>
3895 GeometryInfo<3>::min_cell_refinement_case_for_face_refinement(
3896  const RefinementCase<2> &face_refinement_case,
3897  const unsigned int face_no,
3898  const bool face_orientation,
3899  const bool /*face_flip*/,
3900  const bool face_rotation)
3901 {
3902  const unsigned int dim = 3;
3903  AssertIndexRange(face_refinement_case,
3904  RefinementCase<dim - 1>::isotropic_refinement + 1);
3905  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
3906 
3907  const RefinementCase<2> flip[4] = {RefinementCase<2>::no_refinement,
3908  RefinementCase<2>::cut_y,
3909  RefinementCase<2>::cut_x,
3910  RefinementCase<2>::cut_xy};
3911 
3912  // correct the face_refinement_case for
3913  // face_orientation and face_rotation. for
3914  // face_orientation, 'true' is the default
3915  // value whereas for face_rotation, 'false'
3916  // is standard. If
3917  // <tt>face_rotation==face_orientation</tt>,
3918  // then one of them is non-standard and we
3919  // have to swap cut_x and cut_y, otherwise no
3920  // change is necessary. face_flip has no
3921  // influence. however, in order to keep the
3922  // interface consistent with other functions,
3923  // we still include it as an argument to this
3924  // function
3925  const RefinementCase<dim - 1> std_face_ref =
3926  (face_orientation == face_rotation) ? flip[face_refinement_case] :
3927  face_refinement_case;
3928 
3929  const RefinementCase<dim> face_to_cell[3][4] = {
3930  {RefinementCase<dim>::no_refinement, // faces 0 and 1
3931  RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
3932  RefinementCase<dim>::cut_z,
3933  RefinementCase<dim>::cut_yz},
3934 
3935  {RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are
3936  // "exchanged on faces 2 and 3")
3937  RefinementCase<dim>::cut_z,
3938  RefinementCase<dim>::cut_x,
3939  RefinementCase<dim>::cut_xz},
3940 
3941  {RefinementCase<dim>::no_refinement, // faces 4 and 5
3942  RefinementCase<dim>::cut_x,
3943  RefinementCase<dim>::cut_y,
3944  RefinementCase<dim>::cut_xy}};
3945 
3946  return face_to_cell[face_no / 2][std_face_ref];
3947 }
3948 
3949 
3950 
3951 template <int dim>
3952 inline RefinementCase<dim>
3953 GeometryInfo<dim>::min_cell_refinement_case_for_line_refinement(
3954  const unsigned int)
3955 {
3956  Assert(false, ExcNotImplemented());
3957 
3958  return RefinementCase<dim>::no_refinement;
3959 }
3960 
3961 template <>
3962 inline RefinementCase<1>
3963 GeometryInfo<1>::min_cell_refinement_case_for_line_refinement(
3964  const unsigned int line_no)
3965 {
3966  (void)line_no;
3967  AssertIndexRange(line_no, 1);
3968 
3969  return RefinementCase<1>::cut_x;
3970 }
3971 
3972 
3973 template <>
3974 inline RefinementCase<2>
3975 GeometryInfo<2>::min_cell_refinement_case_for_line_refinement(
3976  const unsigned int line_no)
3977 {
3978  const unsigned int dim = 2;
3979  (void)dim;
3980  AssertIndexRange(line_no, GeometryInfo<dim>::lines_per_cell);
3981 
3982  return (line_no / 2) != 0u ? RefinementCase<2>::cut_x :
3983  RefinementCase<2>::cut_y;
3984 }
3985 
3986 
3987 template <>
3988 inline RefinementCase<3>
3989 GeometryInfo<3>::min_cell_refinement_case_for_line_refinement(
3990  const unsigned int line_no)
3991 {
3992  const unsigned int dim = 3;
3993  AssertIndexRange(line_no, GeometryInfo<dim>::lines_per_cell);
3994 
3995  const RefinementCase<dim> ref_cases[6] = {
3996  RefinementCase<dim>::cut_y, // lines 0 and 1
3997  RefinementCase<dim>::cut_x, // lines 2 and 3
3998  RefinementCase<dim>::cut_y, // lines 4 and 5
3999  RefinementCase<dim>::cut_x, // lines 6 and 7
4000  RefinementCase<dim>::cut_z, // lines 8 and 9
4001  RefinementCase<dim>::cut_z}; // lines 10 and 11
4002 
4003  return ref_cases[line_no / 2];
4004 }
4005 
4006 
4007 
4008 template <>
4009 inline unsigned int
4010 GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
4011  const bool face_orientation,
4012  const bool face_flip,
4013  const bool face_rotation)
4014 {
4015  AssertIndexRange(vertex, GeometryInfo<3>::vertices_per_face);
4016 
4017  // set up a table to make sure that
4018  // we handle non-standard faces correctly
4019  //
4020  // so set up a table that for each vertex (of
4021  // a quad in standard position) describes
4022  // which vertex to take
4023  //
4024  // first index: four vertices 0...3
4025  //
4026  // second index: face_orientation; 0:
4027  // opposite normal, 1: standard
4028  //
4029  // third index: face_flip; 0: standard, 1:
4030  // face rotated by 180 degrees
4031  //
4032  // forth index: face_rotation: 0: standard,
4033  // 1: face rotated by 90 degrees
4034 
4035  const unsigned int vertex_translation[4][2][2][2] = {
4036  {{{0, 2}, // vertex 0, face_orientation=false, face_flip=false,
4037  // face_rotation=false and true
4038  {3, 1}}, // vertex 0, face_orientation=false, face_flip=true,
4039  // face_rotation=false and true
4040  {{0, 1}, // vertex 0, face_orientation=true, face_flip=false,
4041  // face_rotation=false and true
4042  {3, 2}}}, // vertex 0, face_orientation=true, face_flip=true,
4043  // face_rotation=false and true
4044 
4045  {{{2, 3}, // vertex 1 ...
4046  {1, 0}},
4047  {{1, 3}, {2, 0}}},
4048 
4049  {{{1, 0}, // vertex 2 ...
4050  {2, 3}},
4051  {{2, 0}, {1, 3}}},
4052 
4053  {{{3, 1}, // vertex 3 ...
4054  {0, 2}},
4055  {{3, 2}, {0, 1}}}};
4056 
4057  return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
4058 }
4059 
4060 
4061 
4062 template <int dim>
4063 inline unsigned int
4064 GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
4065  const bool,
4066  const bool,
4067  const bool)
4068 {
4069  Assert(dim > 1, ExcImpossibleInDim(dim));
4070  AssertIndexRange(vertex, GeometryInfo<dim>::vertices_per_face);
4071  return vertex;
4072 }
4073 
4074 
4075 
4076 template <>
4077 inline unsigned int
4078 GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
4079  const bool face_orientation,
4080  const bool face_flip,
4081  const bool face_rotation)
4082 {
4083  AssertIndexRange(line, GeometryInfo<3>::lines_per_face);
4084 
4085 
4086  // make sure we handle
4087  // non-standard faces correctly
4088  //
4089  // so set up a table that for each line (of a
4090  // quad) describes which line to take
4091  //
4092  // first index: four lines 0...3
4093  //
4094  // second index: face_orientation; 0:
4095  // opposite normal, 1: standard
4096  //
4097  // third index: face_flip; 0: standard, 1:
4098  // face rotated by 180 degrees
4099  //
4100  // forth index: face_rotation: 0: standard,
4101  // 1: face rotated by 90 degrees
4102 
4103  const unsigned int line_translation[4][2][2][2] = {
4104  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4105  // face_rotation=false and true
4106  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4107  // face_rotation=false and true
4108  {{0, 3}, // line 0, face_orientation=true, face_flip=false,
4109  // face_rotation=false and true
4110  {1, 2}}}, // line 0, face_orientation=true, face_flip=true,
4111  // face_rotation=false and true
4112 
4113  {{{3, 1}, // line 1 ...
4114  {2, 0}},
4115  {{1, 2}, {0, 3}}},
4116 
4117  {{{0, 3}, // line 2 ...
4118  {1, 2}},
4119  {{2, 0}, {3, 1}}},
4120 
4121  {{{1, 2}, // line 3 ...
4122  {0, 3}},
4123  {{3, 1}, {2, 0}}}};
4124 
4125  return line_translation[line][face_orientation][face_flip][face_rotation];
4126 }
4127 
4128 
4129 
4130 template <int dim>
4131 inline unsigned int
4132 GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
4133  const bool,
4134  const bool,
4135  const bool)
4136 {
4137  Assert(false, ExcNotImplemented());
4138  return line;
4139 }
4140 
4141 
4142 
4143 template <>
4144 inline unsigned int
4145 GeometryInfo<2>::standard_to_real_line_vertex(const unsigned int vertex,
4146  const bool line_orientation)
4147 {
4148  return line_orientation ? vertex : (1 - vertex);
4149 }
4150 
4151 
4152 
4153 template <int dim>
4154 inline unsigned int
4155 GeometryInfo<dim>::standard_to_real_line_vertex(const unsigned int vertex,
4156  const bool)
4157 {
4158  Assert(false, ExcNotImplemented());
4159  return vertex;
4160 }
4161 
4162 
4163 
4164 template <>
4165 inline std::array<unsigned int, 2>
4166 GeometryInfo<2>::standard_quad_vertex_to_line_vertex_index(
4167  const unsigned int vertex)
4168 {
4169  return {{vertex % 2, vertex / 2}};
4170 }
4171 
4172 
4173 
4174 template <int dim>
4175 inline std::array<unsigned int, 2>
4176 GeometryInfo<dim>::standard_quad_vertex_to_line_vertex_index(
4177  const unsigned int vertex)
4178 {
4179  Assert(false, ExcNotImplemented());
4180  (void)vertex;
4181  return {{0, 0}};
4182 }
4183 
4184 
4185 
4186 template <>
4187 inline std::array<unsigned int, 2>
4188 GeometryInfo<3>::standard_hex_line_to_quad_line_index(const unsigned int i)
4189 {
4190  // set up a table that for each
4191  // line describes a) from which
4192  // quad to take it, b) which line
4193  // therein it is if the face is
4194  // oriented correctly
4195  static const unsigned int lookup_table[GeometryInfo<3>::lines_per_cell][2] = {
4196  {4, 0}, // take first four lines from bottom face
4197  {4, 1},
4198  {4, 2},
4199  {4, 3},
4200 
4201  {5, 0}, // second four lines from top face
4202  {5, 1},
4203  {5, 2},
4204  {5, 3},
4205 
4206  {0, 0}, // the rest randomly
4207  {1, 0},
4208  {0, 1},
4209  {1, 1}};
4210 
4211  return {{lookup_table[i][0], lookup_table[i][1]}};
4212 }
4213 
4214 
4215 
4216 template <int dim>
4217 inline std::array<unsigned int, 2>
4218 GeometryInfo<dim>::standard_hex_line_to_quad_line_index(const unsigned int line)
4219 {
4220  Assert(false, ExcNotImplemented());
4221  (void)line;
4222  return {{0, 0}};
4223 }
4224 
4225 
4226 
4227 template <>
4228 inline std::array<unsigned int, 2>
4229 GeometryInfo<3>::standard_hex_vertex_to_quad_vertex_index(
4230  const unsigned int vertex)
4231 {
4232  // get the corner indices by asking either the bottom or the top face for its
4233  // vertices. handle non-standard faces by calling the vertex reordering
4234  // function from GeometryInfo
4235 
4236  // bottom face (4) for first four vertices, top face (5) for the rest
4237  return {{4 + vertex / 4, vertex % 4}};
4238 }
4239 
4240 
4241 
4242 template <int dim>
4243 inline std::array<unsigned int, 2>
4244 GeometryInfo<dim>::standard_hex_vertex_to_quad_vertex_index(
4245  const unsigned int vertex)
4246 {
4247  Assert(false, ExcNotImplemented());
4248  (void)vertex;
4249  return {{0, 0}};
4250 }
4251 
4252 
4253 
4254 template <>
4255 inline unsigned int
4256 GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
4257  const bool face_orientation,
4258  const bool face_flip,
4259  const bool face_rotation)
4260 {
4261  AssertIndexRange(line, GeometryInfo<3>::lines_per_face);
4262 
4263 
4264  // make sure we handle
4265  // non-standard faces correctly
4266  //
4267  // so set up a table that for each line (of a
4268  // quad) describes which line to take
4269  //
4270  // first index: four lines 0...3
4271  //
4272  // second index: face_orientation; 0:
4273  // opposite normal, 1: standard
4274  //
4275  // third index: face_flip; 0: standard, 1:
4276  // face rotated by 180 degrees
4277  //
4278  // forth index: face_rotation: 0: standard,
4279  // 1: face rotated by 90 degrees
4280 
4281  const unsigned int line_translation[4][2][2][2] = {
4282  {{{2, 0}, // line 0, face_orientation=false, face_flip=false,
4283  // face_rotation=false and true
4284  {3, 1}}, // line 0, face_orientation=false, face_flip=true,
4285  // face_rotation=false and true
4286  {{0, 2}, // line 0, face_orientation=true, face_flip=false,
4287  // face_rotation=false and true
4288  {1, 3}}}, // line 0, face_orientation=true, face_flip=true,
4289  // face_rotation=false and true
4290 
4291  {{{3, 1}, // line 1 ...
4292  {2, 0}},
4293  {{1, 3}, {0, 2}}},
4294 
4295  {{{0, 3}, // line 2 ...
4296  {1, 2}},
4297  {{2, 1}, {3, 0}}},
4298 
4299  {{{1, 2}, // line 3 ...
4300  {0, 3}},
4301  {{3, 0}, {2, 1}}}};
4302 
4303  return line_translation[line][face_orientation][face_flip][face_rotation];
4304 }
4305 
4306 
4307 
4308 template <int dim>
4309 inline unsigned int
4310 GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
4311  const bool,
4312  const bool,
4313  const bool)
4314 {
4315  Assert(false, ExcNotImplemented());
4316  return line;
4317 }
4318 
4319 
4320 
4321 template <>
4322 inline unsigned int
4323 GeometryInfo<1>::child_cell_on_face(const RefinementCase<1> &,
4324  const unsigned int face,
4325  const unsigned int subface,
4326  const bool,
4327  const bool,
4328  const bool,
4329  const RefinementCase<0> &)
4330 {
4331  (void)subface;
4332  AssertIndexRange(face, faces_per_cell);
4333  AssertIndexRange(subface, max_children_per_face);
4334 
4335  return face;
4336 }
4337 
4338 
4339 
4340 template <>
4341 inline unsigned int
4342 GeometryInfo<2>::child_cell_on_face(const RefinementCase<2> &ref_case,
4343  const unsigned int face,
4344  const unsigned int subface,
4345  const bool /*face_orientation*/,
4346  const bool face_flip,
4347  const bool /*face_rotation*/,
4348  const RefinementCase<1> &)
4349 {
4350  AssertIndexRange(face, faces_per_cell);
4351  AssertIndexRange(subface, max_children_per_face);
4352 
4353  // always return the child adjacent to the specified
4354  // subface. if the face of a cell is not refined, don't
4355  // throw an assertion but deliver the child adjacent to
4356  // the face nevertheless, i.e. deliver the child of
4357  // this cell adjacent to the subface of a possibly
4358  // refined neighbor. this simplifies setting neighbor
4359  // information in execute_refinement.
4360  const unsigned int
4361  subcells[2][RefinementCase<2>::isotropic_refinement][faces_per_cell]
4362  [max_children_per_face] = {
4363  {
4364  // Normal orientation (face_flip = false)
4365  {{0, 0}, {1, 1}, {0, 1}, {0, 1}}, // cut_x
4366  {{0, 1}, {0, 1}, {0, 0}, {1, 1}}, // cut_y
4367  {{0, 2}, {1, 3}, {0, 1}, {2, 3}} // cut_xy, i.e., isotropic
4368  },
4369  {
4370  // Flipped orientation (face_flip = true)
4371  {{0, 0}, {1, 1}, {1, 0}, {1, 0}}, // cut_x
4372  {{1, 0}, {1, 0}, {0, 0}, {1, 1}}, // cut_y
4373  {{2, 0}, {3, 1}, {1, 0}, {3, 2}} // cut_xy, i.e., isotropic
4374  }};
4375 
4376  return subcells[face_flip][ref_case - 1][face][subface];
4377 }
4378 
4379 
4380 
4381 template <>
4382 inline unsigned int
4383 GeometryInfo<3>::child_cell_on_face(const RefinementCase<3> &ref_case,
4384  const unsigned int face,
4385  const unsigned int subface,
4386  const bool face_orientation,
4387  const bool face_flip,
4388  const bool face_rotation,
4389  const RefinementCase<2> &face_ref_case)
4390 {
4391  const unsigned int dim = 3;
4392 
4393  Assert(ref_case > RefinementCase<dim - 1>::no_refinement,
4394  ExcMessage("Cell has no children."));
4395  AssertIndexRange(face, faces_per_cell);
4396  if (!(subface == 0 &&
4397  face_ref_case == RefinementCase<dim - 1>::no_refinement))
4398  {
4399  AssertIndexRange(subface,
4400  GeometryInfo<dim - 1>::n_children(face_ref_case));
4401  }
4402 
4403  // invalid number used for invalid cases,
4404  // e.g. when the children are more refined at
4405  // a given face than the face itself
4406  const unsigned int e = numbers::invalid_unsigned_int;
4407 
4408  // the whole process of finding a child cell
4409  // at a given subface considering the
4410  // possibly anisotropic refinement cases of
4411  // the cell and the face as well as
4412  // orientation, flip and rotation of the face
4413  // is quite complicated. thus, we break it
4414  // down into several steps.
4415 
4416  // first step: convert the given face refine
4417  // case to a face refine case concerning the
4418  // face in standard orientation (, flip and
4419  // rotation). This only affects cut_x and
4420  // cut_y
4421  const RefinementCase<dim - 1> flip[4] = {
4422  RefinementCase<dim - 1>::no_refinement,
4423  RefinementCase<dim - 1>::cut_y,
4424  RefinementCase<dim - 1>::cut_x,
4425  RefinementCase<dim - 1>::cut_xy};
4426  // for face_orientation, 'true' is the
4427  // default value whereas for face_rotation,
4428  // 'false' is standard. If
4429  // <tt>face_rotation==face_orientation</tt>,
4430  // then one of them is non-standard and we
4431  // have to swap cut_x and cut_y, otherwise no
4432  // change is necessary.
4433  const RefinementCase<dim - 1> std_face_ref =
4434  (face_orientation == face_rotation) ? flip[face_ref_case] : face_ref_case;
4435 
4436  // second step: convert the given subface
4437  // index to the one for a standard face
4438  // respecting face_orientation, face_flip and
4439  // face_rotation
4440 
4441  // first index: face_ref_case
4442  // second index: face_orientation
4443  // third index: face_flip
4444  // forth index: face_rotation
4445  // fifth index: subface index
4446  const unsigned int subface_exchange[4][2][2][2][4] = {
4447  // no_refinement (subface 0 stays 0,
4448  // all others are invalid)
4449  {{{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}},
4450  {{{0, e, e, e}, {0, e, e, e}}, {{0, e, e, e}, {0, e, e, e}}}},
4451  // cut_x (here, if the face is only
4452  // rotated OR only falsely oriented,
4453  // then subface 0 of the non-standard
4454  // face does NOT correspond to one of
4455  // the subfaces of a standard
4456  // face. Thus we indicate the subface
4457  // which is located at the lower left
4458  // corner (the origin of the face's
4459  // local coordinate system) with
4460  // '0'. The rest of this issue is
4461  // taken care of using the above
4462  // conversion to a 'standard face
4463  // refine case')
4464  {{{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}},
4465  {{{0, 1, e, e}, {0, 1, e, e}}, {{1, 0, e, e}, {1, 0, e, e}}}},
4466  // cut_y (the same applies as for
4467  // cut_x)
4468  {{{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}},
4469  {{{0, 1, e, e}, {1, 0, e, e}}, {{1, 0, e, e}, {0, 1, e, e}}}},
4470  // cut_xyz: this information is
4471  // identical to the information
4472  // returned by
4473  // GeometryInfo<3>::real_to_standard_face_vertex()
4474  {{{{0, 2, 1, 3}, // face_orientation=false, face_flip=false,
4475  // face_rotation=false, subfaces 0,1,2,3
4476  {2, 3, 0, 1}}, // face_orientation=false, face_flip=false,
4477  // face_rotation=true, subfaces 0,1,2,3
4478  {{3, 1, 2, 0}, // face_orientation=false, face_flip=true,
4479  // face_rotation=false, subfaces 0,1,2,3
4480  {1, 0, 3, 2}}}, // face_orientation=false, face_flip=true,
4481  // face_rotation=true, subfaces 0,1,2,3
4482  {{{0, 1, 2, 3}, // face_orientation=true, face_flip=false,
4483  // face_rotation=false, subfaces 0,1,2,3
4484  {1, 3, 0, 2}}, // face_orientation=true, face_flip=false,
4485  // face_rotation=true, subfaces 0,1,2,3
4486  {{3, 2, 1, 0}, // face_orientation=true, face_flip=true,
4487  // face_rotation=false, subfaces 0,1,2,3
4488  {2, 0, 3, 1}}}}}; // face_orientation=true, face_flip=true,
4489  // face_rotation=true, subfaces 0,1,2,3
4490 
4491  const unsigned int std_subface =
4492  subface_exchange[face_ref_case][face_orientation][face_flip][face_rotation]
4493  [subface];
4494  Assert(std_subface != e, ExcInternalError());
4495 
4496  // third step: these are the children, which
4497  // can be found at the given subfaces of an
4498  // isotropically refined (standard) face
4499  //
4500  // first index: (refinement_case-1)
4501  // second index: face_index
4502  // third index: subface_index (isotropic refinement)
4503  const unsigned int iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell]
4504  [max_children_per_face] = {
4505  // cut_x
4506  {{0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
4507  {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
4508  {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
4509  {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
4510  {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
4511  {0, 1, 0, 1}}, // face 5, subfaces 0,1,2,3
4512  // cut_y
4513  {{0, 1, 0, 1},
4514  {0, 1, 0, 1},
4515  {0, 0, 0, 0},
4516  {1, 1, 1, 1},
4517  {0, 0, 1, 1},
4518  {0, 0, 1, 1}},
4519  // cut_xy
4520  {{0, 2, 0, 2},
4521  {1, 3, 1, 3},
4522  {0, 0, 1, 1},
4523  {2, 2, 3, 3},
4524  {0, 1, 2, 3},
4525  {0, 1, 2, 3}},
4526  // cut_z
4527  {{0, 0, 1, 1},
4528  {0, 0, 1, 1},
4529  {0, 1, 0, 1},
4530  {0, 1, 0, 1},
4531  {0, 0, 0, 0},
4532  {1, 1, 1, 1}},
4533  // cut_xz
4534  {{0, 0, 1, 1},
4535  {2, 2, 3, 3},
4536  {0, 1, 2, 3},
4537  {0, 1, 2, 3},
4538  {0, 2, 0, 2},
4539  {1, 3, 1, 3}},
4540  // cut_yz
4541  {{0, 1, 2, 3},
4542  {0, 1, 2, 3},
4543  {0, 2, 0, 2},
4544  {1, 3, 1, 3},
4545  {0, 0, 1, 1},
4546  {2, 2, 3, 3}},
4547  // cut_xyz
4548  {{0, 2, 4, 6},
4549  {1, 3, 5, 7},
4550  {0, 4, 1, 5},
4551  {2, 6, 3, 7},
4552  {0, 1, 2, 3},
4553  {4, 5, 6, 7}}};
4554 
4555  // forth step: check, whether the given face
4556  // refine case is valid for the given cell
4557  // refine case. this is the case, if the
4558  // given face refine case is at least as
4559  // refined as the face is for the given cell
4560  // refine case
4561 
4562  // note, that we are considering standard
4563  // face refinement cases here and thus must
4564  // not pass the given orientation, flip and
4565  // rotation flags
4566  if ((std_face_ref & face_refinement_case(ref_case, face)) ==
4567  face_refinement_case(ref_case, face))
4568  {
4569  // all is fine. for anisotropic face
4570  // refine cases, select one of the
4571  // isotropic subfaces which neighbors the
4572  // same child
4573 
4574  // first index: (standard) face refine case
4575  // second index: subface index
4576  const unsigned int equivalent_iso_subface[4][4] = {
4577  {0, e, e, e}, // no_refinement
4578  {0, 3, e, e}, // cut_x
4579  {0, 3, e, e}, // cut_y
4580  {0, 1, 2, 3}}; // cut_xy
4581 
4582  const unsigned int equ_std_subface =
4583  equivalent_iso_subface[std_face_ref][std_subface];
4584  Assert(equ_std_subface != e, ExcInternalError());
4585 
4586  return iso_children[ref_case - 1][face][equ_std_subface];
4587  }
4588  else
4589  {
4590  // the face_ref_case was too coarse,
4591  // throw an error
4592  Assert(false,
4593  ExcMessage("The face RefineCase is too coarse "
4594  "for the given cell RefineCase."));
4595  }
4596  // we only get here in case of an error
4597  return e;
4598 }
4599 
4600 
4601 
4602 template <>
4603 inline unsigned int
4604 GeometryInfo<4>::child_cell_on_face(const RefinementCase<4> &,
4605  const unsigned int,
4606  const unsigned int,
4607  const bool,
4608  const bool,
4609  const bool,
4610  const RefinementCase<3> &)
4611 {
4612  Assert(false, ExcNotImplemented());
4613  return numbers::invalid_unsigned_int;
4614 }
4615 
4616 
4617 
4618 template <>
4619 inline unsigned int
4620 GeometryInfo<2>::face_to_cell_lines(const unsigned int face,
4621  const unsigned int line,
4622  const bool,
4623  const bool,
4624  const bool)
4625 {
4626  (void)line;
4627  AssertIndexRange(face, faces_per_cell);
4628  AssertIndexRange(line, lines_per_face);
4629 
4630  // The face is a line itself.
4631  return face;
4632 }
4633 
4634 
4635 
4636 template <>
4637 inline unsigned int
4638 GeometryInfo<3>::face_to_cell_lines(const unsigned int face,
4639  const unsigned int line,
4640  const bool face_orientation,
4641  const bool face_flip,
4642  const bool face_rotation)
4643 {
4644  AssertIndexRange(face, faces_per_cell);
4645  AssertIndexRange(line, lines_per_face);
4646 
4647  const unsigned lines[faces_per_cell][lines_per_face] = {
4648  {8, 10, 0, 4}, // left face
4649  {9, 11, 1, 5}, // right face
4650  {2, 6, 8, 9}, // front face
4651  {3, 7, 10, 11}, // back face
4652  {0, 1, 2, 3}, // bottom face
4653  {4, 5, 6, 7}}; // top face
4654  return lines[face][real_to_standard_face_line(
4655  line, face_orientation, face_flip, face_rotation)];
4656 }
4657 
4658 
4659 
4660 inline unsigned int
4661 GeometryInfo<0>::face_to_cell_lines(const unsigned int,
4662  const unsigned int,
4663  const bool,
4664  const bool,
4665  const bool)
4666 {
4667  Assert(false, ExcNotImplemented());
4668  return numbers::invalid_unsigned_int;
4669 }
4670 
4671 
4672 
4673 template <int dim>
4674 inline unsigned int
4675 GeometryInfo<dim>::face_to_cell_lines(const unsigned int,
4676  const unsigned int,
4677  const bool,
4678  const bool,
4679  const bool)
4680 {
4681  Assert(false, ExcNotImplemented());
4682  return numbers::invalid_unsigned_int;
4683 }
4684 
4685 
4686 
4687 template <int dim>
4688 inline unsigned int
4689 GeometryInfo<dim>::face_to_cell_vertices(const unsigned int face,
4690  const unsigned int vertex,
4691  const bool face_orientation,
4692  const bool face_flip,
4693  const bool face_rotation)
4694 {
4695  return child_cell_on_face(RefinementCase<dim>::isotropic_refinement,
4696  face,
4697  vertex,
4698  face_orientation,
4699  face_flip,
4700  face_rotation);
4701 }
4702 
4703 
4704 
4705 inline unsigned int
4706 GeometryInfo<0>::face_to_cell_vertices(const unsigned int,
4707  const unsigned int,
4708  const bool,
4709  const bool,
4710  const bool)
4711 {
4712  Assert(false, ExcNotImplemented());
4713  return numbers::invalid_unsigned_int;
4714 }
4715 
4716 
4717 
4718 template <int dim>
4719 template <typename Number>
4720 inline Point<dim, Number>
4721 GeometryInfo<dim>::project_to_unit_cell(const Point<dim, Number> &q)
4722 {
4723  Point<dim, Number> p;
4724  for (unsigned int i = 0; i < dim; ++i)
4725  p[i] = std::min(std::max(q[i], Number(0.)), Number(1.));
4726 
4727  return p;
4728 }
4729 
4730 
4731 
4732 template <int dim>
4733 inline double
4734 GeometryInfo<dim>::distance_to_unit_cell(const Point<dim> &p)
4735 {
4736  double result = 0.0;
4737 
4738  for (unsigned int i = 0; i < dim; ++i)
4739  {
4740  result = std::max(result, -p[i]);
4741  result = std::max(result, p[i] - 1.);
4742  }
4743 
4744  return result;
4745 }
4746 
4747 
4748 
4749 template <int dim>
4750 inline double
4751 GeometryInfo<dim>::d_linear_shape_function(const Point<dim> & xi,
4752  const unsigned int i)
4753 {
4754  AssertIndexRange(i, GeometryInfo<dim>::vertices_per_cell);
4755 
4756  switch (dim)
4757  {
4758  case 1:
4759  {
4760  const double x = xi[0];
4761  switch (i)
4762  {
4763  case 0:
4764  return 1 - x;
4765  case 1:
4766  return x;
4767  }
4768  break;
4769  }
4770 
4771  case 2:
4772  {
4773  const double x = xi[0];
4774  const double y = xi[1];
4775  switch (i)
4776  {
4777  case 0:
4778  return (1 - x) * (1 - y);
4779  case 1:
4780  return x * (1 - y);
4781  case 2:
4782  return (1 - x) * y;
4783  case 3:
4784  return x * y;
4785  }
4786  break;
4787  }
4788 
4789  case 3:
4790  {
4791  const double x = xi[0];
4792  const double y = xi[1];
4793  const double z = xi[2];
4794  switch (i)
4795  {
4796  case 0:
4797  return (1 - x) * (1 - y) * (1 - z);
4798  case 1:
4799  return x * (1 - y) * (1 - z);
4800  case 2:
4801  return (1 - x) * y * (1 - z);
4802  case 3:
4803  return x * y * (1 - z);
4804  case 4:
4805  return (1 - x) * (1 - y) * z;
4806  case 5:
4807  return x * (1 - y) * z;
4808  case 6:
4809  return (1 - x) * y * z;
4810  case 7:
4811  return x * y * z;
4812  }
4813  break;
4814  }
4815 
4816  default:
4817  Assert(false, ExcNotImplemented());
4818  }
4819  return -1e9;
4820 }
4821 
4822 
4823 
4824 template <>
4825 Tensor<1, 1> inline GeometryInfo<1>::d_linear_shape_function_gradient(
4826  const Point<1> &,
4827  const unsigned int i)
4828 {
4829  AssertIndexRange(i, GeometryInfo<1>::vertices_per_cell);
4830 
4831  switch (i)
4832  {
4833  case 0:
4834  return Point<1>(-1.);
4835  case 1:
4836  return Point<1>(1.);
4837  }
4838 
4839  return Point<1>(-1e9);
4840 }
4841 
4842 
4843 
4844 template <>
4845 Tensor<1, 2> inline GeometryInfo<2>::d_linear_shape_function_gradient(
4846  const Point<2> & xi,
4847  const unsigned int i)
4848 {
4849  AssertIndexRange(i, GeometryInfo<2>::vertices_per_cell);
4850 
4851  const double x = xi[0];
4852  const double y = xi[1];
4853  switch (i)
4854  {
4855  case 0:
4856  return Point<2>(-(1 - y), -(1 - x));
4857  case 1:
4858  return Point<2>(1 - y, -x);
4859  case 2:
4860  return Point<2>(-y, 1 - x);
4861  case 3:
4862  return Point<2>(y, x);
4863  }
4864  return Point<2>(-1e9, -1e9);
4865 }
4866 
4867 
4868 
4869 template <>
4870 Tensor<1, 3> inline GeometryInfo<3>::d_linear_shape_function_gradient(
4871  const Point<3> & xi,
4872  const unsigned int i)
4873 {
4874  AssertIndexRange(i, GeometryInfo<3>::vertices_per_cell);
4875 
4876  const double x = xi[0];
4877  const double y = xi[1];
4878  const double z = xi[2];
4879  switch (i)
4880  {
4881  case 0:
4882  return Point<3>(-(1 - y) * (1 - z),
4883  -(1 - x) * (1 - z),
4884  -(1 - x) * (1 - y));
4885  case 1:
4886  return Point<3>((1 - y) * (1 - z), -x * (1 - z), -x * (1 - y));
4887  case 2:
4888  return Point<3>(-y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y);
4889  case 3:
4890  return Point<3>(y * (1 - z), x * (1 - z), -x * y);
4891  case 4:
4892  return Point<3>(-(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y));
4893  case 5:
4894  return Point<3>((1 - y) * z, -x * z, x * (1 - y));
4895  case 6:
4896  return Point<3>(-y * z, (1 - x) * z, (1 - x) * y);
4897  case 7:
4898  return Point<3>(y * z, x * z, x * y);
4899  }
4900 
4901  return Point<3>(-1e9, -1e9, -1e9);
4902 }
4903 
4904 
4905 
4906 template <int dim>
4907 inline Tensor<1, dim>
4908 GeometryInfo<dim>::d_linear_shape_function_gradient(const Point<dim> &,
4909  const unsigned int)
4910 {
4911  Assert(false, ExcNotImplemented());
4912  return Tensor<1, dim>();
4913 }
4914 
4915 
4916 unsigned int inline GeometryInfo<0>::n_children(const RefinementCase<0> &)
4917 {
4918  return 0;
4919 }
4920 
4921 
4922 namespace internal
4923 {
4924  namespace GeometryInfoHelper
4925  {
4926  // wedge product of a single
4927  // vector in 2d: we just have to
4928  // rotate it by 90 degrees to the
4929  // right
4930  inline Tensor<1, 2>
4931  wedge_product(const Tensor<1, 2> (&derivative)[1])
4932  {
4933  Tensor<1, 2> result;
4934  result[0] = derivative[0][1];
4935  result[1] = -derivative[0][0];
4936 
4937  return result;
4938  }
4939 
4940 
4941  // wedge product of 2 vectors in
4942  // 3d is the cross product
4943  inline Tensor<1, 3>
4944  wedge_product(const Tensor<1, 3> (&derivative)[2])
4945  {
4946  return cross_product_3d(derivative[0], derivative[1]);
4947  }
4948 
4949 
4950  // wedge product of dim vectors
4951  // in dim-d: that's the
4952  // determinant of the matrix
4953  template <int dim>
4954  inline Tensor<0, dim>
4955  wedge_product(const Tensor<1, dim> (&derivative)[dim])
4956  {
4957  Tensor<2, dim> jacobian;
4958  for (unsigned int i = 0; i < dim; ++i)
4959  jacobian[i] = derivative[i];
4960 
4961  return determinant(jacobian);
4962  }
4963  } // namespace GeometryInfoHelper
4964 } // namespace internal
4965 
4966 
4967 template <int dim>
4968 template <int spacedim>
4969 inline void
4970 GeometryInfo<dim>::alternating_form_at_vertices
4971 # ifndef DEAL_II_CXX14_CONSTEXPR_BUG
4972  (const Point<spacedim> (&vertices)[vertices_per_cell],
4973  Tensor<spacedim - dim, spacedim> (&forms)[vertices_per_cell])
4974 # else
4975  (const Point<spacedim> *vertices, Tensor<spacedim - dim, spacedim> *forms)
4976 # endif
4977 {
4978  // for each of the vertices,
4979  // compute the alternating form
4980  // of the mapped unit
4981  // vectors. consider for
4982  // example the case of a quad
4983  // in spacedim==3: to do so, we
4984  // need to see how the
4985  // infinitesimal vectors
4986  // (d\xi_1,0) and (0,d\xi_2)
4987  // are transformed into
4988  // spacedim-dimensional space
4989  // and then form their cross
4990  // product (i.e. the wedge product
4991  // of two vectors). to this end, note
4992  // that
4993  // \vec x = sum_i \vec v_i phi_i(\vec xi)
4994  // so the transformed vectors are
4995  // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
4996  // and
4997  // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
4998  // which boils down to the columns
4999  // of the 3x2 matrix \grad_\xi x(\xi)
5000  //
5001  // a similar reasoning would
5002  // hold for all dim,spacedim
5003  // pairs -- we only have to
5004  // compute the wedge product of
5005  // the columns of the
5006  // derivatives
5007  for (unsigned int i = 0; i < vertices_per_cell; ++i)
5008  {
5009  Tensor<1, spacedim> derivatives[dim];
5010 
5011  for (unsigned int j = 0; j < vertices_per_cell; ++j)
5012  {
5013  const Tensor<1, dim> grad_phi_j =
5014  d_linear_shape_function_gradient(unit_cell_vertex(i), j);
5015  for (unsigned int l = 0; l < dim; ++l)
5016  derivatives[l] += vertices[j] * grad_phi_j[l];
5017  }
5018 
5019  forms[i] = internal::GeometryInfoHelper::wedge_product(derivatives);
5020  }
5021 }
5022 
5023 #endif // DOXYGEN
5024 DEAL_II_NAMESPACE_CLOSE
5025 
5026 #endif
GeometryPrimitive::GeometryPrimitive
GeometryPrimitive(const unsigned int object_dimension)
GeometryPrimitive::GeometryPrimitive
GeometryPrimitive(const Object object)
Point
Definition: point.h:111
RefinementCase::operator|
RefinementCase operator|(const RefinementCase &r) const
RefinementCase::memory_consumption
static std::size_t memory_consumption()
RefinementCase::operator~
RefinementCase operator~() const
RefinementCase::operator&
RefinementCase operator&(const RefinementCase &r) const
RefinementCase::RefinementCase
RefinementCase(const typename RefinementPossibilities< dim >::Possibilities refinement_case)
RefinementCase::cut_axis
static RefinementCase cut_axis(const unsigned int i)
RefinementCase::value
std::uint8_t value
Definition: geometry_info.h:887
RefinementCase::serialize
void serialize(Archive &ar, const unsigned int version)
RefinementCase::RefinementCase
RefinementCase(const std::uint8_t refinement_case)
Tensor
Definition: tensor.h:503
internal::SubfaceCase::memory_consumption
static constexpr std::size_t memory_consumption()
internal::SubfaceCase::SubfaceCase
SubfaceCase(const typename SubfacePossibilities< dim >::Possibilities subface_possibility)
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
vertices
Point< 3 > vertices[4]
Definition: data_out_base.cc:192
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
GeometryInfo::ExcInvalidCoordinate
static ::ExceptionBase & ExcInvalidCoordinate(double arg1)
GeometryInfo::ExcInvalidSubface
static ::ExceptionBase & ExcInvalidSubface(int arg1, int arg2, int arg3)
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1473
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
internal::SubfaceCase::ExcInvalidSubfaceCase
static ::ExceptionBase & ExcInvalidSubfaceCase(int arg1)
StandardExceptions::ExcImpossibleInDim
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
DataOutBase::eps
@ eps
Definition: data_out_base.h:1621
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
DeclException3
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:555
DeclException1
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
RefinementCase::ExcInvalidRefinementCase
static ::ExceptionBase & ExcInvalidRefinementCase(int arg1)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
iota_view.h
LAPACKSupport::U
static const char U
Definition: lapack_support.h:169
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:190
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:201
std_cxx20::ranges::iota_view
boost::integer_range< IncrementableType > iota_view
Definition: iota_view.h:46
ndarray
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108
GeometryInfo< 0 >::n_children
static unsigned int n_children(const RefinementCase< 0 > &refinement_case)
GeometryInfo< 0 >::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< 0 >::dx_to_deal
static const std::array< unsigned int, vertices_per_cell > dx_to_deal
Definition: geometry_info.h:1432
GeometryInfo< 0 >::vertex_indices
static std::array< unsigned int, vertices_per_cell > vertex_indices()
GeometryInfo< 0 >::face_indices
static std::array< unsigned int, 0 > face_indices()
GeometryInfo< 0 >::face_to_cell_lines
static unsigned int face_to_cell_lines(const unsigned int face, const unsigned int line, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
GeometryInfo< 0 >::ucd_to_deal
static const std::array< unsigned int, vertices_per_cell > ucd_to_deal
Definition: geometry_info.h:1417
GeometryInfo::child_cell_from_point
static unsigned int child_cell_from_point(const Point< dim > &p)
GeometryInfo::standard_to_real_line_vertex
static unsigned int standard_to_real_line_vertex(const unsigned int vertex, const bool line_orientation=true)
GeometryInfo::ucd_to_deal
static constexpr std::array< unsigned int, vertices_per_cell > ucd_to_deal
Definition: geometry_info.h:2102
GeometryInfo::unit_tangential_vectors
static constexpr ndarray< Tensor< 1, dim >, faces_per_cell, dim - 1 > unit_tangential_vectors
Definition: geometry_info.h:2634
GeometryInfo::is_inside_unit_cell
static bool is_inside_unit_cell(const Point< dim > &p, const double eps)
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)
GeometryInfo::distance_to_unit_cell
static double distance_to_unit_cell(const Point< dim > &p)
GeometryInfo::child_cell_on_face
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)
GeometryInfo::quads_per_face
static constexpr unsigned int quads_per_face
Definition: geometry_info.h:2051
GeometryInfo::standard_to_real_face_vertex
static unsigned int standard_to_real_face_vertex(const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
GeometryInfo::face_indices
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
GeometryInfo::opposite_face
static constexpr std::array< unsigned int, faces_per_cell > opposite_face
Definition: geometry_info.h:2642
GeometryInfo::max_children_per_face
static constexpr unsigned int max_children_per_face
Definition: geometry_info.h:2010
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::vertices_per_cell
static constexpr unsigned int vertices_per_cell
Definition: geometry_info.h:2016
GeometryInfo::standard_hex_line_to_quad_line_index
static std::array< unsigned int, 2 > standard_hex_line_to_quad_line_index(const unsigned int line)
GeometryInfo::unit_cell_vertex
static Point< dim > unit_cell_vertex(const unsigned int vertex)
GeometryInfo::dx_to_deal
static constexpr std::array< unsigned int, vertices_per_cell > dx_to_deal
Definition: geometry_info.h:2118
GeometryInfo::vertex_indices
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
GeometryInfo::lines_per_cell
static constexpr unsigned int lines_per_cell
Definition: geometry_info.h:2063
GeometryInfo::line_refinement_case
static RefinementCase< 1 > line_refinement_case(const RefinementCase< dim > &cell_refinement_case, const unsigned int line_no)
GeometryInfo::unit_normal_direction
static constexpr std::array< unsigned int, faces_per_cell > unit_normal_direction
Definition: geometry_info.h:2584
GeometryInfo::min_cell_refinement_case_for_line_refinement
static RefinementCase< dim > min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
GeometryInfo::faces_per_cell
static constexpr unsigned int faces_per_cell
Definition: geometry_info.h:1982
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)
GeometryInfo::subface_ratio
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
GeometryInfo::hexes_per_cell
static constexpr unsigned int hexes_per_cell
Definition: geometry_info.h:2081
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::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::d_linear_shape_function
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
GeometryInfo::unit_normal_vector
static constexpr std::array< Tensor< 1, dim >, faces_per_cell > unit_normal_vector
Definition: geometry_info.h:2617
GeometryInfo::project_to_unit_cell
static Point< dim, Number > project_to_unit_cell(const Point< dim, Number > &p)
GeometryInfo::d_linear_shape_function_gradient
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
GeometryInfo::vertices_per_face
static constexpr unsigned int vertices_per_face
Definition: geometry_info.h:2039
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::line_to_cell_vertices
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
GeometryInfo::is_inside_unit_cell
static bool is_inside_unit_cell(const Point< dim > &p)
GeometryInfo::n_children
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
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::face_to_cell_lines
static unsigned int face_to_cell_lines(const unsigned int face, const unsigned int line, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
GeometryInfo::standard_quad_vertex_to_line_vertex_index
static std::array< unsigned int, 2 > standard_quad_vertex_to_line_vertex_index(const unsigned int vertex)
GeometryInfo::vertex_to_face
static constexpr ndarray< unsigned int, vertices_per_cell, dim > vertex_to_face
Definition: geometry_info.h:2132
GeometryInfo::alternating_form_at_vertices
static void alternating_form_at_vertices(const Point< spacedim >(&vertices)[vertices_per_cell], Tensor< spacedim - dim, spacedim >(&forms)[vertices_per_cell])
GeometryInfo::face_refinement_case
static RefinementCase< dim - 1 > face_refinement_case(const RefinementCase< dim > &cell_refinement_case, const unsigned int face_no, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
GeometryInfo::unit_normal_orientation
static constexpr std::array< int, faces_per_cell > unit_normal_orientation
Definition: geometry_info.h:2603
GeometryInfo::quads_per_cell
static constexpr unsigned int quads_per_cell
Definition: geometry_info.h:2074
GeometryInfo::max_children_per_cell
static constexpr unsigned int max_children_per_cell
Definition: geometry_info.h:1977
GeometryInfo::n_subfaces
static unsigned int n_subfaces(const internal::SubfaceCase< dim > &subface_case)
GeometryInfo::lines_per_face
static constexpr unsigned int lines_per_face
Definition: geometry_info.h:2045
GeometryInfo::standard_hex_vertex_to_quad_vertex_index
static std::array< unsigned int, 2 > standard_hex_vertex_to_quad_vertex_index(const unsigned int vertex)
determinant
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:2745
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:2660
Generated by doxygen 1.9.1