reference_cell.cc

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2020 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#include <deal.II/base/polynomial.h>
#include <deal.II/base/polynomials_barycentric.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/tensor_product_polynomials.h>
 
#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_simplex_p.h>
#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/fe/mapping_p1.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/fe/mapping_q1.h>
 
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/reference_cell.h>
#include <deal.II/grid/tria.h>
 
#include <algorithm>
#include <iostream>
#include <memory>
 
DEAL_II_NAMESPACE_OPEN
 
namespace
{
  namespace VTKCellType
  {
    // Define VTK constants for linear, quadratic and
    // high-order Lagrange geometrices
    enum : unsigned int
    {
      VTK_VERTEX = 1,
      // Linear cells
      VTK_LINE       = 3,
      VTK_TRIANGLE   = 5,
      VTK_QUAD       = 9,
      VTK_TETRA      = 10,
      VTK_HEXAHEDRON = 12,
      VTK_WEDGE      = 13,
      VTK_PYRAMID    = 14,
      // Quadratic cells
      VTK_QUADRATIC_EDGE       = 21,
      VTK_QUADRATIC_TRIANGLE   = 22,
      VTK_QUADRATIC_QUAD       = 23,
      VTK_QUADRATIC_TETRA      = 24,
      VTK_QUADRATIC_HEXAHEDRON = 25,
      VTK_QUADRATIC_WEDGE      = 26,
      VTK_QUADRATIC_PYRAMID    = 27,
      // Lagrange cells
      VTK_LAGRANGE_CURVE         = 68,
      VTK_LAGRANGE_TRIANGLE      = 69,
      VTK_LAGRANGE_QUADRILATERAL = 70,
      VTK_LAGRANGE_TETRAHEDRON   = 71,
      VTK_LAGRANGE_HEXAHEDRON    = 72,
      VTK_LAGRANGE_WEDGE         = 73,
      VTK_LAGRANGE_PYRAMID       = 74,
      // Invalid code
      VTK_INVALID = numbers::invalid_unsigned_int
    };
 
  } // namespace VTKCellType
 
} // namespace
 
 
 
template <int dim>
std::string
ReferenceCell<dim>::to_string() const
{
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        return "Vertex";
      case ReferenceCells::Line:
        return "Line";
      case ReferenceCells::Triangle:
        return "Tri";
      case ReferenceCells::Quadrilateral:
        return "Quad";
      case ReferenceCells::Tetrahedron:
        return "Tet";
      case ReferenceCells::Pyramid:
        return "Pyramid";
      case ReferenceCells::Wedge:
        return "Wedge";
      case ReferenceCells::Hexahedron:
        return "Hex";
      case ReferenceCells::Invalid<dim>:
        return "Invalid";
      default:
        DEAL_II_NOT_IMPLEMENTED();
    }
 
  return "Invalid";
}
 
 
 
namespace
{
  // Return vertex_no-th vertex of the child_no-ith child cell.
  //
  // @note This function is not yet implemented in 3d since it is only called on
  // faces of cells.
  template <int dim>
  Point<dim>
  child_vertex(const ReferenceCell<dim>  reference_cell,
               const unsigned int        child_no,
               const unsigned int        vertex_no,
               const RefinementCase<dim> refinement_case)
  {
    if (dim > 1)
      AssertIndexRange(child_no, reference_cell.n_children(refinement_case));
    AssertIndexRange(vertex_no, reference_cell.n_vertices());
 
    constexpr Point<dim> V0;
    // V isn't used for dim == 0
    [[maybe_unused]] const auto V = [](const unsigned int d) {
      return Point<dim>::unit_vector(d);
    };
 
    switch (reference_cell)
      {
        case ReferenceCells::Vertex:
          return V0;
        case ReferenceCells::Line:
          {
            Assert(refinement_case == RefinementCase<dim>::isotropic_refinement,
                   ExcNotImplemented());
            if constexpr (dim == 1)
              {
                static constexpr ndarray<Point<1>, 2, 2>
                  isotropic_child_vertices = {{
                    {{V0, 0.5 * V(0)}},
                    {{0.5 * V(0), V(0)}},
                  }};
                return isotropic_child_vertices[child_no][vertex_no];
              }
            else
              DEAL_II_ASSERT_UNREACHABLE();
          }
        case ReferenceCells::Triangle:
          {
            Assert(refinement_case == RefinementCase<dim>::isotropic_refinement,
                   ExcNotImplemented());
            if constexpr (dim == 2)
              {
                static constexpr ndarray<Point<2>, 4, 3>
                  isotropic_child_vertices = {{
                    {{V0, 0.5 * V(0), 0.5 * V(1)}},
                    {{0.5 * V(0), V(0), 0.5 * (V(0) + V(1))}},
                    {{0.5 * V(1), 0.5 * V(0) + 0.5 * V(1), V(1)}},
                    {{0.5 * V(0), 0.5 * V(0) + 0.5 * V(1), 0.5 * V(1)}},
                  }};
                return isotropic_child_vertices[child_no][vertex_no];
              }
            else
              DEAL_II_ASSERT_UNREACHABLE();
          }
        case ReferenceCells::Quadrilateral:
          {
            if constexpr (dim == 2)
              {
                static constexpr Point<2> M = 0.5 * (V(0) + V(1));
 
                static constexpr ndarray<Point<2>, 2, 4> cut_x_child_vertices =
                  {{
                    {{V0, 0.5 * V(0), V(1), V(1) + 0.5 * V(0)}},
                    {{0.5 * V(0), V(0), 0.5 * V(0) + V(1), V(0) + V(1)}},
                  }};
 
                static constexpr ndarray<Point<2>, 2, 4> cut_y_child_vertices =
                  {{
                    {{V0, V(0), 0.5 * V(1), V(0) + 0.5 * V(1)}},
                    {{0.5 * V(1), V(0) + 0.5 * V(1), V(1), V(0) + V(1)}},
                  }};
 
                static constexpr ndarray<Point<2>, 4, 4>
                  isotropic_child_vertices = {{
                    {{V0, 0.5 * V(0), 0.5 * V(1), M}},
                    {{0.5 * V(0), V(0), M, V(0) + 0.5 * V(1)}},
                    {{0.5 * V(1), M, V(1), V(1) + 0.5 * V(0)}},
                    {{M, V(0) + 0.5 * V(1), V(1) + 0.5 * V(0), V(0) + V(1)}},
                  }};
 
                switch (refinement_case)
                  {
                    case RefinementCase<2>::cut_x:
                      return cut_x_child_vertices[child_no][vertex_no];
                    case RefinementCase<2>::cut_y:
                      return cut_y_child_vertices[child_no][vertex_no];
                    case RefinementCase<2>::isotropic_refinement:
                      return isotropic_child_vertices[child_no][vertex_no];
                    default:
                      DEAL_II_ASSERT_UNREACHABLE();
                  }
              }
            else
              DEAL_II_ASSERT_UNREACHABLE();
          }
        case ReferenceCells::Tetrahedron:
        case ReferenceCells::Pyramid:
        case ReferenceCells::Wedge:
        case ReferenceCells::Hexahedron:
          DEAL_II_NOT_IMPLEMENTED();
          return {};
        default:
          DEAL_II_ASSERT_UNREACHABLE();
          return {};
      }
  }
} // namespace
 
 
 
template <int dim>
Point<dim>
ReferenceCell<dim>::subface_vertex_location(
  const unsigned int            face_no,
  const unsigned int            subface_no,
  const unsigned int            subface_vertex_no,
  const RefinementCase<dim - 1> face_refinement_case) const
{
  AssertIndexRange(face_no, n_faces());
  if (dim > 1)
    {
      AssertIndexRange(subface_no,
                       face_reference_cell(face_no).n_children(
                         face_refinement_case));
      Assert(face_refinement_case != RefinementCase<dim - 1>::no_refinement,
             ExcMessage("This function may only be called for subfaces."));
    }
  AssertIndexRange(subface_vertex_no,
                   face_reference_cell(face_no).n_vertices());
 
  Point<dim> p;
  for (const unsigned int vertex_no :
       face_reference_cell(face_no).vertex_indices())
    p += face_vertex_location(face_no, vertex_no) *
         face_reference_cell(face_no).d_linear_shape_function(
           child_vertex(face_reference_cell(face_no),
                        subface_no,
                        subface_vertex_no,
                        face_refinement_case),
           vertex_no);
 
  return p;
}
 
 
 
template <int dim>
std::pair<unsigned int, RefinementCase<dim - 1>>
ReferenceCell<dim>::equivalent_refinement_case(
  const types::geometric_orientation combined_face_orientation,
  const internal::SubfaceCase<dim>   subface_case,
  const unsigned int                 subface_no) const
{
  // 1. in 1d subfaces don't exist, but we still support some subface code
  //    (such as QProjector's functions) to enable dimension-independent
  //    programming. To match the convention used by 3d this will always return
  //    (0, isotropic_refinement).
  //
  // 2. historically we have permitted subface calculations in 2d for unrefined
  //    faces. In that case we don't actually need the value of subface_case
  //    since there is only one possible refinement - the returned value has to
  //    be isotropic_refinement.
  //
  // 3. Similarly, in 3d we treat case_none as case_isotropic.
  if constexpr (dim == 1)
    {
      (void)subface_case;
      AssertIndexRange(combined_face_orientation, n_face_orientations(0));
      AssertIndexRange(subface_no, 1);
      return std::make_pair(0, RefinementCase<dim - 1>::isotropic_refinement);
    }
 
  if constexpr (dim == 2)
    {
      (void)subface_case;
      AssertIndexRange(combined_face_orientation, n_face_orientations(0));
      AssertIndexRange(subface_no,
                       face_reference_cell(0).n_isotropic_children());
      return std::make_pair(combined_face_orientation ==
                                numbers::reverse_line_orientation ?
                              1 - subface_no :
                              subface_no,
                            RefinementCase<dim - 1>::isotropic_refinement);
    }
  else if constexpr (dim == 3)
    {
      Assert(*this == ReferenceCells::Hexahedron, ExcNotImplemented());
 
      static const RefinementCase<dim - 1>
        equivalent_refine_case[internal::SubfaceCase<dim>::case_isotropic + 1]
                              [GeometryInfo<3>::max_children_per_face] = {
                                // case_none. there should be only
                                // invalid values here. However, as
                                // this function is also called (in
                                // tests) for cells which have no
                                // refined faces, use isotropic
                                // refinement instead
                                {RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_xy},
                                // case_x
                                {RefinementCase<dim - 1>::cut_x,
                                 RefinementCase<dim - 1>::cut_x,
                                 RefinementCase<dim - 1>::no_refinement,
                                 RefinementCase<dim - 1>::no_refinement},
                                // case_x1y
                                {RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_x,
                                 RefinementCase<dim - 1>::no_refinement},
                                // case_x2y
                                {RefinementCase<dim - 1>::cut_x,
                                 RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::no_refinement},
                                // case_x1y2y
                                {RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_xy},
                                // case_y
                                {RefinementCase<dim - 1>::cut_y,
                                 RefinementCase<dim - 1>::cut_y,
                                 RefinementCase<dim - 1>::no_refinement,
                                 RefinementCase<dim - 1>::no_refinement},
                                // case_y1x
                                {RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_y,
                                 RefinementCase<dim - 1>::no_refinement},
                                // case_y2x
                                {RefinementCase<dim - 1>::cut_y,
                                 RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::no_refinement},
                                // case_y1x2x
                                {RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_xy},
                                // case_xy (case_isotropic)
                                {RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_xy,
                                 RefinementCase<dim - 1>::cut_xy}};
 
      constexpr unsigned int X = numbers::invalid_unsigned_int;
      static const unsigned int
        equivalent_subface_number[internal::SubfaceCase<dim>::case_isotropic +
                                  1][GeometryInfo<3>::max_children_per_face] = {
          // case_none, see above
          {0, 1, 2, 3},
          // case_x
          {0, 1, X, X},
          // case_x1y
          {0, 2, 1, X},
          // case_x2y
          {0, 1, 3, X},
          // case_x1y2y
          {0, 2, 1, 3},
          // case_y
          {0, 1, X, X},
          // case_y1x
          {0, 1, 1, X},
          // case_y2x
          {0, 2, 3, X},
          // case_y1x2x
          {0, 1, 2, 3},
          // case_xy (case_isotropic)
          {0, 1, 2, 3}};
 
      static const RefinementCase<dim - 1> rotated_refinement_case[4] = {
        RefinementCase<dim - 1>::no_refinement,
        RefinementCase<dim - 1>::cut_y,
        RefinementCase<dim - 1>::cut_x,
        RefinementCase<dim - 1>::cut_xy};
      const auto [face_orientation, face_rotation, face_flip] =
        internal::split_face_orientation(combined_face_orientation);
 
      const auto equivalent_refinement_case =
        equivalent_refine_case[subface_case][subface_no];
      const unsigned int equivalent_subface_no =
        equivalent_subface_number[subface_case][subface_no];
      // make sure, that we got a valid subface and RefineCase
      Assert(equivalent_refinement_case != RefinementCase<dim>::no_refinement,
             ExcInternalError());
      Assert(equivalent_subface_no != X, ExcInternalError());
      // now, finally respect non-standard faces
      const RefinementCase<dim - 1> final_refinement_case =
        (face_orientation == face_rotation ?
           rotated_refinement_case[equivalent_refinement_case] :
           equivalent_refinement_case);
 
      const unsigned int final_subface_no =
        GeometryInfo<dim>::child_cell_on_face(RefinementCase<dim>(
                                                final_refinement_case),
                                              /*face_no = */ 4,
                                              equivalent_subface_no,
                                              face_orientation,
                                              face_flip,
                                              face_rotation,
                                              equivalent_refinement_case);
 
      return std::make_pair(final_subface_no, final_refinement_case);
    }
  else
    {
      (void)combined_face_orientation;
      (void)subface_case;
      (void)subface_no;
 
      DEAL_II_NOT_IMPLEMENTED();
      return {};
    }
}
 
 
 
template <int dim>
template <int spacedim>
std::unique_ptr<Mapping<dim, spacedim>>
ReferenceCell<dim>::get_default_mapping(const unsigned int degree) const
{
  if (is_hyper_cube())
    return std::make_unique<MappingQ<dim, spacedim>>(degree);
  else if (is_simplex())
    if (degree == 1)
      return std::make_unique<MappingP1<dim, spacedim>>();
    else
      return std::make_unique<MappingFE<dim, spacedim>>(
        FE_SimplexP<dim, spacedim>(degree));
  else if (*this == ReferenceCells::Pyramid)
    return std::make_unique<MappingFE<dim, spacedim>>(
      FE_PyramidP<dim, spacedim>(degree));
  else if (*this == ReferenceCells::Wedge)
    return std::make_unique<MappingFE<dim, spacedim>>(
      FE_WedgeP<dim, spacedim>(degree));
  else
    {
      DEAL_II_NOT_IMPLEMENTED();
    }
 
  return std::make_unique<MappingQ<dim, spacedim>>(degree);
}
 
 
 
template <int dim>
template <int spacedim>
const Mapping<dim, spacedim> &
ReferenceCell<dim>::get_default_linear_mapping() const
{
  if (is_hyper_cube())
    {
      return StaticMappingQ1<dim, spacedim>::mapping;
    }
  else if (is_simplex())
    {
      static const MappingP1<dim, spacedim> mapping;
      return mapping;
    }
  else if (*this == ReferenceCells::Pyramid)
    {
      static const MappingFE<dim, spacedim> mapping(
        FE_PyramidP<dim, spacedim>(1));
      return mapping;
    }
  else if (*this == ReferenceCells::Wedge)
    {
      static const MappingFE<dim, spacedim> mapping(
        FE_WedgeP<dim, spacedim>(1));
      return mapping;
    }
  else
    {
      DEAL_II_NOT_IMPLEMENTED();
    }
 
  return StaticMappingQ1<dim, spacedim>::mapping; // never reached
}
 
 
 
template <int dim>
Quadrature<dim>
ReferenceCell<dim>::get_gauss_type_quadrature(const unsigned n_points_1d) const
{
  if (is_hyper_cube())
    return QGauss<dim>(n_points_1d);
  else if (is_simplex())
    return QGaussSimplex<dim>(n_points_1d);
  else if (*this == ReferenceCells::Pyramid)
    return QGaussPyramid<dim>(n_points_1d);
  else if (*this == ReferenceCells::Wedge)
    return QGaussWedge<dim>(n_points_1d);
  else
    DEAL_II_NOT_IMPLEMENTED();
 
  return Quadrature<dim>(); // never reached
}
 
 
 
template <int dim>
const Quadrature<dim> &
ReferenceCell<dim>::get_nodal_type_quadrature() const
{
  // A function that is used to fill a quadrature object of the
  // desired type the first time we encounter a particular
  // reference cell
  const auto create_quadrature = [](const ReferenceCell<dim> &reference_cell) {
    std::vector<Point<dim>> vertices(reference_cell.n_vertices());
    for (const unsigned int v : reference_cell.vertex_indices())
      vertices[v] = reference_cell.vertex(v);
 
    return Quadrature<dim>(vertices);
  };
 
  if (is_hyper_cube())
    {
      static const Quadrature<dim> quadrature = create_quadrature(*this);
      return quadrature;
    }
  else if (is_simplex())
    {
      static const Quadrature<dim> quadrature = create_quadrature(*this);
      return quadrature;
    }
  else if (*this == ReferenceCells::Pyramid)
    {
      static const Quadrature<dim> quadrature = create_quadrature(*this);
      return quadrature;
    }
  else if (*this == ReferenceCells::Wedge)
    {
      static const Quadrature<dim> quadrature = create_quadrature(*this);
      return quadrature;
    }
  else
    DEAL_II_NOT_IMPLEMENTED();
 
  static const Quadrature<dim> dummy;
  return dummy; // never reached
}
 
 
 
template <int dim>
unsigned int
ReferenceCell<dim>::exodusii_vertex_to_deal_vertex(
  const unsigned int vertex_n) const
{
  AssertIndexRange(vertex_n, n_vertices());
 
  switch (this->kind)
    {
      case ReferenceCells::Line:
      case ReferenceCells::Triangle:
        return vertex_n;
      case ReferenceCells::Quadrilateral:
        {
          constexpr std::array<unsigned int, 4> exodus_to_deal{{0, 1, 3, 2}};
          return exodus_to_deal[vertex_n];
        }
      case ReferenceCells::Tetrahedron:
        return vertex_n;
      case ReferenceCells::Hexahedron:
        {
          constexpr std::array<unsigned int, 8> exodus_to_deal{
            {0, 1, 3, 2, 4, 5, 7, 6}};
          return exodus_to_deal[vertex_n];
        }
      case ReferenceCells::Wedge:
        {
          constexpr std::array<unsigned int, 6> exodus_to_deal{
            {2, 1, 0, 5, 4, 3}};
          return exodus_to_deal[vertex_n];
        }
      case ReferenceCells::Pyramid:
        {
          constexpr std::array<unsigned int, 5> exodus_to_deal{{0, 1, 3, 2, 4}};
          return exodus_to_deal[vertex_n];
        }
      default:
        DEAL_II_NOT_IMPLEMENTED();
    }
 
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int dim>
unsigned int
ReferenceCell<dim>::exodusii_face_to_deal_face(const unsigned int face_n) const
{
  AssertIndexRange(face_n, n_faces());
 
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        return 0;
      case ReferenceCells::Line:
      case ReferenceCells::Triangle:
        return face_n;
      case ReferenceCells::Quadrilateral:
        {
          constexpr std::array<unsigned int, 4> exodus_to_deal{{2, 1, 3, 0}};
          return exodus_to_deal[face_n];
        }
      case ReferenceCells::Tetrahedron:
        {
          constexpr std::array<unsigned int, 4> exodus_to_deal{{1, 3, 2, 0}};
          return exodus_to_deal[face_n];
        }
      case ReferenceCells::Hexahedron:
        {
          constexpr std::array<unsigned int, 6> exodus_to_deal{
            {2, 1, 3, 0, 4, 5}};
          return exodus_to_deal[face_n];
        }
      case ReferenceCells::Wedge:
        {
          constexpr std::array<unsigned int, 6> exodus_to_deal{{3, 4, 2, 0, 1}};
          return exodus_to_deal[face_n];
        }
      case ReferenceCells::Pyramid:
        {
          constexpr std::array<unsigned int, 5> exodus_to_deal{{3, 2, 4, 1, 0}};
          return exodus_to_deal[face_n];
        }
      default:
        DEAL_II_NOT_IMPLEMENTED();
    }
 
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int dim>
unsigned int
ReferenceCell<dim>::unv_vertex_to_deal_vertex(const unsigned int vertex_n) const
{
  AssertIndexRange(vertex_n, n_vertices());
  // Information on this file format isn't easy to find - the documents here
  //
  // https://www.ceas3.uc.edu/sdrluff/
  //
  // Don't actually explain anything about the sections we care about (2412) in
  // any detail. For node numbering I worked backwards from what is actually in
  // our test files (since that's supposed to work), which all use some
  // non-standard clockwise numbering scheme which starts at the bottom right
  // vertex.
  if (*this == ReferenceCells::Line)
    {
      return vertex_n;
    }
  else if (*this == ReferenceCells::Quadrilateral)
    {
      constexpr std::array<unsigned int, 4> unv_to_deal{{1, 0, 2, 3}};
      return unv_to_deal[vertex_n];
    }
  else if (*this == ReferenceCells::Hexahedron)
    {
      constexpr std::array<unsigned int, 8> unv_to_deal{
        {6, 7, 5, 4, 2, 3, 1, 0}};
      return unv_to_deal[vertex_n];
    }
 
  DEAL_II_NOT_IMPLEMENTED();
 
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int dim>
unsigned int
ReferenceCell<dim>::vtk_linear_type() const
{
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        return VTKCellType::VTK_VERTEX;
      case ReferenceCells::Line:
        return VTKCellType::VTK_LINE;
      case ReferenceCells::Triangle:
        return VTKCellType::VTK_TRIANGLE;
      case ReferenceCells::Quadrilateral:
        return VTKCellType::VTK_QUAD;
      case ReferenceCells::Tetrahedron:
        return VTKCellType::VTK_TETRA;
      case ReferenceCells::Pyramid:
        return VTKCellType::VTK_PYRAMID;
      case ReferenceCells::Wedge:
        return VTKCellType::VTK_WEDGE;
      case ReferenceCells::Hexahedron:
        return VTKCellType::VTK_HEXAHEDRON;
      case ReferenceCells::Invalid<dim>:
        return VTKCellType::VTK_INVALID;
      default:
        DEAL_II_NOT_IMPLEMENTED();
    }
 
  return VTKCellType::VTK_INVALID;
}
 
 
 
template <int dim>
unsigned int
ReferenceCell<dim>::vtk_quadratic_type() const
{
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        return VTKCellType::VTK_VERTEX;
      case ReferenceCells::Line:
        return VTKCellType::VTK_QUADRATIC_EDGE;
      case ReferenceCells::Triangle:
        return VTKCellType::VTK_QUADRATIC_TRIANGLE;
      case ReferenceCells::Quadrilateral:
        return VTKCellType::VTK_QUADRATIC_QUAD;
      case ReferenceCells::Tetrahedron:
        return VTKCellType::VTK_QUADRATIC_TETRA;
      case ReferenceCells::Pyramid:
        return VTKCellType::VTK_QUADRATIC_PYRAMID;
      case ReferenceCells::Wedge:
        return VTKCellType::VTK_QUADRATIC_WEDGE;
      case ReferenceCells::Hexahedron:
        return VTKCellType::VTK_QUADRATIC_HEXAHEDRON;
      case ReferenceCells::Invalid<dim>:
        return VTKCellType::VTK_INVALID;
      default:
        DEAL_II_NOT_IMPLEMENTED();
    }
 
  return VTKCellType::VTK_INVALID;
}
 
 
 
template <int dim>
unsigned int
ReferenceCell<dim>::vtk_lagrange_type() const
{
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        return VTKCellType::VTK_VERTEX;
      case ReferenceCells::Line:
        return VTKCellType::VTK_LAGRANGE_CURVE;
      case ReferenceCells::Triangle:
        return VTKCellType::VTK_LAGRANGE_TRIANGLE;
      case ReferenceCells::Quadrilateral:
        return VTKCellType::VTK_LAGRANGE_QUADRILATERAL;
      case ReferenceCells::Tetrahedron:
        return VTKCellType::VTK_LAGRANGE_TETRAHEDRON;
      case ReferenceCells::Pyramid:
        return VTKCellType::VTK_LAGRANGE_PYRAMID;
      case ReferenceCells::Wedge:
        return VTKCellType::VTK_LAGRANGE_WEDGE;
      case ReferenceCells::Hexahedron:
        return VTKCellType::VTK_LAGRANGE_HEXAHEDRON;
      case ReferenceCells::Invalid<dim>:
        return VTKCellType::VTK_INVALID;
      default:
        DEAL_II_NOT_IMPLEMENTED();
    }
 
  return VTKCellType::VTK_INVALID;
}
 
 
 
template <int dim>
unsigned int
ReferenceCell<dim>::vtk_lexicographic_to_node_index(
  const std::array<unsigned, dim> &node_indices,
  const std::array<unsigned, dim> &nodes_per_direction,
  const bool                       legacy_format) const
{
  (void)legacy_format;
 
  if constexpr (dim == 0)
    {
      DEAL_II_NOT_IMPLEMENTED();
      return 0;
    }
  else if constexpr (dim == 1)
    {
      const unsigned int i = node_indices[0];
 
      const bool ibdy = (i == 0 || i == nodes_per_direction[0]);
      // How many boundaries do we lie on at once?
      const int nbdy = (ibdy ? 1 : 0);
 
      if (nbdy == 1) // Vertex DOF
        { // ijk is a corner node. Return the proper index (somewhere in [0,7]):
          return i ? 1 : 0;
        }
 
      const int offset = 2;
      return (i - 1) + offset;
    }
  else if constexpr (dim == 2)
    {
      Assert(*this == ReferenceCells::Quadrilateral, ExcNotImplemented());
 
      const unsigned int i = node_indices[0];
      const unsigned int j = node_indices[1];
 
      const bool ibdy = (i == 0 || i == nodes_per_direction[0]);
      const bool jbdy = (j == 0 || j == nodes_per_direction[1]);
      // How many boundaries do we lie on at once?
      const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0);
 
      if (nbdy == 2) // Vertex DOF
        { // ijk is a corner node. Return the proper index (somewhere in [0,3]):
          return (i != 0u ? (j != 0u ? 2 : 1) : (j != 0u ? 3 : 0));
        }
 
      int offset = 4;
      if (nbdy == 1) // Edge DOF
        {
          if (!ibdy)
            { // On i axis
              return (i - 1) +
                     (j != 0u ? nodes_per_direction[0] - 1 +
                                  nodes_per_direction[1] - 1 :
                                0) +
                     offset;
            }
 
          if (!jbdy)
            { // On j axis
              return (j - 1) +
                     (i != 0u ? nodes_per_direction[0] - 1 :
                                2 * (nodes_per_direction[0] - 1) +
                                  nodes_per_direction[1] - 1) +
                     offset;
            }
        }
 
      offset += 2 * (nodes_per_direction[0] - 1 + nodes_per_direction[1] - 1);
      // nbdy == 0: Face DOF
      return offset + (i - 1) + (nodes_per_direction[0] - 1) * ((j - 1));
    }
  else if constexpr (dim == 3)
    {
      Assert(*this == ReferenceCells::Hexahedron, ExcNotImplemented());
 
      const unsigned int i = node_indices[0];
      const unsigned int j = node_indices[1];
      const unsigned int k = node_indices[2];
 
      const bool ibdy = (i == 0 || i == nodes_per_direction[0]);
      const bool jbdy = (j == 0 || j == nodes_per_direction[1]);
      const bool kbdy = (k == 0 || k == nodes_per_direction[2]);
      // How many boundaries do we lie on at once?
      const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0) + (kbdy ? 1 : 0);
 
      if (nbdy == 3) // Vertex DOF
        { // ijk is a corner node. Return the proper index (somewhere in [0,7]):
          return (i != 0u ? (j != 0u ? 2 : 1) : (j != 0u ? 3 : 0)) +
                 (k != 0u ? 4 : 0);
        }
 
      int offset = 8;
      if (nbdy == 2) // Edge DOF
        {
          if (!ibdy)
            { // On i axis
              return (i - 1) +
                     (j != 0u ? nodes_per_direction[0] - 1 +
                                  nodes_per_direction[1] - 1 :
                                0) +
                     (k != 0u ? 2 * (nodes_per_direction[0] - 1 +
                                     nodes_per_direction[1] - 1) :
                                0) +
                     offset;
            }
          if (!jbdy)
            { // On j axis
              return (j - 1) +
                     (i != 0u ? nodes_per_direction[0] - 1 :
                                2 * (nodes_per_direction[0] - 1) +
                                  nodes_per_direction[1] - 1) +
                     (k != 0u ? 2 * (nodes_per_direction[0] - 1 +
                                     nodes_per_direction[1] - 1) :
                                0) +
                     offset;
            }
          // !kbdy, On k axis
          offset +=
            4 * (nodes_per_direction[0] - 1) + 4 * (nodes_per_direction[1] - 1);
          if (legacy_format)
            return (k - 1) +
                   (nodes_per_direction[2] - 1) *
                     (i != 0u ? (j != 0u ? 3 : 1) : (j != 0u ? 2 : 0)) +
                   offset;
          else
            return (k - 1) +
                   (nodes_per_direction[2] - 1) *
                     (i != 0u ? (j != 0u ? 2 : 1) : (j != 0u ? 3 : 0)) +
                   offset;
        }
 
      offset += 4 * (nodes_per_direction[0] - 1 + nodes_per_direction[1] - 1 +
                     nodes_per_direction[2] - 1);
      if (nbdy == 1) // Face DOF
        {
          if (ibdy) // On i-normal face
            {
              return (j - 1) + ((nodes_per_direction[1] - 1) * (k - 1)) +
                     (i != 0u ? (nodes_per_direction[1] - 1) *
                                  (nodes_per_direction[2] - 1) :
                                0) +
                     offset;
            }
          offset +=
            2 * (nodes_per_direction[1] - 1) * (nodes_per_direction[2] - 1);
          if (jbdy) // On j-normal face
            {
              return (i - 1) + ((nodes_per_direction[0] - 1) * (k - 1)) +
                     (j != 0u ? (nodes_per_direction[2] - 1) *
                                  (nodes_per_direction[0] - 1) :
                                0) +
                     offset;
            }
          offset +=
            2 * (nodes_per_direction[2] - 1) * (nodes_per_direction[0] - 1);
          // kbdy, On k-normal face
          return (i - 1) + ((nodes_per_direction[0] - 1) * (j - 1)) +
                 (k != 0u ? (nodes_per_direction[0] - 1) *
                              (nodes_per_direction[1] - 1) :
                            0) +
                 offset;
        }
 
      // nbdy == 0: Body DOF
      offset +=
        2 * ((nodes_per_direction[1] - 1) * (nodes_per_direction[2] - 1) +
             (nodes_per_direction[2] - 1) * (nodes_per_direction[0] - 1) +
             (nodes_per_direction[0] - 1) * (nodes_per_direction[1] - 1));
      return offset + (i - 1) +
             (nodes_per_direction[0] - 1) *
               ((j - 1) + (nodes_per_direction[1] - 1) * ((k - 1)));
    }
  else
    {
      DEAL_II_NOT_IMPLEMENTED();
      return 0;
    }
}
 
 
template <int dim>
unsigned int
ReferenceCell<dim>::vtk_vertex_to_deal_vertex(
  const unsigned int vertex_index) const
{
  AssertIndexRange(vertex_index, n_vertices());
 
  // For some of the following, deal.II uses the same ordering as VTK
  // and in that case, we only need to return 'vertex_index' (i.e.,
  // use the identity mapping). For some others, we need to translate.
  //
  // For the ordering, see the VTK manual (for example at
  // http://www.princeton.edu/~efeibush/viscourse/vtk.pdf, page 9).
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        return vertex_index;
      case ReferenceCells::Line:
        return vertex_index;
      case ReferenceCells::Triangle:
        return vertex_index;
      case ReferenceCells::Quadrilateral:
        {
          static constexpr std::array<unsigned int, 4> index_translation_table =
            {{0, 1, 3, 2}};
          return index_translation_table[vertex_index];
        }
      case ReferenceCells::Tetrahedron:
        return vertex_index;
      case ReferenceCells::Pyramid:
        {
          static constexpr std::array<unsigned int, 5> index_translation_table =
            {{0, 1, 3, 2, 4}};
          return index_translation_table[vertex_index];
        }
      case ReferenceCells::Wedge:
        return vertex_index;
      case ReferenceCells::Hexahedron:
        {
          static constexpr std::array<unsigned int, 8> index_translation_table =
            {{0, 1, 3, 2, 4, 5, 7, 6}};
          return index_translation_table[vertex_index];
        }
      case ReferenceCells::Invalid<dim>:
        {
          DEAL_II_NOT_IMPLEMENTED();
          return numbers::invalid_unsigned_int;
        }
      default:
        DEAL_II_NOT_IMPLEMENTED();
    }
 
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int dim>
unsigned int
ReferenceCell<dim>::gmsh_element_type() const
{
  /*
    From the GMSH documentation:
 
    elm-type
    defines the geometrical type of the n-th element:
 
    1
    Line (2 nodes).
 
    2
    Triangle (3 nodes).
 
    3
    Quadrangle (4 nodes).
 
    4
    Tetrahedron (4 nodes).
 
    5
    Hexahedron (8 nodes).
 
    6
    Prism (6 nodes).
 
    7
    Pyramid (5 nodes).
 
    8
    Second order line (3 nodes: 2 associated with the vertices and 1 with the
    edge).
 
    9
    Second order triangle (6 nodes: 3 associated with the vertices and 3 with
    the edges).
 
    10 Second order quadrangle (9 nodes: 4 associated with the
    vertices, 4 with the edges and 1 with the face).
 
    11 Second order tetrahedron (10 nodes: 4 associated with the vertices and 6
    with the edges).
 
    12 Second order hexahedron (27 nodes: 8 associated with the vertices, 12
    with the edges, 6 with the faces and 1 with the volume).
 
    13 Second order prism (18 nodes: 6 associated with the vertices, 9 with the
    edges and 3 with the quadrangular faces).
 
    14 Second order pyramid (14 nodes: 5 associated with the vertices, 8 with
    the edges and 1 with the quadrangular face).
 
    15 Point (1 node).
  */
 
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        return 15;
      case ReferenceCells::Line:
        return 1;
      case ReferenceCells::Triangle:
        return 2;
      case ReferenceCells::Quadrilateral:
        return 3;
      case ReferenceCells::Tetrahedron:
        return 4;
      case ReferenceCells::Pyramid:
        return 7;
      case ReferenceCells::Wedge:
        return 6;
      case ReferenceCells::Hexahedron:
        return 5;
      case ReferenceCells::Invalid<dim>:
      default:
        DEAL_II_NOT_IMPLEMENTED();
    }
 
  return numbers::invalid_unsigned_int;
}
 
 
 
namespace
{
  // Compute the nearest point to @p on the line segment and the square of its
  // distance to @p.
  template <int dim>
  std::pair<Point<dim>, double>
  project_to_line(const Point<dim> &x0,
                  const Point<dim> &x1,
                  const Point<dim> &p)
  {
    Assert(x0 != x1, ExcInternalError());
    // t is the convex combination coefficient (x = (1 - t) * x0 + t * x1)
    // defining the position of the closest point on the line (not line segment)
    // to p passing through x0 and x1. This formula is equivalent to the
    // standard 'project a vector onto another vector', where each vector is
    // shifted to start at x0.
    const double t = ((x1 - x0) * (p - x0)) / ((x1 - x0).norm_square());
 
    // Only consider points between x0 and x1
    if (t <= 0)
      return std::make_pair(x0, x0.distance_square(p));
    else if (t <= 1)
      {
        const auto p2 = x0 + t * (x1 - x0);
        return std::make_pair(p2, p2.distance_square(p));
      }
    else
      return std::make_pair(x1, x1.distance_square(p));
  }
 
 
 
  // template base-case
  template <int dim>
  std::pair<Point<dim>, double>
  project_to_quad(const std::array<Point<dim>, 3> & /*vertices*/,
                  const Point<dim> & /*p*/,
                  const ReferenceCell<dim - 1> /*face_reference_cell*/)
  {
    DEAL_II_ASSERT_UNREACHABLE();
    return std::make_pair(Point<dim>(),
                          std::numeric_limits<double>::signaling_NaN());
  }

  template <>
  std::pair<Point<3>, double>
  project_to_quad(const std::array<Point<3>, 3> &vertices,
                  const Point<3>                &p,
                  const ReferenceCell<2>         face_reference_cell)
  {
    Assert(face_reference_cell == ReferenceCells::Triangle ||
             face_reference_cell == ReferenceCells::Quadrilateral,
           ExcNotImplemented());
 
    // Make the problem slightly easier by shifting everything to avoid a point
    // at the origin (this way we can invert the matrix of vertices). Use 2.0 so
    // that the bottom left vertex of a Pyramid is now at x = 1.
    std::array<Point<3>, 3> shifted_vertices = vertices;
    const Tensor<1, 3>      shift{{2.0, 2.0, 2.0}};
    for (Point<3> &shifted_vertex : shifted_vertices)
      shifted_vertex += shift;
    const Point<3> shifted_p = p + shift;
 
    // As we are projecting onto a face of a reference cell, the vectors
    // describing its local coordinate system should be orthogonal. We don't
    // know which of the three vectors computed from p are mutually orthogonal
    // for triangles so that case requires an extra check.
    Tensor<1, 3>   e0;
    Tensor<1, 3>   e1;
    const Point<3> vertex = shifted_vertices[0];
    // Triangles are difficult because of two cases:
    // 1. the top face of a Tetrahedron, which does not have a right angle
    // 2. wedges and pyramids, whose faces do not lie on the reference simplex
    //
    // Deal with both by creating a locally orthogonal (but not necessarily
    // orthonormal) coordinate system and testing if the projected point is in
    // the triangle by expressing it as a convex combination of the vertices.
    if (face_reference_cell == ReferenceCells::Triangle)
      {
        e0 = shifted_vertices[1] - shifted_vertices[0];
        e1 = shifted_vertices[2] - shifted_vertices[0];
        e1 -= (e0 * e1) * e0 / (e0.norm_square());
      }
    else
      {
        e0 = shifted_vertices[1] - shifted_vertices[0];
        e1 = shifted_vertices[2] - shifted_vertices[0];
      }
    Assert(std::abs(e0 * e1) <= 1e-14, ExcInternalError());
    // the quadrilaterals on pyramids and wedges don't necessarily have edge
    // lengths of 1 so we cannot skip the denominator
    const double   c0 = e0 * (shifted_p - vertex) / e0.norm_square();
    const double   c1 = e1 * (shifted_p - vertex) / e1.norm_square();
    const Point<3> projected_shifted_p = vertex + c0 * e0 + c1 * e1;
 
    bool in_quad = false;
    if (face_reference_cell == ReferenceCells::Triangle)
      {
        Tensor<2, 3> shifted_vertex_matrix;
        for (unsigned int i = 0; i < 3; ++i)
          shifted_vertex_matrix[i] = shifted_vertices[i];
        const Tensor<1, 3> combination_coordinates =
          invert(transpose(shifted_vertex_matrix)) * projected_shifted_p;
        bool is_convex_combination = true;
        for (unsigned int i = 0; i < 3; ++i)
          is_convex_combination = is_convex_combination &&
                                  (0.0 <= combination_coordinates[i]) &&
                                  (combination_coordinates[i] <= 1.0);
        in_quad = is_convex_combination;
      }
    else
      in_quad = (0.0 <= c0 && c0 <= 1.0 && 0.0 <= c1 && c1 <= 1.0);
 
    if (in_quad)
      return std::make_pair(projected_shifted_p - shift,
                            shifted_p.distance_square(projected_shifted_p));
    else
      return std::make_pair(Point<3>(), std::numeric_limits<double>::max());
  }
} // namespace
 
 
 
template <int dim>
Point<dim>
ReferenceCell<dim>::closest_point(const Point<dim> &p) const
{
  // Handle simple cases first:
  if constexpr (dim == 0)
    return p;
  if (contains_point(p, 0.0))
    return p;
 
  if constexpr (dim == 1)
    return project_to_line(vertex(0), vertex(1), p).first;
  else
    {
      // Find the closest vertex so that we only need to check adjacent faces
      // and lines.
      Point<dim>   result;
      unsigned int closest_vertex_no        = 0;
      double closest_vertex_distance_square = vertex(0).distance_square(p);
      for (unsigned int i = 1; i < n_vertices(); ++i)
        {
          const double new_vertex_distance_square =
            vertex(i).distance_square(p);
          if (new_vertex_distance_square < closest_vertex_distance_square)
            {
              closest_vertex_no              = i;
              closest_vertex_distance_square = new_vertex_distance_square;
            }
        }
 
      double min_distance_square = std::numeric_limits<double>::max();
      if constexpr (dim == 2)
        {
          for (const unsigned int face_no :
               faces_for_given_vertex(closest_vertex_no))
            {
              const Point<dim> v0 = vertex(line_to_cell_vertices(face_no, 0));
              const Point<dim> v1 = vertex(line_to_cell_vertices(face_no, 1));
 
              auto pair = project_to_line(v0, v1, p);
              if (pair.second < min_distance_square)
                {
                  result              = pair.first;
                  min_distance_square = pair.second;
                }
            }
        }
      else
        {
          // Check faces and then lines.
          //
          // For reference cells with sloped faces (i.e., all 3D shapes except
          // Hexahedra), we might be able to do a valid normal projection to a
          // face with a different slope which is on the 'other side' of the
          // reference cell. To catch that case we have to unconditionally check
          // lines after checking faces.
          //
          // For pyramids the closest vertex might not be on the closest face:
          // for example, the origin is closest to vertex 4 which is not on the
          // bottom plane. Get around that by just checking all faces for
          // pyramids.
          const std::array<unsigned int, 5> all_pyramid_faces{{0, 1, 2, 3, 4}};
          const auto                       &faces =
            *this == ReferenceCells::Pyramid ?
                                    ArrayView<const unsigned int>(all_pyramid_faces) :
                                    faces_for_given_vertex(closest_vertex_no);
          for (const unsigned int face_no : faces)
            {
              auto face_cell = face_reference_cell(face_no);
              // We only need the first three points since for quads the last
              // point is redundant
              std::array<Point<dim>, 3> vertices;
              for (unsigned int vertex_no = 0; vertex_no < 3; ++vertex_no)
                vertices[vertex_no] = vertex(face_to_cell_vertices(
                  face_no, vertex_no, numbers::default_geometric_orientation));
 
              auto pair = project_to_quad(vertices, p, face_cell);
              if (pair.second < min_distance_square)
                {
                  result              = pair.first;
                  min_distance_square = pair.second;
                }
            }
 
          for (const unsigned int face_no :
               faces_for_given_vertex(closest_vertex_no))
            {
              auto face_cell = face_reference_cell(face_no);
              for (const unsigned int face_line_no : face_cell.line_indices())
                {
                  const auto cell_line_no =
                    face_to_cell_lines(face_no,
                                       face_line_no,
                                       numbers::default_geometric_orientation);
                  const auto v0 =
                    vertex(line_to_cell_vertices(cell_line_no, 0));
                  const auto v1 =
                    vertex(line_to_cell_vertices(cell_line_no, 1));
                  auto pair = project_to_line(v0, v1, p);
                  if (pair.second < min_distance_square)
                    {
                      result              = pair.first;
                      min_distance_square = pair.second;
                    }
                }
            }
        }
 
      Assert(min_distance_square < std::numeric_limits<double>::max(),
             ExcInternalError());
 
      // If necessary, slightly adjust the computed point so that it is closer
      // to being on the surface of the reference cell. Due to roundoff it is
      // difficult to place points on sloped surfaces (e.g., for Pyramids) so
      // this check isn't perfect but does improve the accuracy of the projected
      // point.
      if (!contains_point(result, 0.0))
        {
          constexpr unsigned int x_index = 0;
          constexpr unsigned int y_index = (dim >= 2 ? 1 : 0);
          constexpr unsigned int z_index = (dim >= 3 ? 2 : 0);
          switch (this->kind)
            {
              case ReferenceCells::Vertex:
                DEAL_II_ASSERT_UNREACHABLE();
                break;
                // the bounds for each dimension of a hypercube are mutually
                // independent:
              case ReferenceCells::Line:
              case ReferenceCells::Quadrilateral:
              case ReferenceCells::Hexahedron:
                for (unsigned int d = 0; d < dim; ++d)
                  result[d] = std::clamp(result[d], 0.0, 1.0);
                // simplices can use the standard definition of a simplex:
                break;
              case ReferenceCells::Triangle:
                result[x_index] = std::clamp(result[x_index], 0.0, 1.0);
                result[y_index] =
                  std::clamp(result[y_index], 0.0, 1.0 - result[x_index]);
                break;
              case ReferenceCells::Tetrahedron:
                result[x_index] = std::clamp(result[x_index], 0.0, 1.0);
                result[y_index] =
                  std::clamp(result[y_index], 0.0, 1.0 - result[x_index]);
                result[z_index] =
                  std::clamp(result[z_index],
                             0.0,
                             1.0 - result[x_index] - result[y_index]);
                break;
              // wedges and pyramids are more ad-hoc:
              case ReferenceCells::Wedge:
                result[x_index] = std::clamp(result[x_index], 0.0, 1.0);
                result[y_index] =
                  std::clamp(result[y_index], 0.0, 1.0 - result[x_index]);
                result[z_index] = std::clamp(result[z_index], 0.0, 1.0);
                break;
              case ReferenceCells::Pyramid:
                {
                  result[x_index] = std::clamp(result[x_index], -1.0, 1.0);
                  result[y_index] = std::clamp(result[y_index], -1.0, 1.0);
                  // It suffices to transform everything to the first quadrant
                  // to adjust z:
                  const auto x_abs = std::abs(result[x_index]);
                  const auto y_abs = std::abs(result[y_index]);
 
                  if (y_abs <= x_abs)
                    result[z_index] =
                      std::clamp(result[z_index], 0.0, 1.0 - x_abs);
                  else
                    result[z_index] =
                      std::clamp(result[z_index], 0.0, 1.0 - y_abs);
                }
                break;
              default:
                DEAL_II_NOT_IMPLEMENTED();
            }
        }
      // We should be within 4 * eps of the cell by this point. The roundoff
      // error comes from, e.g., computing (1 - x) + x when moving points onto
      // the top of a Pyramid.
      Assert(contains_point(result,
                            4.0 * std::numeric_limits<double>::epsilon()),
             ExcInternalError());
 
      return result;
    }
}
 
 
 
template <int dim>
std::ostream &
operator<<(std::ostream &out, const ReferenceCell<dim> &reference_cell)
{
  AssertThrow(out.fail() == false, ExcIO());
 
  // Output as an integer to avoid outputting it as a character with
  // potentially non-printing value:
  out << static_cast<unsigned int>(reference_cell.kind);
  return out;
}
 
 
 
template <int dim>
std::istream &
operator>>(std::istream &in, ReferenceCell<dim> &reference_cell)
{
  AssertThrow(in.fail() == false, ExcIO());
 
  // Read the information as an integer and convert it to the correct type
  unsigned int value;
  in >> value;
  reference_cell.kind = static_cast<decltype(reference_cell.kind)>(value);
 
  // Ensure that the object we read is valid
  Assert(
    (reference_cell == ReferenceCells::Vertex) ||
      (reference_cell == ReferenceCells::Line) ||
      (reference_cell == ReferenceCells::Triangle) ||
      (reference_cell == ReferenceCells::Quadrilateral) ||
      (reference_cell == ReferenceCells::Tetrahedron) ||
      (reference_cell == ReferenceCells::Hexahedron) ||
      (reference_cell == ReferenceCells::Wedge) ||
      (reference_cell == ReferenceCells::Pyramid) ||
      (reference_cell == ReferenceCells::Invalid<dim>),
    ExcMessage(
      "The reference cell kind just read does not correspond to one of the "
      "valid choices. There must be an error."));
 
  return in;
}
 
 
#include "grid/reference_cell.inst"
 
#ifndef DOXYGEN
// Here also instantiate select members for dimension 0 (vertices) in addition
// to the general class instantiation for dim=1,2,3 in the .inst file :
template std::string
ReferenceCell<0>::to_string() const;
 
template unsigned int
ReferenceCell<0>::vtk_linear_type() const;
 
template unsigned int
ReferenceCell<0>::vtk_quadratic_type() const;
 
template unsigned int
ReferenceCell<0>::vtk_lagrange_type() const;
 
template unsigned int
ReferenceCell<0>::vtk_vertex_to_deal_vertex(
  const unsigned int vertex_index) const;
 
template Quadrature<0>
ReferenceCell<0>::get_gauss_type_quadrature(const unsigned n_points_1d) const;
 
template const Quadrature<0> &
ReferenceCell<0>::get_nodal_type_quadrature() const;
 
template std::ostream &
operator<<(std::ostream &out, const ReferenceCell<0> &reference_cell);
 
template std::istream &
operator>>(std::istream &in, ReferenceCell<0> &reference_cell);
#endif
 
DEAL_II_NAMESPACE_CLOSE