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,
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,
119template <
int dim,
int spacedim>
120std::unique_ptr<Mapping<dim, spacedim>>
126 return std::make_unique<MappingQ<dim, spacedim>>(degree);
131 return std::make_unique<MappingFE<dim, spacedim>>(
134 return std::make_unique<MappingFE<dim, spacedim>>(
141 return std::make_unique<MappingQ<dim, spacedim>>(degree);
146template <
int dim,
int spacedim>
215 const auto create_quadrature = [](
const ReferenceCell &reference_cell) {
216 std::vector<Point<dim>>
vertices(reference_cell.n_vertices());
217 for (
const unsigned int v : reference_cell.vertex_indices())
218 vertices[v] = reference_cell.vertex<dim>(v);
264 constexpr std::array<unsigned int, 4> exodus_to_deal{{0, 1, 3, 2}};
265 return exodus_to_deal[vertex_n];
271 constexpr std::array<unsigned int, 8> exodus_to_deal{
272 {0, 1, 3, 2, 4, 5, 7, 6}};
273 return exodus_to_deal[vertex_n];
277 constexpr std::array<unsigned int, 6> exodus_to_deal{
279 return exodus_to_deal[vertex_n];
283 constexpr std::array<unsigned int, 5> exodus_to_deal{{0, 1, 3, 2, 4}};
284 return exodus_to_deal[vertex_n];
309 constexpr std::array<unsigned int, 4> exodus_to_deal{{2, 1, 3, 0}};
310 return exodus_to_deal[face_n];
314 constexpr std::array<unsigned int, 4> exodus_to_deal{{1, 3, 2, 0}};
315 return exodus_to_deal[face_n];
319 constexpr std::array<unsigned int, 6> exodus_to_deal{
321 return exodus_to_deal[face_n];
325 constexpr std::array<unsigned int, 6> exodus_to_deal{{3, 4, 2, 0, 1}};
326 return exodus_to_deal[face_n];
330 constexpr std::array<unsigned int, 5> exodus_to_deal{{3, 2, 4, 1, 0}};
331 return exodus_to_deal[face_n];
361 constexpr std::array<unsigned int, 4> unv_to_deal{{1, 0, 2, 3}};
362 return unv_to_deal[vertex_n];
366 constexpr std::array<unsigned int, 8> unv_to_deal{
367 {6, 7, 5, 4, 2, 3, 1, 0}};
368 return unv_to_deal[vertex_n];
384 return VTKCellType::VTK_VERTEX;
386 return VTKCellType::VTK_LINE;
388 return VTKCellType::VTK_TRIANGLE;
390 return VTKCellType::VTK_QUAD;
392 return VTKCellType::VTK_TETRA;
394 return VTKCellType::VTK_PYRAMID;
396 return VTKCellType::VTK_WEDGE;
398 return VTKCellType::VTK_HEXAHEDRON;
400 return VTKCellType::VTK_INVALID;
405 return VTKCellType::VTK_INVALID;
416 return VTKCellType::VTK_VERTEX;
418 return VTKCellType::VTK_QUADRATIC_EDGE;
420 return VTKCellType::VTK_QUADRATIC_TRIANGLE;
422 return VTKCellType::VTK_QUADRATIC_QUAD;
424 return VTKCellType::VTK_QUADRATIC_TETRA;
426 return VTKCellType::VTK_QUADRATIC_PYRAMID;
428 return VTKCellType::VTK_QUADRATIC_WEDGE;
430 return VTKCellType::VTK_QUADRATIC_HEXAHEDRON;
432 return VTKCellType::VTK_INVALID;
437 return VTKCellType::VTK_INVALID;
448 return VTKCellType::VTK_VERTEX;
450 return VTKCellType::VTK_LAGRANGE_CURVE;
452 return VTKCellType::VTK_LAGRANGE_TRIANGLE;
454 return VTKCellType::VTK_LAGRANGE_QUADRILATERAL;
456 return VTKCellType::VTK_LAGRANGE_TETRAHEDRON;
458 return VTKCellType::VTK_LAGRANGE_PYRAMID;
460 return VTKCellType::VTK_LAGRANGE_WEDGE;
462 return VTKCellType::VTK_LAGRANGE_HEXAHEDRON;
464 return VTKCellType::VTK_INVALID;
469 return VTKCellType::VTK_INVALID;
476ReferenceCell::vtk_lexicographic_to_node_index<0>(
477 const std::array<unsigned, 0> &,
478 const std::array<unsigned, 0> &,
489ReferenceCell::vtk_lexicographic_to_node_index<1>(
490 const std::array<unsigned, 1> &,
491 const std::array<unsigned, 1> &,
506ReferenceCell::vtk_lexicographic_to_node_index<2>(
507 const std::array<unsigned, 2> &node_indices,
508 const std::array<unsigned, 2> &nodes_per_direction,
513 const unsigned int i = node_indices[0];
514 const unsigned int j = node_indices[1];
516 const bool ibdy = (i == 0 || i == nodes_per_direction[0]);
517 const bool jbdy = (j == 0 || j == nodes_per_direction[1]);
519 const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0);
523 return (i != 0u ? (j != 0u ? 2 : 1) : (j != 0u ? 3 : 0));
533 nodes_per_direction[0] - 1 + nodes_per_direction[1] - 1 :
541 (i != 0u ? nodes_per_direction[0] - 1 :
542 2 * (nodes_per_direction[0] - 1) +
543 nodes_per_direction[1] - 1) +
548 offset += 2 * (nodes_per_direction[0] - 1 + nodes_per_direction[1] - 1);
550 return offset + (i - 1) + (nodes_per_direction[0] - 1) * ((j - 1));
567ReferenceCell::vtk_lexicographic_to_node_index<3>(
568 const std::array<unsigned, 3> &node_indices,
569 const std::array<unsigned, 3> &nodes_per_direction,
570 const bool legacy_format)
const
574 const unsigned int i = node_indices[0];
575 const unsigned int j = node_indices[1];
576 const unsigned int k = node_indices[2];
578 const bool ibdy = (i == 0 || i == nodes_per_direction[0]);
579 const bool jbdy = (j == 0 || j == nodes_per_direction[1]);
580 const bool kbdy = (k == 0 || k == nodes_per_direction[2]);
582 const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0) + (kbdy ? 1 : 0);
586 return (i != 0u ? (j != 0u ? 2 : 1) : (j != 0u ? 3 : 0)) +
597 nodes_per_direction[0] - 1 + nodes_per_direction[1] - 1 :
599 (k != 0u ? 2 * (nodes_per_direction[0] - 1 +
600 nodes_per_direction[1] - 1) :
607 (i != 0u ? nodes_per_direction[0] - 1 :
608 2 * (nodes_per_direction[0] - 1) +
609 nodes_per_direction[1] - 1) +
610 (k != 0u ? 2 * (nodes_per_direction[0] - 1 +
611 nodes_per_direction[1] - 1) :
617 4 * (nodes_per_direction[0] - 1) + 4 * (nodes_per_direction[1] - 1);
620 (nodes_per_direction[2] - 1) *
621 (i != 0u ? (j != 0u ? 3 : 1) : (j != 0u ? 2 : 0)) +
625 (nodes_per_direction[2] - 1) *
626 (i != 0u ? (j != 0u ? 2 : 1) : (j != 0u ? 3 : 0)) +
630 offset += 4 * (nodes_per_direction[0] - 1 + nodes_per_direction[1] - 1 +
631 nodes_per_direction[2] - 1);
636 return (j - 1) + ((nodes_per_direction[1] - 1) * (k - 1)) +
637 (i != 0u ? (nodes_per_direction[1] - 1) *
638 (nodes_per_direction[2] - 1) :
642 offset += 2 * (nodes_per_direction[1] - 1) * (nodes_per_direction[2] - 1);
645 return (i - 1) + ((nodes_per_direction[0] - 1) * (k - 1)) +
646 (j != 0u ? (nodes_per_direction[2] - 1) *
647 (nodes_per_direction[0] - 1) :
651 offset += 2 * (nodes_per_direction[2] - 1) * (nodes_per_direction[0] - 1);
653 return (i - 1) + ((nodes_per_direction[0] - 1) * (j - 1)) +
655 (nodes_per_direction[0] - 1) * (nodes_per_direction[1] - 1) :
661 offset += 2 * ((nodes_per_direction[1] - 1) * (nodes_per_direction[2] - 1) +
662 (nodes_per_direction[2] - 1) * (nodes_per_direction[0] - 1) +
663 (nodes_per_direction[0] - 1) * (nodes_per_direction[1] - 1));
664 return offset + (i - 1) +
665 (nodes_per_direction[0] - 1) *
666 ((j - 1) + (nodes_per_direction[1] - 1) * ((k - 1)));
692 static constexpr std::array<unsigned int, 4> index_translation_table =
694 return index_translation_table[vertex_index];
700 static constexpr std::array<unsigned int, 5> index_translation_table =
702 return index_translation_table[vertex_index];
708 static constexpr std::array<unsigned int, 8> index_translation_table =
709 {{0, 1, 3, 2, 4, 5, 7, 6}};
710 return index_translation_table[vertex_index];
815 std::pair<Point<dim>,
double>
826 const double t = ((x1 - x0) * (p - x0)) / ((x1 - x0).norm_square());
833 const auto p2 = x0 + t * (x1 - x0);
834 return std::make_pair(p2, p2.distance_square(p));
842 std::pair<Point<dim>,
double>
843 project_to_quad(
const std::array<
Point<dim>, 3> & ,
849 std::numeric_limits<double>::signaling_NaN());
870 std::pair<Point<3>,
double>
882 std::array<Point<3>, 3> shifted_vertices =
vertices;
884 for (
Point<3> &shifted_vertex : shifted_vertices)
885 shifted_vertex +=
shift;
894 const Point<3> vertex = shifted_vertices[0];
904 e0 = shifted_vertices[1] - shifted_vertices[0];
905 e1 = shifted_vertices[2] - shifted_vertices[0];
910 e0 = shifted_vertices[1] - shifted_vertices[0];
911 e1 = shifted_vertices[2] - shifted_vertices[0];
916 const double c0 = e0 * (shifted_p - vertex) / e0.
norm_square();
917 const double c1 = e1 * (shifted_p - vertex) / e1.norm_square();
918 const Point<3> projected_shifted_p = vertex + c0 * e0 + c1 * e1;
920 bool in_quad =
false;
924 for (
unsigned int i = 0; i < 3; ++i)
925 shifted_vertex_matrix[i] = shifted_vertices[i];
928 bool is_convex_combination =
true;
929 for (
unsigned int i = 0; i < 3; ++i)
930 is_convex_combination = is_convex_combination &&
931 (0.0 <= combination_coordinates[i]) &&
932 (combination_coordinates[i] <= 1.0);
933 in_quad = is_convex_combination;
936 in_quad = (0.0 <= c0 && c0 <= 1.0 && 0.0 <= c1 && c1 <= 1.0);
939 return std::make_pair(projected_shifted_p - shift,
942 return std::make_pair(
Point<3>(), std::numeric_limits<double>::max());
960 return project_to_line(vertex<dim>(0), vertex<dim>(1), p).first;
965 unsigned int closest_vertex_no = 0;
966 double closest_vertex_distance_square = vertex<dim>(0).
distance_square(p);
967 for (
unsigned int i = 1; i <
n_vertices(); ++i)
969 const double new_vertex_distance_square =
970 vertex<dim>(i).distance_square(p);
971 if (new_vertex_distance_square < closest_vertex_distance_square)
973 closest_vertex_no = i;
974 closest_vertex_distance_square = new_vertex_distance_square;
978 double min_distance_square = std::numeric_limits<double>::max();
981 for (
const unsigned int face_no :
987 auto pair = project_to_line(
v0,
v1, p);
988 if (pair.second < min_distance_square)
991 min_distance_square = pair.second;
1008 const std::array<unsigned int, 5> all_pyramid_faces{{0, 1, 2, 3, 4}};
1012 for (
const unsigned int face_no : faces)
1017 std::array<Point<dim>, 3>
vertices;
1018 for (
unsigned int vertex_no = 0; vertex_no < 3; ++vertex_no)
1022 auto pair = project_to_quad(
vertices, p, face_cell);
1023 if (pair.second < min_distance_square)
1025 result = pair.first;
1026 min_distance_square = pair.second;
1030 for (
const unsigned int face_no :
1034 for (
const unsigned int face_line_no : face_cell.line_indices())
1036 const auto cell_line_no =
1044 auto pair = project_to_line(
v0,
v1, p);
1045 if (pair.second < min_distance_square)
1047 result = pair.first;
1048 min_distance_square = pair.second;
1054 Assert(min_distance_square < std::numeric_limits<double>::max(),
1063 constexpr unsigned int x_index = 0;
1064 constexpr unsigned int y_index = (dim >= 2 ? 1 : 0);
1065 constexpr unsigned int z_index = (dim >= 3 ? 2 : 0);
1076 for (
unsigned int d = 0; d < dim; ++d)
1077 result[d] = std::clamp(result[d], 0.0, 1.0);
1081 result[x_index] = std::clamp(result[x_index], 0.0, 1.0);
1083 std::clamp(result[y_index], 0.0, 1.0 - result[x_index]);
1086 result[x_index] = std::clamp(result[x_index], 0.0, 1.0);
1088 std::clamp(result[y_index], 0.0, 1.0 - result[x_index]);
1090 std::clamp(result[z_index],
1092 1.0 - result[x_index] - result[y_index]);
1096 result[x_index] = std::clamp(result[x_index], 0.0, 1.0);
1098 std::clamp(result[y_index], 0.0, 1.0 - result[x_index]);
1099 result[z_index] = std::clamp(result[z_index], 0.0, 1.0);
1103 result[x_index] = std::clamp(result[x_index], -1.0, 1.0);
1104 result[y_index] = std::clamp(result[y_index], -1.0, 1.0);
1107 const auto x_abs =
std::abs(result[x_index]);
1108 const auto y_abs =
std::abs(result[y_index]);
1111 result[z_index] = std::clamp(result[z_index], 0.0, 1.0 - x_abs);
1113 result[z_index] = std::clamp(result[z_index], 0.0, 1.0 - y_abs);
1138 out << static_cast<unsigned int>(reference_cell.kind);
1152 reference_cell.kind =
static_cast<decltype(reference_cell.kind)
>(value);
1166 "The reference cell kind just read does not correspond to one of the "
1167 "valid choices. There must be an error."));
1179#include "reference_cell.inst"
Abstract base class for mapping classes.
constexpr numbers::NumberTraits< Number >::real_type distance_square(const Point< dim, Number > &p) const
ArrayView< const unsigned int > faces_for_given_vertex(const unsigned int vertex_index) const
unsigned int n_vertices() const
Quadrature< dim > get_gauss_type_quadrature(const unsigned n_points_1d) const
bool is_hyper_cube() const
unsigned int vtk_linear_type() const
unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const unsigned char face_orientation) const
static constexpr unsigned char default_combined_face_orientation()
const Quadrature< dim > & get_nodal_type_quadrature() const
static constexpr ndarray< unsigned int, 2, 2 > line_vertex_permutations
bool contains_point(const Point< dim > &p, const double tolerance=0) const
unsigned int gmsh_element_type() const
static constexpr ndarray< unsigned int, 8, 4 > quadrilateral_vertex_permutations
unsigned int n_faces() const
unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex) 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
unsigned int vtk_vertex_to_deal_vertex(const unsigned int vertex_index) 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
ReferenceCell face_reference_cell(const unsigned int face_no) 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
unsigned int face_to_cell_lines(const unsigned int face, const unsigned int line, const unsigned char face_orientation) const
Point< dim > closest_point(const Point< dim > &p) const
static constexpr ndarray< unsigned int, 6, 3 > triangle_vertex_permutations
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
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
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
std::istream & operator>>(std::istream &in, ReferenceCell &reference_cell)
std::ostream & operator<<(std::ostream &out, const ReferenceCell &reference_cell)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)