Reference documentation for deal.II version Git 191d06ed00 2021-05-11 21:15:49 -0400
\(\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\}}\)
reference_cell.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2020 - 2021 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 
20 
24 #include <deal.II/fe/fe_wedge_p.h>
25 #include <deal.II/fe/mapping_fe.h>
26 #include <deal.II/fe/mapping_q1.h>
28 
31 #include <deal.II/grid/tria.h>
32 
33 #include <memory>
34 
36 
37 namespace
38 {
39  namespace VTKCellType
40  {
41  // Define VTK constants for linear, quadratic and
42  // high-order Lagrange geometrices
43  enum
44  {
45  VTK_VERTEX = 1,
46  // Linear cells
47  VTK_LINE = 3,
48  VTK_TRIANGLE = 5,
49  VTK_QUAD = 9,
50  VTK_TETRA = 10,
51  VTK_HEXAHEDRON = 12,
52  VTK_WEDGE = 13,
53  VTK_PYRAMID = 14,
54  // Quadratic cells
55  VTK_QUADRATIC_EDGE = 21,
56  VTK_QUADRATIC_TRIANGLE = 22,
57  VTK_QUADRATIC_QUAD = 23,
58  VTK_QUADRATIC_TETRA = 24,
59  VTK_QUADRATIC_HEXAHEDRON = 25,
60  VTK_QUADRATIC_WEDGE = 26,
61  VTK_QUADRATIC_PYRAMID = 27,
62  // Lagrange cells
63  VTK_LAGRANGE_CURVE = 68,
64  VTK_LAGRANGE_TRIANGLE = 69,
65  VTK_LAGRANGE_QUADRILATERAL = 70,
66  VTK_LAGRANGE_TETRAHEDRON = 71,
67  VTK_LAGRANGE_HEXAHEDRON = 72,
68  VTK_LAGRANGE_WEDGE = 73,
69  VTK_LAGRANGE_PYRAMID = 74,
70  // Invalid code
71  VTK_INVALID = static_cast<unsigned int>(-1)
72  };
73 
74  } // namespace VTKCellType
75 
76 } // namespace
77 
78 
79 std::string
81 {
82  if (*this == ReferenceCells::Vertex)
83  return "Vertex";
84  else if (*this == ReferenceCells::Line)
85  return "Line";
86  else if (*this == ReferenceCells::Triangle)
87  return "Tri";
88  else if (*this == ReferenceCells::Quadrilateral)
89  return "Quad";
90  else if (*this == ReferenceCells::Tetrahedron)
91  return "Tet";
92  else if (*this == ReferenceCells::Pyramid)
93  return "Pyramid";
94  else if (*this == ReferenceCells::Wedge)
95  return "Wedge";
96  else if (*this == ReferenceCells::Hexahedron)
97  return "Hex";
98  else if (*this == ReferenceCells::Invalid)
99  return "Invalid";
100 
101  Assert(false, ExcNotImplemented());
102 
103  return "Invalid";
104 }
105 
106 
107 
108 template <int dim, int spacedim>
109 std::unique_ptr<Mapping<dim, spacedim>>
110 ReferenceCell::get_default_mapping(const unsigned int degree) const
111 {
112  AssertDimension(dim, get_dimension());
113 
114  if (is_hyper_cube())
115  return std::make_unique<MappingQGeneric<dim, spacedim>>(degree);
116  else if (is_simplex())
117  return std::make_unique<MappingFE<dim, spacedim>>(
119  else if (*this == ReferenceCells::Pyramid)
120  return std::make_unique<MappingFE<dim, spacedim>>(
122  else if (*this == ReferenceCells::Wedge)
123  return std::make_unique<MappingFE<dim, spacedim>>(
124  FE_WedgeP<dim, spacedim>(degree));
125  else
126  {
127  Assert(false, ExcNotImplemented());
128  }
129 
130  return std::make_unique<MappingQGeneric<dim, spacedim>>(degree);
131 }
132 
133 
134 
135 template <int dim, int spacedim>
138 {
139  AssertDimension(dim, get_dimension());
140 
141  if (is_hyper_cube())
142  {
144  }
145  else if (is_simplex())
146  {
147  static const MappingFE<dim, spacedim> mapping(
149  return mapping;
150  }
151  else if (*this == ReferenceCells::Pyramid)
152  {
153  static const MappingFE<dim, spacedim> mapping(
155  return mapping;
156  }
157  else if (*this == ReferenceCells::Wedge)
158  {
159  static const MappingFE<dim, spacedim> mapping(
161  return mapping;
162  }
163  else
164  {
165  Assert(false, ExcNotImplemented());
166  }
167 
168  return StaticMappingQ1<dim, spacedim>::mapping; // never reached
169 }
170 
171 
172 
173 template <int dim>
175 ReferenceCell::get_gauss_type_quadrature(const unsigned n_points_1D) const
176 {
177  AssertDimension(dim, get_dimension());
178 
179  if (is_hyper_cube())
180  return QGauss<dim>(n_points_1D);
181  else if (is_simplex())
182  return QGaussSimplex<dim>(n_points_1D);
183  else if (*this == ReferenceCells::Pyramid)
184  return QGaussPyramid<dim>(n_points_1D);
185  else if (*this == ReferenceCells::Wedge)
186  return QGaussWedge<dim>(n_points_1D);
187  else
188  Assert(false, ExcNotImplemented());
189 
190  return Quadrature<dim>(); // never reached
191 }
192 
193 
194 
195 template <int dim>
196 const Quadrature<dim> &
198 {
199  AssertDimension(dim, get_dimension());
200 
201  // A function that is used to fill a quadrature object of the
202  // desired type the first time we encounter a particular
203  // reference cell
204  const auto create_quadrature = [](const ReferenceCell &reference_cell) {
205  Triangulation<dim> tria;
207 
208  return Quadrature<dim>(tria.get_vertices());
209  };
210 
211  if (is_hyper_cube())
212  {
213  static const Quadrature<dim> quadrature = create_quadrature(*this);
214  return quadrature;
215  }
216  else if (is_simplex())
217  {
218  static const Quadrature<dim> quadrature = create_quadrature(*this);
219  return quadrature;
220  }
221  else if (*this == ReferenceCells::Pyramid)
222  {
223  static const Quadrature<dim> quadrature = create_quadrature(*this);
224  return quadrature;
225  }
226  else if (*this == ReferenceCells::Wedge)
227  {
228  static const Quadrature<dim> quadrature = create_quadrature(*this);
229  return quadrature;
230  }
231  else
232  Assert(false, ExcNotImplemented());
233 
234  static const Quadrature<dim> dummy;
235  return dummy; // never reached
236 }
237 
238 
239 
240 unsigned int
241 ReferenceCell::exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const
242 {
243  AssertIndexRange(vertex_n, n_vertices());
244 
245  if (*this == ReferenceCells::Line)
246  {
247  return vertex_n;
248  }
249  else if (*this == ReferenceCells::Triangle)
250  {
251  return vertex_n;
252  }
253  else if (*this == ReferenceCells::Quadrilateral)
254  {
255  constexpr std::array<unsigned int, 4> exodus_to_deal{{0, 1, 3, 2}};
256  return exodus_to_deal[vertex_n];
257  }
258  else if (*this == ReferenceCells::Tetrahedron)
259  {
260  return vertex_n;
261  }
262  else if (*this == ReferenceCells::Hexahedron)
263  {
264  constexpr std::array<unsigned int, 8> exodus_to_deal{
265  {0, 1, 3, 2, 4, 5, 7, 6}};
266  return exodus_to_deal[vertex_n];
267  }
268  else if (*this == ReferenceCells::Wedge)
269  {
270  constexpr std::array<unsigned int, 6> exodus_to_deal{{2, 1, 0, 5, 4, 3}};
271  return exodus_to_deal[vertex_n];
272  }
273  else if (*this == ReferenceCells::Pyramid)
274  {
275  constexpr std::array<unsigned int, 5> exodus_to_deal{{0, 1, 3, 2, 4}};
276  return exodus_to_deal[vertex_n];
277  }
278 
279  Assert(false, ExcNotImplemented());
280 
282 }
283 
284 
285 
286 unsigned int
287 ReferenceCell::exodusii_face_to_deal_face(const unsigned int face_n) const
288 {
289  AssertIndexRange(face_n, n_faces());
290 
291  if (*this == ReferenceCells::Vertex)
292  {
293  return 0;
294  }
295  if (*this == ReferenceCells::Line)
296  {
297  return face_n;
298  }
299  else if (*this == ReferenceCells::Triangle)
300  {
301  return face_n;
302  }
303  else if (*this == ReferenceCells::Quadrilateral)
304  {
305  constexpr std::array<unsigned int, 4> exodus_to_deal{{2, 1, 3, 0}};
306  return exodus_to_deal[face_n];
307  }
308  else if (*this == ReferenceCells::Tetrahedron)
309  {
310  constexpr std::array<unsigned int, 4> exodus_to_deal{{1, 3, 2, 0}};
311  return exodus_to_deal[face_n];
312  }
313  else if (*this == ReferenceCells::Hexahedron)
314  {
315  constexpr std::array<unsigned int, 6> exodus_to_deal{{2, 1, 3, 0, 4, 5}};
316  return exodus_to_deal[face_n];
317  }
318  else if (*this == ReferenceCells::Wedge)
319  {
320  constexpr std::array<unsigned int, 6> exodus_to_deal{{3, 4, 2, 0, 1}};
321  return exodus_to_deal[face_n];
322  }
323  else if (*this == ReferenceCells::Pyramid)
324  {
325  constexpr std::array<unsigned int, 5> exodus_to_deal{{3, 2, 4, 1, 0}};
326  return exodus_to_deal[face_n];
327  }
328 
329  Assert(false, ExcNotImplemented());
330 
332 }
333 
334 
335 
336 unsigned int
337 ReferenceCell::unv_vertex_to_deal_vertex(const unsigned int vertex_n) const
338 {
339  AssertIndexRange(vertex_n, n_vertices());
340  // Information on this file format isn't easy to find - the documents here
341  //
342  // https://www.ceas3.uc.edu/sdrluff/
343  //
344  // Don't actually explain anything about the sections we care about (2412) in
345  // any detail. For node numbering I worked backwards from what is actually in
346  // our test files (since that's supposed to work), which all use some
347  // non-standard clockwise numbering scheme which starts at the bottom right
348  // vertex.
349  if (*this == ReferenceCells::Line)
350  {
351  return vertex_n;
352  }
353  else if (*this == ReferenceCells::Quadrilateral)
354  {
355  constexpr std::array<unsigned int, 4> unv_to_deal{{1, 0, 2, 3}};
356  return unv_to_deal[vertex_n];
357  }
358  else if (*this == ReferenceCells::Hexahedron)
359  {
360  constexpr std::array<unsigned int, 8> unv_to_deal{
361  {6, 7, 5, 4, 2, 3, 1, 0}};
362  return unv_to_deal[vertex_n];
363  }
364 
365  Assert(false, ExcNotImplemented());
366 
368 }
369 
370 
371 
372 unsigned int
374 {
375  if (*this == ReferenceCells::Vertex)
376  return VTKCellType::VTK_VERTEX;
377  else if (*this == ReferenceCells::Line)
378  return VTKCellType::VTK_LINE;
379  else if (*this == ReferenceCells::Triangle)
380  return VTKCellType::VTK_TRIANGLE;
381  else if (*this == ReferenceCells::Quadrilateral)
382  return VTKCellType::VTK_QUAD;
383  else if (*this == ReferenceCells::Tetrahedron)
384  return VTKCellType::VTK_TETRA;
385  else if (*this == ReferenceCells::Pyramid)
386  return VTKCellType::VTK_PYRAMID;
387  else if (*this == ReferenceCells::Wedge)
388  return VTKCellType::VTK_WEDGE;
389  else if (*this == ReferenceCells::Hexahedron)
390  return VTKCellType::VTK_HEXAHEDRON;
391  else if (*this == ReferenceCells::Invalid)
392  return VTKCellType::VTK_INVALID;
393 
394  Assert(false, ExcNotImplemented());
395 
396  return VTKCellType::VTK_INVALID;
397 }
398 
399 
400 
401 unsigned int
403 {
404  if (*this == ReferenceCells::Vertex)
405  return VTKCellType::VTK_VERTEX;
406  else if (*this == ReferenceCells::Line)
407  return VTKCellType::VTK_QUADRATIC_EDGE;
408  else if (*this == ReferenceCells::Triangle)
409  return VTKCellType::VTK_QUADRATIC_TRIANGLE;
410  else if (*this == ReferenceCells::Quadrilateral)
411  return VTKCellType::VTK_QUADRATIC_QUAD;
412  else if (*this == ReferenceCells::Tetrahedron)
413  return VTKCellType::VTK_QUADRATIC_TETRA;
414  else if (*this == ReferenceCells::Pyramid)
415  return VTKCellType::VTK_QUADRATIC_PYRAMID;
416  else if (*this == ReferenceCells::Wedge)
417  return VTKCellType::VTK_QUADRATIC_WEDGE;
418  else if (*this == ReferenceCells::Hexahedron)
419  return VTKCellType::VTK_QUADRATIC_HEXAHEDRON;
420  else if (*this == ReferenceCells::Invalid)
421  return VTKCellType::VTK_INVALID;
422 
423  Assert(false, ExcNotImplemented());
424 
425  return VTKCellType::VTK_INVALID;
426 }
427 
428 
429 
430 unsigned int
432 {
433  if (*this == ReferenceCells::Vertex)
434  return VTKCellType::VTK_VERTEX;
435  else if (*this == ReferenceCells::Line)
436  return VTKCellType::VTK_LAGRANGE_CURVE;
437  else if (*this == ReferenceCells::Triangle)
438  return VTKCellType::VTK_LAGRANGE_TRIANGLE;
439  else if (*this == ReferenceCells::Quadrilateral)
440  return VTKCellType::VTK_LAGRANGE_QUADRILATERAL;
441  else if (*this == ReferenceCells::Tetrahedron)
442  return VTKCellType::VTK_LAGRANGE_TETRAHEDRON;
443  else if (*this == ReferenceCells::Pyramid)
444  return VTKCellType::VTK_LAGRANGE_PYRAMID;
445  else if (*this == ReferenceCells::Wedge)
446  return VTKCellType::VTK_LAGRANGE_WEDGE;
447  else if (*this == ReferenceCells::Hexahedron)
448  return VTKCellType::VTK_LAGRANGE_HEXAHEDRON;
449  else if (*this == ReferenceCells::Invalid)
450  return VTKCellType::VTK_INVALID;
451 
452  Assert(false, ExcNotImplemented());
453 
454  return VTKCellType::VTK_INVALID;
455 }
456 
457 #include "reference_cell.inst"
458 
constexpr const ReferenceCell Pyramid
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
unsigned int vtk_quadratic_type() const
unsigned int exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const
unsigned int exodusii_face_to_deal_face(const unsigned int face_n) const
unsigned int vtk_lagrange_type() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
constexpr const ReferenceCell Triangle
std::string to_string() const
const Quadrature< dim > & get_nodal_type_quadrature() const
constexpr const ReferenceCell Line
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Tetrahedron
#define Assert(cond, exc)
Definition: exceptions.h:1465
Abstract base class for mapping classes.
Definition: mapping.h:303
Quadrature< dim > get_gauss_type_quadrature(const unsigned n_points_1D) const
const std::vector< Point< spacedim > > & get_vertices() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
constexpr const ReferenceCell Hexahedron
const Mapping< dim, spacedim > & get_default_linear_mapping() const
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
constexpr const ReferenceCell Invalid
std::unique_ptr< Mapping< dim, spacedim > > get_default_mapping(const unsigned int degree) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
constexpr const ReferenceCell Quadrilateral
constexpr const ReferenceCell Vertex
static ::ExceptionBase & ExcNotImplemented()
unsigned int vtk_linear_type() const
unsigned int unv_vertex_to_deal_vertex(const unsigned int vertex_n) const