Reference documentation for deal.II version GIT 3801df8983 2022-05-22 23:30:02+00:00
\(\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 - 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 
20 
24 #include <deal.II/fe/fe_wedge_p.h>
25 #include <deal.II/fe/mapping_fe.h>
26 #include <deal.II/fe/mapping_q.h>
27 #include <deal.II/fe/mapping_q1.h>
28 
31 #include <deal.II/grid/tria.h>
32 
33 #include <iostream>
34 #include <memory>
35 
37 
38 namespace
39 {
40  namespace VTKCellType
41  {
42  // Define VTK constants for linear, quadratic and
43  // high-order Lagrange geometrices
44  enum
45  {
46  VTK_VERTEX = 1,
47  // Linear cells
48  VTK_LINE = 3,
49  VTK_TRIANGLE = 5,
50  VTK_QUAD = 9,
51  VTK_TETRA = 10,
52  VTK_HEXAHEDRON = 12,
53  VTK_WEDGE = 13,
54  VTK_PYRAMID = 14,
55  // Quadratic cells
56  VTK_QUADRATIC_EDGE = 21,
57  VTK_QUADRATIC_TRIANGLE = 22,
58  VTK_QUADRATIC_QUAD = 23,
59  VTK_QUADRATIC_TETRA = 24,
60  VTK_QUADRATIC_HEXAHEDRON = 25,
61  VTK_QUADRATIC_WEDGE = 26,
62  VTK_QUADRATIC_PYRAMID = 27,
63  // Lagrange cells
64  VTK_LAGRANGE_CURVE = 68,
65  VTK_LAGRANGE_TRIANGLE = 69,
66  VTK_LAGRANGE_QUADRILATERAL = 70,
67  VTK_LAGRANGE_TETRAHEDRON = 71,
68  VTK_LAGRANGE_HEXAHEDRON = 72,
69  VTK_LAGRANGE_WEDGE = 73,
70  VTK_LAGRANGE_PYRAMID = 74,
71  // Invalid code
72  VTK_INVALID = static_cast<unsigned int>(-1)
73  };
74 
75  } // namespace VTKCellType
76 
77 } // namespace
78 
79 
80 std::string
82 {
83  if (*this == ReferenceCells::Vertex)
84  return "Vertex";
85  else if (*this == ReferenceCells::Line)
86  return "Line";
87  else if (*this == ReferenceCells::Triangle)
88  return "Tri";
89  else if (*this == ReferenceCells::Quadrilateral)
90  return "Quad";
91  else if (*this == ReferenceCells::Tetrahedron)
92  return "Tet";
93  else if (*this == ReferenceCells::Pyramid)
94  return "Pyramid";
95  else if (*this == ReferenceCells::Wedge)
96  return "Wedge";
97  else if (*this == ReferenceCells::Hexahedron)
98  return "Hex";
99  else if (*this == ReferenceCells::Invalid)
100  return "Invalid";
101 
102  Assert(false, ExcNotImplemented());
103 
104  return "Invalid";
105 }
106 
107 
108 
109 template <int dim, int spacedim>
110 std::unique_ptr<Mapping<dim, spacedim>>
111 ReferenceCell::get_default_mapping(const unsigned int degree) const
112 {
114 
115  if (is_hyper_cube())
116  return std::make_unique<MappingQ<dim, spacedim>>(degree);
117  else if (is_simplex())
118  return std::make_unique<MappingFE<dim, spacedim>>(
120  else if (*this == ReferenceCells::Pyramid)
121  return std::make_unique<MappingFE<dim, spacedim>>(
123  else if (*this == ReferenceCells::Wedge)
124  return std::make_unique<MappingFE<dim, spacedim>>(
125  FE_WedgeP<dim, spacedim>(degree));
126  else
127  {
128  Assert(false, ExcNotImplemented());
129  }
130 
131  return std::make_unique<MappingQ<dim, spacedim>>(degree);
132 }
133 
134 
135 
136 template <int dim, int spacedim>
139 {
141 
142  if (is_hyper_cube())
143  {
145  }
146  else if (is_simplex())
147  {
148  static const MappingFE<dim, spacedim> mapping(
150  return mapping;
151  }
152  else if (*this == ReferenceCells::Pyramid)
153  {
154  static const MappingFE<dim, spacedim> mapping(
156  return mapping;
157  }
158  else if (*this == ReferenceCells::Wedge)
159  {
160  static const MappingFE<dim, spacedim> mapping(
162  return mapping;
163  }
164  else
165  {
166  Assert(false, ExcNotImplemented());
167  }
168 
169  return StaticMappingQ1<dim, spacedim>::mapping; // never reached
170 }
171 
172 
173 
174 template <int dim>
176 ReferenceCell::get_gauss_type_quadrature(const unsigned n_points_1D) const
177 {
179 
180  if (is_hyper_cube())
181  return QGauss<dim>(n_points_1D);
182  else if (is_simplex())
183  return QGaussSimplex<dim>(n_points_1D);
184  else if (*this == ReferenceCells::Pyramid)
185  return QGaussPyramid<dim>(n_points_1D);
186  else if (*this == ReferenceCells::Wedge)
187  return QGaussWedge<dim>(n_points_1D);
188  else
189  Assert(false, ExcNotImplemented());
190 
191  return Quadrature<dim>(); // never reached
192 }
193 
194 
195 
196 template <int dim>
197 const Quadrature<dim> &
199 {
201 
202  // A function that is used to fill a quadrature object of the
203  // desired type the first time we encounter a particular
204  // reference cell
205  const auto create_quadrature = [](const ReferenceCell &reference_cell) {
206  std::vector<Point<dim>> vertices(reference_cell.n_vertices());
207  for (const unsigned int v : reference_cell.vertex_indices())
208  vertices[v] = reference_cell.vertex<dim>(v);
209 
210  return Quadrature<dim>(vertices);
211  };
212 
213  if (is_hyper_cube())
214  {
215  static const Quadrature<dim> quadrature = create_quadrature(*this);
216  return quadrature;
217  }
218  else if (is_simplex())
219  {
220  static const Quadrature<dim> quadrature = create_quadrature(*this);
221  return quadrature;
222  }
223  else if (*this == ReferenceCells::Pyramid)
224  {
225  static const Quadrature<dim> quadrature = create_quadrature(*this);
226  return quadrature;
227  }
228  else if (*this == ReferenceCells::Wedge)
229  {
230  static const Quadrature<dim> quadrature = create_quadrature(*this);
231  return quadrature;
232  }
233  else
234  Assert(false, ExcNotImplemented());
235 
236  static const Quadrature<dim> dummy;
237  return dummy; // never reached
238 }
239 
240 
241 
242 unsigned int
243 ReferenceCell::exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const
244 {
245  AssertIndexRange(vertex_n, n_vertices());
246 
247  if (*this == ReferenceCells::Line)
248  {
249  return vertex_n;
250  }
251  else if (*this == ReferenceCells::Triangle)
252  {
253  return vertex_n;
254  }
255  else if (*this == ReferenceCells::Quadrilateral)
256  {
257  constexpr std::array<unsigned int, 4> exodus_to_deal{{0, 1, 3, 2}};
258  return exodus_to_deal[vertex_n];
259  }
260  else if (*this == ReferenceCells::Tetrahedron)
261  {
262  return vertex_n;
263  }
264  else if (*this == ReferenceCells::Hexahedron)
265  {
266  constexpr std::array<unsigned int, 8> exodus_to_deal{
267  {0, 1, 3, 2, 4, 5, 7, 6}};
268  return exodus_to_deal[vertex_n];
269  }
270  else if (*this == ReferenceCells::Wedge)
271  {
272  constexpr std::array<unsigned int, 6> exodus_to_deal{{2, 1, 0, 5, 4, 3}};
273  return exodus_to_deal[vertex_n];
274  }
275  else if (*this == ReferenceCells::Pyramid)
276  {
277  constexpr std::array<unsigned int, 5> exodus_to_deal{{0, 1, 3, 2, 4}};
278  return exodus_to_deal[vertex_n];
279  }
280 
281  Assert(false, ExcNotImplemented());
282 
284 }
285 
286 
287 
288 unsigned int
289 ReferenceCell::exodusii_face_to_deal_face(const unsigned int face_n) const
290 {
291  AssertIndexRange(face_n, n_faces());
292 
293  if (*this == ReferenceCells::Vertex)
294  {
295  return 0;
296  }
297  if (*this == ReferenceCells::Line)
298  {
299  return face_n;
300  }
301  else if (*this == ReferenceCells::Triangle)
302  {
303  return face_n;
304  }
305  else if (*this == ReferenceCells::Quadrilateral)
306  {
307  constexpr std::array<unsigned int, 4> exodus_to_deal{{2, 1, 3, 0}};
308  return exodus_to_deal[face_n];
309  }
310  else if (*this == ReferenceCells::Tetrahedron)
311  {
312  constexpr std::array<unsigned int, 4> exodus_to_deal{{1, 3, 2, 0}};
313  return exodus_to_deal[face_n];
314  }
315  else if (*this == ReferenceCells::Hexahedron)
316  {
317  constexpr std::array<unsigned int, 6> exodus_to_deal{{2, 1, 3, 0, 4, 5}};
318  return exodus_to_deal[face_n];
319  }
320  else if (*this == ReferenceCells::Wedge)
321  {
322  constexpr std::array<unsigned int, 6> exodus_to_deal{{3, 4, 2, 0, 1}};
323  return exodus_to_deal[face_n];
324  }
325  else if (*this == ReferenceCells::Pyramid)
326  {
327  constexpr std::array<unsigned int, 5> exodus_to_deal{{3, 2, 4, 1, 0}};
328  return exodus_to_deal[face_n];
329  }
330 
331  Assert(false, ExcNotImplemented());
332 
334 }
335 
336 
337 
338 unsigned int
339 ReferenceCell::unv_vertex_to_deal_vertex(const unsigned int vertex_n) const
340 {
341  AssertIndexRange(vertex_n, n_vertices());
342  // Information on this file format isn't easy to find - the documents here
343  //
344  // https://www.ceas3.uc.edu/sdrluff/
345  //
346  // Don't actually explain anything about the sections we care about (2412) in
347  // any detail. For node numbering I worked backwards from what is actually in
348  // our test files (since that's supposed to work), which all use some
349  // non-standard clockwise numbering scheme which starts at the bottom right
350  // vertex.
351  if (*this == ReferenceCells::Line)
352  {
353  return vertex_n;
354  }
355  else if (*this == ReferenceCells::Quadrilateral)
356  {
357  constexpr std::array<unsigned int, 4> unv_to_deal{{1, 0, 2, 3}};
358  return unv_to_deal[vertex_n];
359  }
360  else if (*this == ReferenceCells::Hexahedron)
361  {
362  constexpr std::array<unsigned int, 8> unv_to_deal{
363  {6, 7, 5, 4, 2, 3, 1, 0}};
364  return unv_to_deal[vertex_n];
365  }
366 
367  Assert(false, ExcNotImplemented());
368 
370 }
371 
372 
373 
374 unsigned int
376 {
377  if (*this == ReferenceCells::Vertex)
378  return VTKCellType::VTK_VERTEX;
379  else if (*this == ReferenceCells::Line)
380  return VTKCellType::VTK_LINE;
381  else if (*this == ReferenceCells::Triangle)
382  return VTKCellType::VTK_TRIANGLE;
383  else if (*this == ReferenceCells::Quadrilateral)
384  return VTKCellType::VTK_QUAD;
385  else if (*this == ReferenceCells::Tetrahedron)
386  return VTKCellType::VTK_TETRA;
387  else if (*this == ReferenceCells::Pyramid)
388  return VTKCellType::VTK_PYRAMID;
389  else if (*this == ReferenceCells::Wedge)
390  return VTKCellType::VTK_WEDGE;
391  else if (*this == ReferenceCells::Hexahedron)
392  return VTKCellType::VTK_HEXAHEDRON;
393  else if (*this == ReferenceCells::Invalid)
394  return VTKCellType::VTK_INVALID;
395 
396  Assert(false, ExcNotImplemented());
397 
398  return VTKCellType::VTK_INVALID;
399 }
400 
401 
402 
403 unsigned int
405 {
406  if (*this == ReferenceCells::Vertex)
407  return VTKCellType::VTK_VERTEX;
408  else if (*this == ReferenceCells::Line)
409  return VTKCellType::VTK_QUADRATIC_EDGE;
410  else if (*this == ReferenceCells::Triangle)
411  return VTKCellType::VTK_QUADRATIC_TRIANGLE;
412  else if (*this == ReferenceCells::Quadrilateral)
413  return VTKCellType::VTK_QUADRATIC_QUAD;
414  else if (*this == ReferenceCells::Tetrahedron)
415  return VTKCellType::VTK_QUADRATIC_TETRA;
416  else if (*this == ReferenceCells::Pyramid)
417  return VTKCellType::VTK_QUADRATIC_PYRAMID;
418  else if (*this == ReferenceCells::Wedge)
419  return VTKCellType::VTK_QUADRATIC_WEDGE;
420  else if (*this == ReferenceCells::Hexahedron)
421  return VTKCellType::VTK_QUADRATIC_HEXAHEDRON;
422  else if (*this == ReferenceCells::Invalid)
423  return VTKCellType::VTK_INVALID;
424 
425  Assert(false, ExcNotImplemented());
426 
427  return VTKCellType::VTK_INVALID;
428 }
429 
430 
431 
432 unsigned int
434 {
435  if (*this == ReferenceCells::Vertex)
436  return VTKCellType::VTK_VERTEX;
437  else if (*this == ReferenceCells::Line)
438  return VTKCellType::VTK_LAGRANGE_CURVE;
439  else if (*this == ReferenceCells::Triangle)
440  return VTKCellType::VTK_LAGRANGE_TRIANGLE;
441  else if (*this == ReferenceCells::Quadrilateral)
442  return VTKCellType::VTK_LAGRANGE_QUADRILATERAL;
443  else if (*this == ReferenceCells::Tetrahedron)
444  return VTKCellType::VTK_LAGRANGE_TETRAHEDRON;
445  else if (*this == ReferenceCells::Pyramid)
446  return VTKCellType::VTK_LAGRANGE_PYRAMID;
447  else if (*this == ReferenceCells::Wedge)
448  return VTKCellType::VTK_LAGRANGE_WEDGE;
449  else if (*this == ReferenceCells::Hexahedron)
450  return VTKCellType::VTK_LAGRANGE_HEXAHEDRON;
451  else if (*this == ReferenceCells::Invalid)
452  return VTKCellType::VTK_INVALID;
453 
454  Assert(false, ExcNotImplemented());
455 
456  return VTKCellType::VTK_INVALID;
457 }
458 
459 
460 
461 unsigned int
463 {
464  /*
465  From the GMSH documentation:
466 
467  elm-type
468  defines the geometrical type of the n-th element:
469 
470  1
471  Line (2 nodes).
472 
473  2
474  Triangle (3 nodes).
475 
476  3
477  Quadrangle (4 nodes).
478 
479  4
480  Tetrahedron (4 nodes).
481 
482  5
483  Hexahedron (8 nodes).
484 
485  6
486  Prism (6 nodes).
487 
488  7
489  Pyramid (5 nodes).
490 
491  8
492  Second order line (3 nodes: 2 associated with the vertices and 1 with the
493  edge).
494 
495  9
496  Second order triangle (6 nodes: 3 associated with the vertices and 3 with
497  the edges).
498 
499  10 Second order quadrangle (9 nodes: 4 associated with the
500  vertices, 4 with the edges and 1 with the face).
501 
502  11 Second order tetrahedron (10 nodes: 4 associated with the vertices and 6
503  with the edges).
504 
505  12 Second order hexahedron (27 nodes: 8 associated with the vertices, 12
506  with the edges, 6 with the faces and 1 with the volume).
507 
508  13 Second order prism (18 nodes: 6 associated with the vertices, 9 with the
509  edges and 3 with the quadrangular faces).
510 
511  14 Second order pyramid (14 nodes: 5 associated with the vertices, 8 with
512  the edges and 1 with the quadrangular face).
513 
514  15 Point (1 node).
515  */
516 
517  if (*this == ReferenceCells::Vertex)
518  return 15;
519  else if (*this == ReferenceCells::Line)
520  return 1;
521  else if (*this == ReferenceCells::Triangle)
522  return 2;
523  else if (*this == ReferenceCells::Quadrilateral)
524  return 3;
525  else if (*this == ReferenceCells::Tetrahedron)
526  return 4;
527  else if (*this == ReferenceCells::Pyramid)
528  return 7;
529  else if (*this == ReferenceCells::Wedge)
530  return 6;
531  else if (*this == ReferenceCells::Hexahedron)
532  return 5;
533  else if (*this == ReferenceCells::Invalid)
534  {
535  Assert(false, ExcNotImplemented());
537  }
538 
539  Assert(false, ExcNotImplemented());
540 
542 }
543 
544 
545 
546 std::ostream &
547 operator<<(std::ostream &out, const ReferenceCell &reference_cell)
548 {
549  AssertThrow(out.fail() == false, ExcIO());
550 
551  // Output as an integer to avoid outputting it as a character with
552  // potentially non-printing value:
553  out << static_cast<unsigned int>(reference_cell.kind);
554  return out;
555 }
556 
557 
558 
559 std::istream &
561 {
562  AssertThrow(in.fail() == false, ExcIO());
563 
564  // Read the information as an integer and convert it to the correct type
565  unsigned int value;
566  in >> value;
567  reference_cell.kind = static_cast<decltype(reference_cell.kind)>(value);
568 
569  // Ensure that the object we read is valid
570  Assert(
580  ExcMessage(
581  "The reference cell kind just read does not correspond to one of the valid choices. There must be an error."));
582 
583  return in;
584 }
585 
586 
587 
588 #include "reference_cell.inst"
589 
OutputOperator< VectorType > & operator<<(OutputOperator< VectorType > &out, unsigned int step)
Definition: operator.h:165
std::istream & operator>>(std::istream &in, Point< dim, Number > &p)
Definition: point.h:688
unsigned int n_vertices() const
bool is_hyper_cube() const
unsigned int vtk_linear_type() const
Quadrature< dim > get_gauss_type_quadrature(const unsigned n_points_1D) const
unsigned int gmsh_element_type() const
unsigned int n_faces() const
unsigned int exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const
unsigned int unv_vertex_to_deal_vertex(const unsigned int vertex_n) const
std::unique_ptr< Mapping< dim, spacedim > > get_default_mapping(const unsigned int degree) const
unsigned int vtk_quadratic_type() const
unsigned int vtk_lagrange_type() const
unsigned int get_dimension() const
unsigned int exodusii_face_to_deal_face(const unsigned int face_n) const
const Mapping< dim, spacedim > & get_default_linear_mapping() const
std::string to_string() const
bool is_simplex() const
const Quadrature< dim > & get_nodal_type_quadrature() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:416
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:417
Point< 3 > vertices[4]
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Quadrilateral
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Pyramid
constexpr const ReferenceCell Invalid
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Hexahedron
constexpr const ReferenceCell Vertex
constexpr const ReferenceCell Line
static const unsigned int invalid_unsigned_int
Definition: types.h:201