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\}}\)
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
28
31#include <deal.II/grid/tria.h>
32
33#include <iostream>
34#include <memory>
35
37
38namespace
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
80std::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
109template <int dim, int spacedim>
110std::unique_ptr<Mapping<dim, spacedim>>
111ReferenceCell::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>>(
126 else
127 {
128 Assert(false, ExcNotImplemented());
129 }
130
131 return std::make_unique<MappingQ<dim, spacedim>>(degree);
132}
133
134
135
136template <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
174template <int dim>
176ReferenceCell::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
196template <int dim>
197const 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
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
242unsigned int
243ReferenceCell::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
288unsigned int
289ReferenceCell::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
338unsigned int
339ReferenceCell::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
374unsigned 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
403unsigned 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
432unsigned 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
461unsigned 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
546std::ostream &
547operator<<(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
559std::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(
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:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Point< 3 > vertices[4]
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
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