reference_cell.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2020 - 2022 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------
 
#ifndef dealii_tria_reference_cell_h
#define dealii_tria_reference_cell_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/array_view.h>
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/ndarray.h>
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/utilities.h>
 
#include <iosfwd>
#include <string>
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations
#ifndef DOXYGEN
template <int dim, int spacedim>
class Mapping;
 
template <int dim>
class Quadrature;
 
class ReferenceCell;
#endif
 
 
namespace internal
{
  namespace ReferenceCell
  {
    DEAL_II_CONSTEXPR ::ReferenceCell
                      make_reference_cell_from_int(const std::uint8_t kind);
 
  } // namespace ReferenceCell
} // namespace internal
 
 
 
class ReferenceCell
{
public:
  static ReferenceCell
  n_vertices_to_type(const int dim, const unsigned int n_vertices);
 
  DEAL_II_CONSTEXPR
  ReferenceCell();
 
  bool
  is_hyper_cube() const;
 
  bool
  is_simplex() const;
 
  unsigned int
  get_dimension() const;
 
  template <int dim>
  double
  d_linear_shape_function(const Point<dim> &xi, const unsigned int i) const;
 
  template <int dim>
  Tensor<1, dim>
  d_linear_shape_function_gradient(const Point<dim> & xi,
                                   const unsigned int i) const;
 
  template <int dim, int spacedim = dim>
  std::unique_ptr<Mapping<dim, spacedim>>
  get_default_mapping(const unsigned int degree) const;
 
  template <int dim, int spacedim = dim>
  const Mapping<dim, spacedim> &
  get_default_linear_mapping() const;
 
  template <int dim>
  Quadrature<dim>
  get_gauss_type_quadrature(const unsigned n_points_1D) const;
 
  template <int dim>
  Quadrature<dim>
  get_midpoint_quadrature() const;
 
  template <int dim>
  const Quadrature<dim> &
  get_nodal_type_quadrature() const;
 
  unsigned int
  n_vertices() const;
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  vertex_indices() const;
 
  template <int dim>
  Point<dim>
  vertex(const unsigned int v) const;
 
  unsigned int
  n_lines() const;
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  line_indices() const;
 
  unsigned int
  n_faces() const;
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  face_indices() const;
 
  unsigned int
  n_isotropic_children() const;
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  isotropic_child_indices() const;
 
  ReferenceCell
  face_reference_cell(const unsigned int face_no) const;
 
  unsigned int
  child_cell_on_face(const unsigned int  face,
                     const unsigned int  subface,
                     const unsigned char face_orientation = 1) const;
 
  std::array<unsigned int, 2>
  standard_vertex_to_face_and_vertex_index(const unsigned int vertex) const;
 
  std::array<unsigned int, 2>
  standard_line_to_face_and_line_index(const unsigned int line) const;
 
  unsigned int
  face_to_cell_lines(const unsigned int  face,
                     const unsigned int  line,
                     const unsigned char face_orientation) const;
 
  unsigned int
  face_to_cell_vertices(const unsigned int  face,
                        const unsigned int  vertex,
                        const unsigned char face_orientation) const;
 
  unsigned int
  standard_to_real_face_vertex(const unsigned int  vertex,
                               const unsigned int  face,
                               const unsigned char face_orientation) const;
 
  unsigned int
  standard_to_real_face_line(const unsigned int  line,
                             const unsigned int  face,
                             const unsigned char face_orientation) const;
 
  bool
  standard_vs_true_line_orientation(const unsigned int  line,
                                    const unsigned char face_orientation,
                                    const bool          line_orientation) const;
 
  double
  volume() const;
 
  template <int dim>
  Point<dim>
  barycenter() const;
 
  template <int dim>
  bool
  contains_point(const Point<dim> &p, const double tolerance = 0) const;
 
  /*
   * Return @f$i@f$-th unit tangential vector of a face of the reference cell.
   * The vectors are arranged such that the
   * cross product between the two vectors returns the unit normal vector.
   *
   * @pre @f$i@f$ must be between zero and `dim-1`.
   */
  template <int dim>
  Tensor<1, dim>
  unit_tangential_vectors(const unsigned int face_no,
                          const unsigned int i) const;
 
  template <int dim>
  Tensor<1, dim>
  unit_normal_vectors(const unsigned int face_no) const;
 
  template <typename T, std::size_t N>
  unsigned char
  compute_orientation(const std::array<T, N> &vertices_0,
                      const std::array<T, N> &vertices_1) const;
 
  template <typename T, std::size_t N>
  std::array<T, N>
  permute_according_orientation(const std::array<T, N> &vertices,
                                const unsigned int      orientation) const;
 
  ArrayView<const unsigned int>
  faces_for_given_vertex(const unsigned int vertex_index) const;
 
  unsigned int
  exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const;
 
  unsigned int
  exodusii_face_to_deal_face(const unsigned int face_n) const;
 
  unsigned int
  unv_vertex_to_deal_vertex(const unsigned int vertex_n) const;
 
  unsigned int
  vtk_linear_type() const;
 
  unsigned int
  vtk_quadratic_type() const;
 
  unsigned int
  vtk_lagrange_type() const;
 
  unsigned int
  gmsh_element_type() const;
 
  std::string
  to_string() const;
 
  constexpr operator std::uint8_t() const;
 
  constexpr bool
  operator==(const ReferenceCell &type) const;
 
  constexpr bool
  operator!=(const ReferenceCell &type) const;
 
  template <class Archive>
  void
  serialize(Archive &archive, const unsigned int /*version*/);
 
private:
  std::uint8_t kind;
 
  constexpr ReferenceCell(const std::uint8_t kind);
 
  friend DEAL_II_CONSTEXPR ReferenceCell
  internal::ReferenceCell::make_reference_cell_from_int(const std::uint8_t);
 
  friend std::ostream &
  operator<<(std::ostream &out, const ReferenceCell &reference_cell);
 
  friend std::istream &
  operator>>(std::istream &in, ReferenceCell &reference_cell);
};
 
 
std::ostream &
operator<<(std::ostream &out, const ReferenceCell &reference_cell);
 
std::istream &
operator>>(std::istream &in, ReferenceCell &reference_cell);
 
 
inline constexpr ReferenceCell::ReferenceCell(const std::uint8_t kind)
  : kind(kind)
{}
 
 
 
inline constexpr ReferenceCell::operator std::uint8_t() const
{
  return kind;
}
 
 
 
inline constexpr bool
ReferenceCell::operator==(const ReferenceCell &type) const
{
  return kind == type.kind;
}
 
 
 
inline constexpr bool
ReferenceCell::operator!=(const ReferenceCell &type) const
{
  return kind != type.kind;
}
 
 
 
namespace internal
{
  namespace ReferenceCell
  {
    inline DEAL_II_CONSTEXPR ::ReferenceCell
    make_reference_cell_from_int(const std::uint8_t kind)
    {
      // Make sure these are the only indices from which objects can be
      // created.
      Assert((kind == static_cast<std::uint8_t>(-1)) || (kind < 8),
             ExcInternalError());
 
      // Call the private constructor, which we can from here because this
      // function is a 'friend'.
      return {kind};
    }
  } // namespace ReferenceCell
} // namespace internal
 
 
 
namespace ReferenceCells
{
  DEAL_II_CONSTEXPR const ReferenceCell Vertex =
    internal::ReferenceCell::make_reference_cell_from_int(0);
  DEAL_II_CONSTEXPR const ReferenceCell Line =
    internal::ReferenceCell::make_reference_cell_from_int(1);
  DEAL_II_CONSTEXPR const ReferenceCell Triangle =
    internal::ReferenceCell::make_reference_cell_from_int(2);
  DEAL_II_CONSTEXPR const ReferenceCell Quadrilateral =
    internal::ReferenceCell::make_reference_cell_from_int(3);
  DEAL_II_CONSTEXPR const ReferenceCell Tetrahedron =
    internal::ReferenceCell::make_reference_cell_from_int(4);
  DEAL_II_CONSTEXPR const ReferenceCell Pyramid =
    internal::ReferenceCell::make_reference_cell_from_int(5);
  DEAL_II_CONSTEXPR const ReferenceCell Wedge =
    internal::ReferenceCell::make_reference_cell_from_int(6);
  DEAL_II_CONSTEXPR const ReferenceCell Hexahedron =
    internal::ReferenceCell::make_reference_cell_from_int(7);
  DEAL_II_CONSTEXPR const ReferenceCell Invalid =
    internal::ReferenceCell::make_reference_cell_from_int(
      static_cast<std::uint8_t>(-1));
 
  template <int dim>
  constexpr const ReferenceCell &
  get_simplex();
 
  template <int dim>
  constexpr const ReferenceCell &
  get_hypercube();
} // namespace ReferenceCells
 
 
 
inline DEAL_II_CONSTEXPR
ReferenceCell::ReferenceCell()
  : ReferenceCell(ReferenceCells::Invalid)
{}
 
 
 
template <class Archive>
inline void
ReferenceCell::serialize(Archive &archive, const unsigned int /*version*/)
{
  archive &kind;
}
 
 
 
inline ArrayView<const unsigned int>
ReferenceCell::faces_for_given_vertex(const unsigned int vertex) const
{
  if (*this == ReferenceCells::Line)
    {
      AssertIndexRange(vertex, GeometryInfo<1>::vertices_per_cell);
      return {&GeometryInfo<2>::vertex_to_face[vertex][0], 1};
    }
  else if (*this == ReferenceCells::Quadrilateral)
    {
      AssertIndexRange(vertex, GeometryInfo<2>::vertices_per_cell);
      return {&GeometryInfo<2>::vertex_to_face[vertex][0], 2};
    }
  else if (*this == ReferenceCells::Hexahedron)
    {
      AssertIndexRange(vertex, GeometryInfo<3>::vertices_per_cell);
      return {&GeometryInfo<3>::vertex_to_face[vertex][0], 3};
    }
  else if (*this == ReferenceCells::Triangle)
    {
      AssertIndexRange(vertex, 3);
      static const ndarray<unsigned int, 3, 2> table = {
        {{{0, 2}}, {{0, 1}}, {{1, 2}}}};
 
      return table[vertex];
    }
  else if (*this == ReferenceCells::Tetrahedron)
    {
      AssertIndexRange(vertex, 4);
      static const ndarray<unsigned int, 4, 3> table = {
        {{{0, 1, 2}}, {{0, 1, 3}}, {{0, 2, 3}}, {{1, 2, 3}}}};
 
      return table[vertex];
    }
  else if (*this == ReferenceCells::Wedge)
    {
      AssertIndexRange(vertex, 6);
      static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 4}},
                                                         {{0, 2, 3}},
                                                         {{0, 3, 4}},
                                                         {{1, 2, 4}},
                                                         {{1, 2, 3}},
                                                         {{1, 3, 4}}}};
 
      return table[vertex];
    }
  else if (*this == ReferenceCells::Pyramid)
    {
      AssertIndexRange(vertex, 5);
      static const unsigned int X = numbers::invalid_unsigned_int;
      static const ndarray<unsigned int, 5, 4> table = {{{{0, 1, 3, X}},
                                                         {{0, 2, 3, X}},
                                                         {{0, 1, 4, X}},
                                                         {{0, 2, 4, X}},
                                                         {{1, 2, 3, 4}}}};
 
      return {&table[vertex][0], vertex == 4 ? 4u : 3u};
    }
 
  Assert(false, ExcNotImplemented());
 
  return {};
}
 
 
 
inline bool
ReferenceCell::is_hyper_cube() const
{
  return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
          *this == ReferenceCells::Quadrilateral ||
          *this == ReferenceCells::Hexahedron);
}
 
 
 
inline bool
ReferenceCell::is_simplex() const
{
  return (*this == ReferenceCells::Vertex || *this == ReferenceCells::Line ||
          *this == ReferenceCells::Triangle ||
          *this == ReferenceCells::Tetrahedron);
}
 
 
 
inline unsigned int
ReferenceCell::get_dimension() const
{
  if (*this == ReferenceCells::Vertex)
    return 0;
  else if (*this == ReferenceCells::Line)
    return 1;
  else if ((*this == ReferenceCells::Triangle) ||
           (*this == ReferenceCells::Quadrilateral))
    return 2;
  else if ((*this == ReferenceCells::Tetrahedron) ||
           (*this == ReferenceCells::Pyramid) ||
           (*this == ReferenceCells::Wedge) ||
           (*this == ReferenceCells::Hexahedron))
    return 3;
 
  Assert(false, ExcNotImplemented());
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int dim>
Quadrature<dim>
ReferenceCell::get_midpoint_quadrature() const
{
  return Quadrature<dim>(std::vector<Point<dim>>({barycenter<dim>()}),
                         std::vector<double>({volume()}));
}
 
 
 
inline unsigned int
ReferenceCell::n_vertices() const
{
  if (*this == ReferenceCells::Vertex)
    return 1;
  else if (*this == ReferenceCells::Line)
    return 2;
  else if (*this == ReferenceCells::Triangle)
    return 3;
  else if (*this == ReferenceCells::Quadrilateral)
    return 4;
  else if (*this == ReferenceCells::Tetrahedron)
    return 4;
  else if (*this == ReferenceCells::Pyramid)
    return 5;
  else if (*this == ReferenceCells::Wedge)
    return 6;
  else if (*this == ReferenceCells::Hexahedron)
    return 8;
 
  Assert(false, ExcNotImplemented());
  return numbers::invalid_unsigned_int;
}
 
 
 
inline unsigned int
ReferenceCell::n_lines() const
{
  if (*this == ReferenceCells::Vertex)
    return 0;
  else if (*this == ReferenceCells::Line)
    return 1;
  else if (*this == ReferenceCells::Triangle)
    return 3;
  else if (*this == ReferenceCells::Quadrilateral)
    return 4;
  else if (*this == ReferenceCells::Tetrahedron)
    return 6;
  else if (*this == ReferenceCells::Pyramid)
    return 7;
  else if (*this == ReferenceCells::Wedge)
    return 9;
  else if (*this == ReferenceCells::Hexahedron)
    return 12;
 
  Assert(false, ExcNotImplemented());
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int dim>
Point<dim>
ReferenceCell::vertex(const unsigned int v) const
{
  AssertDimension(dim, get_dimension());
  AssertIndexRange(v, n_vertices());
 
  if ((dim == 0) && (*this == ReferenceCells::Vertex))
    {
      return Point<dim>(0);
    }
  else if ((dim == 1) && (*this == ReferenceCells::Line))
    {
      static const Point<dim> vertices[2] = {
        Point<dim>(),              // the origin
        Point<dim>::unit_vector(0) // unit point along x-axis
      };
      return vertices[v];
    }
  else if ((dim == 2) && (*this == ReferenceCells::Quadrilateral))
    {
      static const Point<dim> vertices[4] = {
        // First the two points on the x-axis
        Point<dim>(),
        Point<dim>::unit_vector(0),
        // Then these two points shifted in the y-direction
        Point<dim>() + Point<dim>::unit_vector(1),
        Point<dim>::unit_vector(0) + Point<dim>::unit_vector(1)};
      return vertices[v];
    }
  else if ((dim == 3) && (*this == ReferenceCells::Hexahedron))
    {
      static const Point<dim> vertices[8] = {
        // First the two points on the x-axis
        Point<dim>(),
        Point<dim>::unit_vector(0),
        // Then these two points shifted in the y-direction
        Point<dim>() + Point<dim>::unit_vector(1),
        Point<dim>::unit_vector(0) + Point<dim>::unit_vector(1),
        // And now all four points shifted in the z-direction
        Point<dim>() + Point<dim>::unit_vector(2),
        Point<dim>::unit_vector(0) + Point<dim>::unit_vector(2),
        Point<dim>() + Point<dim>::unit_vector(1) + Point<dim>::unit_vector(2),
        Point<dim>::unit_vector(0) + Point<dim>::unit_vector(1) +
          Point<dim>::unit_vector(2)};
      return vertices[v];
    }
  else if ((dim == 2) && (*this == ReferenceCells::Triangle))
    {
      static const Point<dim> vertices[3] = {
        Point<dim>(),               // the origin
        Point<dim>::unit_vector(0), // unit point along x-axis
        Point<dim>::unit_vector(1)  // unit point along y-axis
      };
      return vertices[v];
    }
  else if ((dim == 3) && (*this == ReferenceCells::Tetrahedron))
    {
      static const Point<dim> vertices[4] = {
        Point<dim>(),               // the origin
        Point<dim>::unit_vector(0), // unit point along x-axis
        Point<dim>::unit_vector(1), // unit point along y-axis
        Point<dim>::unit_vector(2)  // unit point along z-axis
      };
      return vertices[v];
    }
  else if ((dim == 3) && (*this == ReferenceCells::Pyramid))
    {
      static const Point<dim> vertices[5] = {Point<dim>{-1.0, -1.0, 0.0},
                                             Point<dim>{+1.0, -1.0, 0.0},
                                             Point<dim>{-1.0, +1.0, 0.0},
                                             Point<dim>{+1.0, +1.0, 0.0},
                                             Point<dim>{+0.0, +0.0, 1.0}};
      return vertices[v];
    }
  else if ((dim == 3) && (*this == ReferenceCells::Wedge))
    {
      static const Point<dim> vertices[6] = {
        // First the three points on the triangular base of the wedge:
        Point<dim>(),
        Point<dim>::unit_vector(0),
        Point<dim>::unit_vector(1),
        // And now everything shifted in the z-direction again
        Point<dim>() + Point<dim>::unit_vector(2),
        Point<dim>::unit_vector(0) + Point<dim>::unit_vector(2),
        Point<dim>::unit_vector(1) + Point<dim>::unit_vector(2)};
      return vertices[v];
    }
  else
    {
      Assert(false, ExcNotImplemented());
      return Point<dim>();
    }
}
 
 
inline unsigned int
ReferenceCell::n_faces() const
{
  if (*this == ReferenceCells::Vertex)
    return 0;
  else if (*this == ReferenceCells::Line)
    return 2;
  else if (*this == ReferenceCells::Triangle)
    return 3;
  else if (*this == ReferenceCells::Quadrilateral)
    return 4;
  else if (*this == ReferenceCells::Tetrahedron)
    return 4;
  else if (*this == ReferenceCells::Pyramid)
    return 5;
  else if (*this == ReferenceCells::Wedge)
    return 5;
  else if (*this == ReferenceCells::Hexahedron)
    return 6;
 
  Assert(false, ExcNotImplemented());
  return numbers::invalid_unsigned_int;
}
 
 
 
inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
ReferenceCell::face_indices() const
{
  return {0U, n_faces()};
}
 
 
 
inline unsigned int
ReferenceCell::n_isotropic_children() const
{
  if (*this == ReferenceCells::Vertex)
    return 0;
  else if (*this == ReferenceCells::Line)
    return 2;
  else if (*this == ReferenceCells::Triangle)
    return 4;
  else if (*this == ReferenceCells::Quadrilateral)
    return 4;
  else if (*this == ReferenceCells::Tetrahedron)
    return 8;
  else if (*this == ReferenceCells::Pyramid)
    {
      // We haven't yet decided how to refine pyramids. Update
      // this when we have
      Assert(false, ExcNotImplemented());
      return numbers::invalid_unsigned_int;
    }
  else if (*this == ReferenceCells::Wedge)
    return 8;
  else if (*this == ReferenceCells::Hexahedron)
    return 8;
 
  Assert(false, ExcNotImplemented());
  return numbers::invalid_unsigned_int;
}
 
 
 
inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
ReferenceCell::isotropic_child_indices() const
{
  return {0U, n_isotropic_children()};
}
 
 
 
inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
ReferenceCell::vertex_indices() const
{
  return {0U, n_vertices()};
}
 
 
 
inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
ReferenceCell::line_indices() const
{
  return {0U, n_lines()};
}
 
 
 
inline ReferenceCell
ReferenceCell::face_reference_cell(const unsigned int face_no) const
{
  AssertIndexRange(face_no, n_faces());
 
  if (*this == ReferenceCells::Vertex)
    return ReferenceCells::Invalid;
  else if (*this == ReferenceCells::Line)
    return ReferenceCells::Vertex;
  else if (*this == ReferenceCells::Triangle)
    return ReferenceCells::Line;
  else if (*this == ReferenceCells::Quadrilateral)
    return ReferenceCells::Line;
  else if (*this == ReferenceCells::Tetrahedron)
    return ReferenceCells::Triangle;
  else if (*this == ReferenceCells::Pyramid)
    {
      if (face_no == 0)
        return ReferenceCells::Quadrilateral;
      else
        return ReferenceCells::Triangle;
    }
  else if (*this == ReferenceCells::Wedge)
    {
      if (face_no > 1)
        return ReferenceCells::Quadrilateral;
      else
        return ReferenceCells::Triangle;
    }
  else if (*this == ReferenceCells::Hexahedron)
    return ReferenceCells::Quadrilateral;
 
  Assert(false, ExcNotImplemented());
  return ReferenceCells::Invalid;
}
 
 
 
inline unsigned int
ReferenceCell::child_cell_on_face(
  const unsigned int  face,
  const unsigned int  subface,
  const unsigned char face_orientation_raw) const
{
  AssertIndexRange(face, n_faces());
  AssertIndexRange(subface, face_reference_cell(face).n_isotropic_children());
 
  if (*this == ReferenceCells::Vertex)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Line)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Triangle)
    {
      static const ndarray<unsigned int, 3, 2> subcells = {
        {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
 
      return subcells[face][subface];
    }
  else if (*this == ReferenceCells::Quadrilateral)
    {
      const bool face_orientation = Utilities::get_bit(face_orientation_raw, 0);
      const bool face_flip        = Utilities::get_bit(face_orientation_raw, 2);
      const bool face_rotation    = Utilities::get_bit(face_orientation_raw, 1);
 
      return GeometryInfo<2>::child_cell_on_face(
        RefinementCase<2>(RefinementPossibilities<2>::isotropic_refinement),
        face,
        subface,
        face_orientation,
        face_flip,
        face_rotation);
    }
  else if (*this == ReferenceCells::Tetrahedron)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Pyramid)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Wedge)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Hexahedron)
    {
      const bool face_orientation = Utilities::get_bit(face_orientation_raw, 0);
      const bool face_flip        = Utilities::get_bit(face_orientation_raw, 2);
      const bool face_rotation    = Utilities::get_bit(face_orientation_raw, 1);
 
      return GeometryInfo<3>::child_cell_on_face(
        RefinementCase<3>(RefinementPossibilities<3>::isotropic_refinement),
        face,
        subface,
        face_orientation,
        face_flip,
        face_rotation);
    }
 
  Assert(false, ExcNotImplemented());
  return {};
}
 
 
 
inline std::array<unsigned int, 2>
ReferenceCell::standard_vertex_to_face_and_vertex_index(
  const unsigned int vertex) const
{
  AssertIndexRange(vertex, n_vertices());
  // Work around a GCC warning at higher optimization levels by making all of
  // these tables the same size
  constexpr unsigned int X = numbers::invalid_unsigned_int;
 
  if (*this == ReferenceCells::Vertex)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Line)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Triangle)
    {
      static const ndarray<unsigned int, 6, 2> table = {
        {{{0, 0}}, {{0, 1}}, {{1, 1}}, {{X, X}}, {{X, X}}, {{X, X}}}};
 
      return table[vertex];
    }
  else if (*this == ReferenceCells::Quadrilateral)
    {
      return GeometryInfo<2>::standard_quad_vertex_to_line_vertex_index(vertex);
    }
  else if (*this == ReferenceCells::Tetrahedron)
    {
      static const ndarray<unsigned int, 6, 2> table = {
        {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 2}}, {{X, X}}, {{X, X}}}};
 
      return table[vertex];
    }
  else if (*this == ReferenceCells::Pyramid)
    {
      static const ndarray<unsigned int, 6, 2> table = {
        {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{0, 3}}, {{1, 2}}, {{X, X}}}};
 
      return table[vertex];
    }
  else if (*this == ReferenceCells::Wedge)
    {
      static const ndarray<unsigned int, 6, 2> table = {
        {{{0, 1}}, {{0, 0}}, {{0, 2}}, {{1, 0}}, {{1, 1}}, {{1, 2}}}};
 
      return table[vertex];
    }
  else if (*this == ReferenceCells::Hexahedron)
    {
      return GeometryInfo<3>::standard_hex_vertex_to_quad_vertex_index(vertex);
    }
 
  Assert(false, ExcNotImplemented());
  return {};
}
 
 
 
inline std::array<unsigned int, 2>
ReferenceCell::standard_line_to_face_and_line_index(
  const unsigned int line) const
{
  AssertIndexRange(line, n_lines());
 
  if (*this == ReferenceCells::Vertex)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Line)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Triangle)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Quadrilateral)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Tetrahedron)
    {
      static const std::array<unsigned int, 2> table[6] = {
        {{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}};
 
      return table[line];
    }
  else if (*this == ReferenceCells::Pyramid)
    {
      static const std::array<unsigned int, 2> table[8] = {{{0, 0}},
                                                           {{0, 1}},
                                                           {{0, 2}},
                                                           {{0, 3}},
                                                           {{1, 2}},
                                                           {{2, 1}},
                                                           {{1, 1}},
                                                           {{2, 2}}};
 
      return table[line];
    }
  else if (*this == ReferenceCells::Wedge)
    {
      static const std::array<unsigned int, 2> table[9] = {{{0, 0}},
                                                           {{0, 2}},
                                                           {{0, 1}},
                                                           {{1, 0}},
                                                           {{1, 1}},
                                                           {{1, 2}},
                                                           {{2, 0}},
                                                           {{2, 1}},
                                                           {{3, 1}}};
 
      return table[line];
    }
  else if (*this == ReferenceCells::Hexahedron)
    {
      return GeometryInfo<3>::standard_hex_line_to_quad_line_index(line);
    }
 
  Assert(false, ExcNotImplemented());
  return {};
}
 
 
 
inline unsigned int
ReferenceCell::face_to_cell_lines(const unsigned int  face,
                                  const unsigned int  line,
                                  const unsigned char face_orientation) const
{
  AssertIndexRange(face, n_faces());
  AssertIndexRange(line, face_reference_cell(face).n_lines());
 
  static constexpr unsigned int X = numbers::invalid_unsigned_int;
 
  if (*this == ReferenceCells::Vertex)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Line)
    {
      return GeometryInfo<1>::face_to_cell_lines(
        face,
        line,
        Utilities::get_bit(face_orientation, 0),
        Utilities::get_bit(face_orientation, 2),
        Utilities::get_bit(face_orientation, 1));
    }
  else if (*this == ReferenceCells::Triangle)
    {
      return face;
    }
  else if (*this == ReferenceCells::Quadrilateral)
    {
      return GeometryInfo<2>::face_to_cell_lines(
        face,
        line,
        Utilities::get_bit(face_orientation, 0),
        Utilities::get_bit(face_orientation, 2),
        Utilities::get_bit(face_orientation, 1));
    }
  else if (*this == ReferenceCells::Tetrahedron)
    {
      const static ndarray<unsigned int, 4, 3> table = {
        {{{0, 1, 2}}, {{0, 3, 4}}, {{2, 5, 3}}, {{1, 4, 5}}}};
 
      return table[face]
                  [standard_to_real_face_line(line, face, face_orientation)];
    }
  else if (*this == ReferenceCells::Pyramid)
    {
      static const ndarray<unsigned int, 5, 4> table = {{{{0, 1, 2, 3}},
                                                         {{0, 6, 4, X}},
                                                         {{1, 5, 7, X}},
                                                         {{2, 4, 5, X}},
                                                         {{3, 7, 6, 2}}}};
 
      return table[face]
                  [standard_to_real_face_line(line, face, face_orientation)];
    }
  else if (*this == ReferenceCells::Wedge)
    {
      static const ndarray<unsigned int, 5, 4> table = {{{{0, 2, 1, X}},
                                                         {{3, 4, 5, X}},
                                                         {{6, 7, 0, 3}},
                                                         {{7, 8, 1, 4}},
                                                         {{8, 6, 5, 2}}}};
 
      return table[face]
                  [standard_to_real_face_line(line, face, face_orientation)];
    }
  else if (*this == ReferenceCells::Hexahedron)
    {
      return GeometryInfo<3>::face_to_cell_lines(
        face,
        line,
        Utilities::get_bit(face_orientation, 0),
        Utilities::get_bit(face_orientation, 2),
        Utilities::get_bit(face_orientation, 1));
    }
 
  Assert(false, ExcNotImplemented());
  return numbers::invalid_unsigned_int;
}
 
 
 
inline unsigned int
ReferenceCell::face_to_cell_vertices(const unsigned int  face,
                                     const unsigned int  vertex,
                                     const unsigned char face_orientation) const
{
  AssertIndexRange(face, n_faces());
  AssertIndexRange(vertex, face_reference_cell(face).n_vertices());
 
  if (*this == ReferenceCells::Vertex)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Line)
    {
      return GeometryInfo<1>::face_to_cell_vertices(
        face,
        vertex,
        Utilities::get_bit(face_orientation, 0),
        Utilities::get_bit(face_orientation, 2),
        Utilities::get_bit(face_orientation, 1));
    }
  else if (*this == ReferenceCells::Triangle)
    {
      static const ndarray<unsigned int, 3, 2> table = {
        {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
 
      return table[face][face_orientation != 0u ? vertex : (1 - vertex)];
    }
  else if (*this == ReferenceCells::Quadrilateral)
    {
      return GeometryInfo<2>::face_to_cell_vertices(
        face,
        vertex,
        Utilities::get_bit(face_orientation, 0),
        Utilities::get_bit(face_orientation, 2),
        Utilities::get_bit(face_orientation, 1));
    }
  else if (*this == ReferenceCells::Tetrahedron)
    {
      static const ndarray<unsigned int, 4, 3> table = {
        {{{0, 1, 2}}, {{1, 0, 3}}, {{0, 2, 3}}, {{2, 1, 3}}}};
 
      return table[face][standard_to_real_face_vertex(
        vertex, face, face_orientation)];
    }
  else if (*this == ReferenceCells::Pyramid)
    {
      constexpr auto X = numbers::invalid_unsigned_int;
      static const ndarray<unsigned int, 5, 4> table = {{{{0, 1, 2, 3}},
                                                         {{0, 2, 4, X}},
                                                         {{3, 1, 4, X}},
                                                         {{1, 0, 4, X}},
                                                         {{2, 3, 4, X}}}};
 
      return table[face][standard_to_real_face_vertex(
        vertex, face, face_orientation)];
    }
  else if (*this == ReferenceCells::Wedge)
    {
      constexpr auto X = numbers::invalid_unsigned_int;
      static const ndarray<unsigned int, 6, 4> table = {{{{1, 0, 2, X}},
                                                         {{3, 4, 5, X}},
                                                         {{0, 1, 3, 4}},
                                                         {{1, 2, 4, 5}},
                                                         {{2, 0, 5, 3}}}};
 
      return table[face][standard_to_real_face_vertex(
        vertex, face, face_orientation)];
    }
  else if (*this == ReferenceCells::Hexahedron)
    {
      return GeometryInfo<3>::face_to_cell_vertices(
        face,
        vertex,
        Utilities::get_bit(face_orientation, 0),
        Utilities::get_bit(face_orientation, 2),
        Utilities::get_bit(face_orientation, 1));
    }
 
  Assert(false, ExcNotImplemented());
  return numbers::invalid_unsigned_int;
}
 
 
 
inline unsigned int
ReferenceCell::standard_to_real_face_vertex(
  const unsigned int  vertex,
  const unsigned int  face,
  const unsigned char face_orientation) const
{
  AssertIndexRange(face, n_faces());
  AssertIndexRange(vertex, face_reference_cell(face).n_vertices());
 
  if (*this == ReferenceCells::Vertex)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Line)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Triangle)
    {
      static const ndarray<unsigned int, 2, 2> table = {{{{1, 0}}, {{0, 1}}}};
 
      return table[face_orientation][vertex];
    }
  else if (*this == ReferenceCells::Quadrilateral)
    {
      return GeometryInfo<2>::standard_to_real_line_vertex(vertex,
                                                           face_orientation !=
                                                             0u);
    }
  else if (*this == ReferenceCells::Tetrahedron)
    {
      static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 1}},
                                                         {{0, 1, 2}},
                                                         {{2, 1, 0}},
                                                         {{1, 2, 0}},
                                                         {{1, 0, 2}},
                                                         {{2, 0, 1}}}};
 
      return table[face_orientation][vertex];
    }
  else if (*this == ReferenceCells::Pyramid)
    {
      if (face == 0) // The quadrilateral face
        {
          return GeometryInfo<3>::standard_to_real_face_vertex(
            vertex,
            Utilities::get_bit(face_orientation, 0),
            Utilities::get_bit(face_orientation, 2),
            Utilities::get_bit(face_orientation, 1));
        }
      else // One of the triangular faces
        {
          static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 1}},
                                                             {{0, 1, 2}},
                                                             {{2, 1, 0}},
                                                             {{1, 2, 0}},
                                                             {{1, 0, 2}},
                                                             {{2, 0, 1}}}};
 
          return table[face_orientation][vertex];
        }
    }
  else if (*this == ReferenceCells::Wedge)
    {
      if (face > 1) // One of the quadrilateral faces
        {
          return GeometryInfo<3>::standard_to_real_face_vertex(
            vertex,
            Utilities::get_bit(face_orientation, 0),
            Utilities::get_bit(face_orientation, 2),
            Utilities::get_bit(face_orientation, 1));
        }
      else // One of the triangular faces
        {
          static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 1}},
                                                             {{0, 1, 2}},
                                                             {{2, 1, 0}},
                                                             {{1, 2, 0}},
                                                             {{1, 0, 2}},
                                                             {{2, 0, 1}}}};
 
          return table[face_orientation][vertex];
        }
    }
  else if (*this == ReferenceCells::Hexahedron)
    {
      return GeometryInfo<3>::standard_to_real_face_vertex(
        vertex,
        Utilities::get_bit(face_orientation, 0),
        Utilities::get_bit(face_orientation, 2),
        Utilities::get_bit(face_orientation, 1));
    }
 
  Assert(false, ExcNotImplemented());
  return numbers::invalid_unsigned_int;
}
 
 
 
inline unsigned int
ReferenceCell::standard_to_real_face_line(
  const unsigned int  line,
  const unsigned int  face,
  const unsigned char face_orientation) const
{
  AssertIndexRange(face, n_faces());
  AssertIndexRange(line, face_reference_cell(face).n_lines());
 
  if (*this == ReferenceCells::Vertex)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Line)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Triangle)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Quadrilateral)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Tetrahedron)
    {
      static const ndarray<unsigned int, 6, 3> table = {{{{2, 1, 0}},
                                                         {{0, 1, 2}},
                                                         {{1, 0, 2}},
                                                         {{1, 2, 0}},
                                                         {{0, 2, 1}},
                                                         {{2, 0, 1}}}};
 
      return table[face_orientation][line];
    }
  else if (*this == ReferenceCells::Pyramid)
    {
      if (face == 0) // The quadrilateral face
        {
          return GeometryInfo<3>::standard_to_real_face_line(
            line,
            Utilities::get_bit(face_orientation, 0),
            Utilities::get_bit(face_orientation, 2),
            Utilities::get_bit(face_orientation, 1));
        }
      else // One of the triangular faces
        {
          static const ndarray<unsigned int, 6, 3> table = {{{{2, 1, 0}},
                                                             {{0, 1, 2}},
                                                             {{1, 0, 2}},
                                                             {{1, 2, 0}},
                                                             {{0, 2, 1}},
                                                             {{2, 0, 1}}}};
 
          return table[face_orientation][line];
        }
    }
  else if (*this == ReferenceCells::Wedge)
    {
      if (face > 1) // One of the quadrilateral faces
        {
          return GeometryInfo<3>::standard_to_real_face_line(
            line,
            Utilities::get_bit(face_orientation, 0),
            Utilities::get_bit(face_orientation, 2),
            Utilities::get_bit(face_orientation, 1));
        }
      else // One of the triangular faces
        {
          static const ndarray<unsigned int, 6, 3> table = {{{{2, 1, 0}},
                                                             {{0, 1, 2}},
                                                             {{1, 0, 2}},
                                                             {{1, 2, 0}},
                                                             {{0, 2, 1}},
                                                             {{2, 0, 1}}}};
 
          return table[face_orientation][line];
        }
    }
  else if (*this == ReferenceCells::Hexahedron)
    {
      return GeometryInfo<3>::standard_to_real_face_line(
        line,
        Utilities::get_bit(face_orientation, 0),
        Utilities::get_bit(face_orientation, 2),
        Utilities::get_bit(face_orientation, 1));
    }
 
  Assert(false, ExcNotImplemented());
  return numbers::invalid_unsigned_int;
}
 
 
 
namespace ReferenceCells
{
  template <int dim>
  inline constexpr const ReferenceCell &
  get_simplex()
  {
    switch (dim)
      {
        case 0:
          return ReferenceCells::Vertex;
        case 1:
          return ReferenceCells::Line;
        case 2:
          return ReferenceCells::Triangle;
        case 3:
          return ReferenceCells::Tetrahedron;
        default:
          Assert(false, ExcNotImplemented());
          return ReferenceCells::Invalid;
      }
  }
 
 
 
  template <int dim>
  inline constexpr const ReferenceCell &
  get_hypercube()
  {
    switch (dim)
      {
        case 0:
          return ReferenceCells::Vertex;
        case 1:
          return ReferenceCells::Line;
        case 2:
          return ReferenceCells::Quadrilateral;
        case 3:
          return ReferenceCells::Hexahedron;
        default:
          Assert(false, ExcNotImplemented());
          return ReferenceCells::Invalid;
      }
  }
} // namespace ReferenceCells
 
 
inline ReferenceCell
ReferenceCell::n_vertices_to_type(const int dim, const unsigned int n_vertices)
{
  AssertIndexRange(dim, 4);
  AssertIndexRange(n_vertices, 9);
 
  const auto                                X     = ReferenceCells::Invalid;
  static const ndarray<ReferenceCell, 4, 9> table = {
    {// dim 0
     {{X, ReferenceCells::Vertex, X, X, X, X, X, X, X}},
     // dim 1
     {{X, X, ReferenceCells::Line, X, X, X, X, X, X}},
     // dim 2
     {{X,
       X,
       X,
       ReferenceCells::Triangle,
       ReferenceCells::Quadrilateral,
       X,
       X,
       X,
       X}},
     // dim 3
     {{X,
       X,
       X,
       X,
       ReferenceCells::Tetrahedron,
       ReferenceCells::Pyramid,
       ReferenceCells::Wedge,
       X,
       ReferenceCells::Hexahedron}}}};
  Assert(table[dim][n_vertices] != ReferenceCells::Invalid,
         ExcMessage("The combination of dim = " + std::to_string(dim) +
                    " and n_vertices = " + std::to_string(n_vertices) +
                    " does not correspond to a known reference cell type."));
  return table[dim][n_vertices];
}
 
 
 
template <int dim>
inline double
ReferenceCell::d_linear_shape_function(const Point<dim> & xi,
                                       const unsigned int i) const
{
  AssertDimension(dim, get_dimension());
  if (*this == ReferenceCells::get_hypercube<dim>())
    return GeometryInfo<dim>::d_linear_shape_function(xi, i);
 
  if (*this ==
      ReferenceCells::Triangle) // see also
                                // BarycentricPolynomials<2>::compute_value
    {
      switch (i)
        {
          case 0:
            return 1.0 - xi[std::min(0, dim - 1)] - xi[std::min(1, dim - 1)];
          case 1:
            return xi[std::min(0, dim - 1)];
          case 2:
            return xi[std::min(1, dim - 1)];
        }
    }
 
  if (*this ==
      ReferenceCells::Tetrahedron) // see also
                                   // BarycentricPolynomials<3>::compute_value
    {
      switch (i)
        {
          case 0:
            return 1.0 - xi[std::min(0, dim - 1)] - xi[std::min(1, dim - 1)] -
                   xi[std::min(2, dim - 1)];
          case 1:
            return xi[std::min(0, dim - 1)];
          case 2:
            return xi[std::min(1, dim - 1)];
          case 3:
            return xi[std::min(2, dim - 1)];
        }
    }
 
  if (*this ==
      ReferenceCells::Wedge) // see also
                             // ScalarLagrangePolynomialWedge::compute_value
    {
      return ReferenceCell(ReferenceCells::Triangle)
               .d_linear_shape_function<2>(Point<2>(xi[std::min(0, dim - 1)],
                                                    xi[std::min(1, dim - 1)]),
                                           i % 3) *
             ReferenceCell(ReferenceCells::Line)
               .d_linear_shape_function<1>(Point<1>(xi[std::min(2, dim - 1)]),
                                           i / 3);
    }
 
  if (*this ==
      ReferenceCells::Pyramid) // see also
                               // ScalarLagrangePolynomialPyramid::compute_value
    {
      const double Q14 = 0.25;
      double       ration;
 
      const double r = xi[std::min(0, dim - 1)];
      const double s = xi[std::min(1, dim - 1)];
      const double t = xi[std::min(2, dim - 1)];
 
      if (fabs(t - 1.0) > 1.0e-14)
        {
          ration = (r * s * t) / (1.0 - t);
        }
      else
        {
          ration = 0.0;
        }
 
      if (i == 0)
        return Q14 * ((1.0 - r) * (1.0 - s) - t + ration);
      if (i == 1)
        return Q14 * ((1.0 + r) * (1.0 - s) - t - ration);
      if (i == 2)
        return Q14 * ((1.0 - r) * (1.0 + s) - t - ration);
      if (i == 3)
        return Q14 * ((1.0 + r) * (1.0 + s) - t + ration);
      else
        return t;
    }
 
  Assert(false, ExcNotImplemented());
 
  return 0.0;
}
 
 
 
template <int dim>
inline Tensor<1, dim>
ReferenceCell::d_linear_shape_function_gradient(const Point<dim> & xi,
                                                const unsigned int i) const
{
  AssertDimension(dim, get_dimension());
  if (*this == ReferenceCells::get_hypercube<dim>())
    return GeometryInfo<dim>::d_linear_shape_function_gradient(xi, i);
 
  if (*this ==
      ReferenceCells::Triangle) // see also
                                // BarycentricPolynomials<2>::compute_grad
    {
      switch (i)
        {
          case 0:
            return Point<dim>(-1.0, -1.0);
          case 1:
            return Point<dim>(+1.0, +0.0);
          case 2:
            return Point<dim>(+0.0, +1.0);
        }
    }
 
  Assert(false, ExcNotImplemented());
 
  return Point<dim>(+0.0, +0.0, +0.0);
}
 
 
 
inline double
ReferenceCell::volume() const
{
  if (*this == ReferenceCells::Vertex)
    return 0;
  else if (*this == ReferenceCells::Line)
    return 1;
  else if (*this == ReferenceCells::Triangle)
    return 1. / 2.;
  else if (*this == ReferenceCells::Quadrilateral)
    return 1;
  else if (*this == ReferenceCells::Tetrahedron)
    return 1. / 6.;
  else if (*this == ReferenceCells::Wedge)
    return 1. / 2.;
  else if (*this == ReferenceCells::Pyramid)
    return 4. / 3.;
  else if (*this == ReferenceCells::Hexahedron)
    return 1;
 
  Assert(false, ExcNotImplemented());
  return 0.0;
}
 
 
 
template <int dim>
inline Point<dim>
ReferenceCell::barycenter() const
{
  AssertDimension(dim, get_dimension());
 
  if (*this == ReferenceCells::Vertex)
    return Point<dim>();
  else if (*this == ReferenceCells::Line)
    return Point<dim>(1. / 2.);
  else if (*this == ReferenceCells::Triangle)
    return Point<dim>(1. / 3., 1. / 3.);
  else if (*this == ReferenceCells::Quadrilateral)
    return Point<dim>(1. / 2., 1. / 2.);
  else if (*this == ReferenceCells::Tetrahedron)
    return Point<dim>(1. / 4., 1. / 4., 1. / 4.);
  else if (*this == ReferenceCells::Wedge)
    return Point<dim>(1. / 3, 1. / 3, 1. / 2.);
  else if (*this == ReferenceCells::Pyramid)
    return Point<dim>(0, 0, 1. / 4.);
  else if (*this == ReferenceCells::Hexahedron)
    return Point<dim>(1. / 2., 1. / 2., 1. / 2.);
 
  Assert(false, ExcNotImplemented());
  return Point<dim>();
}
 
 
 
template <int dim>
inline bool
ReferenceCell::contains_point(const Point<dim> &p, const double tolerance) const
{
  AssertDimension(dim, get_dimension());
 
  if (*this == ReferenceCells::Vertex)
    {
      // Vertices are special cases in that they do not actually
      // have coordinates. Error out if this function is called
      // with a vertex:
      Assert(false,
             ExcMessage("Vertices are zero-dimensional objects and "
                        "as a consequence have no coordinates. You "
                        "cannot meaningfully ask whether a point is "
                        "inside a vertex (within a certain tolerance) "
                        "without coordinate values."));
      return false;
    }
  else if (*this == ReferenceCells::get_hypercube<dim>())
    {
      for (unsigned int d = 0; d < dim; ++d)
        if ((p[d] < -tolerance) || (p[d] > 1 + tolerance))
          return false;
      return true;
    }
  else if (*this == ReferenceCells::get_simplex<dim>())
    {
      // First make sure that we are in the first quadrant or octant
      for (unsigned int d = 0; d < dim; ++d)
        if (p[d] < -tolerance)
          return false;
 
      // Now we also need to make sure that we are below the diagonal line
      // or plane that delineates the simplex. This diagonal is given by
      // sum(p[d])<=1, and a diagonal a distance eps away is given by
      // sum(p[d])<=1+eps*sqrt(d). (For example, the point at (1,1) is a
      // distance of 1/sqrt(2) away from the diagonal. That is, its
      // sum satisfies
      //   sum(p[d]) = 2 <= 1 + (1/sqrt(2)) * sqrt(2)
      // in other words, it satisfies the predicate with eps=1/sqrt(2).)
      double sum = 0;
      for (unsigned int d = 0; d < dim; ++d)
        sum += p[d];
      return (sum <= 1 + tolerance * std::sqrt(1. * dim));
    }
  else if (*this == ReferenceCells::Wedge)
    {
      Assert(false, ExcNotImplemented());
    }
  else if (*this == ReferenceCells::Pyramid)
    {
      Assert(false, ExcNotImplemented());
    }
 
  Assert(false, ExcNotImplemented());
 
  return false;
}
 
 
 
template <int dim>
inline Tensor<1, dim>
ReferenceCell::unit_tangential_vectors(const unsigned int face_no,
                                       const unsigned int i) const
{
  AssertDimension(dim, get_dimension());
  AssertIndexRange(i, dim - 1);
 
  if (*this == ReferenceCells::get_hypercube<dim>())
    {
      AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
      return GeometryInfo<dim>::unit_tangential_vectors[face_no][i];
    }
  else if (*this == ReferenceCells::Triangle)
    {
      AssertIndexRange(face_no, 3);
      static const std::array<Tensor<1, dim>, 3> table = {
        {Point<dim>(1, 0),
         Point<dim>(-std::sqrt(0.5), +std::sqrt(0.5)),
         Point<dim>(0, -1)}};
 
      return table[face_no];
    }
  else if (*this == ReferenceCells::Tetrahedron)
    {
      AssertIndexRange(face_no, 4);
      static const ndarray<Tensor<1, dim>, 4, 2> table = {
        {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
         {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
         {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
         {{Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
                      +std::pow(1.0 / 3.0, 1.0 / 4.0),
                      0),
           Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
                      0,
                      +std::pow(1.0 / 3.0, 1.0 / 4.0))}}}};
 
      return table[face_no][i];
    }
  else if (*this == ReferenceCells::Wedge)
    {
      AssertIndexRange(face_no, 5);
      static const ndarray<Tensor<1, dim>, 5, 2> table = {
        {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
         {{Point<dim>(1, 0, 0), Point<dim>(0, 1, 0)}},
         {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
         {{Point<dim>(-1 / std::sqrt(2.0), +1 / std::sqrt(2.0), 0),
           Point<dim>(0, 0, 1)}},
         {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}}}};
 
      return table[face_no][i];
    }
  else if (*this == ReferenceCells::Pyramid)
    {
      AssertIndexRange(face_no, 5);
      static const ndarray<Tensor<1, dim>, 5, 2> table = {
        {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
         {{Point<dim>(+1.0 / sqrt(2.0), 0, +1.0 / sqrt(2.0)),
           Point<dim>(0, 1, 0)}},
         {{Point<dim>(+1.0 / sqrt(2.0), 0, -1.0 / sqrt(2.0)),
           Point<dim>(0, 1, 0)}},
         {{Point<dim>(1, 0, 0),
           Point<dim>(0, +1.0 / sqrt(2.0), +1.0 / sqrt(2.0))}},
         {{Point<dim>(1, 0, 0),
           Point<dim>(0, +1.0 / sqrt(2.0), -1.0 / sqrt(2.0))}}}};
 
      return table[face_no][i];
    }
 
  Assert(false, ExcNotImplemented());
 
  return {};
}
 
 
 
template <int dim>
inline Tensor<1, dim>
ReferenceCell::unit_normal_vectors(const unsigned int face_no) const
{
  AssertDimension(dim, this->get_dimension());
 
  if (is_hyper_cube())
    {
      AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
      return GeometryInfo<dim>::unit_normal_vector[face_no];
    }
  else if (dim == 2)
    {
      Assert(*this == ReferenceCells::Triangle, ExcInternalError());
 
      // Return the rotated vector
      return cross_product_2d(unit_tangential_vectors<dim>(face_no, 0));
    }
  else if (dim == 3)
    {
      return cross_product_3d(unit_tangential_vectors<dim>(face_no, 0),
                              unit_tangential_vectors<dim>(face_no, 1));
    }
 
  Assert(false, ExcNotImplemented());
 
  return {};
}
 
 
 
inline bool
ReferenceCell::standard_vs_true_line_orientation(
  const unsigned int  line,
  const unsigned char face_orientation_raw,
  const bool          line_orientation) const
{
  if (*this == ReferenceCells::Hexahedron)
    {
      static const bool bool_table[2][2][2][2] = {
        {{{true, false},    // lines 0/1, face_orientation=false,
                            // face_flip=false, face_rotation=false and true
          {false, true}},   // lines 0/1, face_orientation=false,
                            // face_flip=true, face_rotation=false and true
         {{true, true},     // lines 0/1, face_orientation=true,
                            // face_flip=false, face_rotation=false and true
          {false, false}}}, // lines 0/1, face_orientation=true,
                            // face_flip=true, face_rotation=false and true
 
        {{{true, true}, // lines 2/3 ...
          {false, false}},
         {{true, false}, {false, true}}}};
 
      const bool face_orientation = Utilities::get_bit(face_orientation_raw, 0);
      const bool face_flip        = Utilities::get_bit(face_orientation_raw, 2);
      const bool face_rotation    = Utilities::get_bit(face_orientation_raw, 1);
 
      return (line_orientation ==
              bool_table[line / 2][face_orientation][face_flip][face_rotation]);
    }
  else
    // TODO: This might actually be wrong for some of the other
    // kinds of objects. We should check this
    return true;
}
 
 
 
namespace internal
{
  template <typename T, std::size_t N>
  class NoPermutation : public ::ExceptionBase
  {
  public:
    NoPermutation(const ::ReferenceCell &entity_type,
                  const std::array<T, N> &     vertices_0,
                  const std::array<T, N> &     vertices_1)
      : entity_type(entity_type)
      , vertices_0(vertices_0)
      , vertices_1(vertices_1)
    {}
 
    virtual ~NoPermutation() noexcept override = default;
 
    virtual void
    print_info(std::ostream &out) const override
    {
      out << '[';
 
      const unsigned int n_vertices = entity_type.n_vertices();
 
      for (unsigned int i = 0; i < n_vertices; ++i)
        {
          out << vertices_0[i];
          if (i + 1 != n_vertices)
            out << ',';
        }
 
      out << "] is not a permutation of [";
 
      for (unsigned int i = 0; i < n_vertices; ++i)
        {
          out << vertices_1[i];
          if (i + 1 != n_vertices)
            out << ',';
        }
 
      out << "]." << std::endl;
    }
 
    const ::ReferenceCell entity_type;
 
    const std::array<T, N> vertices_0;
 
    const std::array<T, N> vertices_1;
  };
} // namespace internal
 
 
 
template <typename T, std::size_t N>
inline unsigned char
ReferenceCell::compute_orientation(const std::array<T, N> &vertices_0,
                                   const std::array<T, N> &vertices_1) const
{
  AssertIndexRange(n_vertices(), N + 1);
  if (*this == ReferenceCells::Line)
    {
      const std::array<T, 2> i{{vertices_0[0], vertices_0[1]}};
      const std::array<T, 2> j{{vertices_1[0], vertices_1[1]}};
 
      // line_orientation=true
      if (i == std::array<T, 2>{{j[0], j[1]}})
        return 1;
 
      // line_orientation=false
      if (i == std::array<T, 2>{{j[1], j[0]}})
        return 0;
    }
  else if (*this == ReferenceCells::Triangle)
    {
      const std::array<T, 3> i{{vertices_0[0], vertices_0[1], vertices_0[2]}};
      const std::array<T, 3> j{{vertices_1[0], vertices_1[1], vertices_1[2]}};
 
      // face_orientation=true, face_rotation=false, face_flip=false
      if (i == std::array<T, 3>{{j[0], j[1], j[2]}})
        return 1;
 
      // face_orientation=true, face_rotation=true, face_flip=false
      if (i == std::array<T, 3>{{j[1], j[2], j[0]}})
        return 3;
 
      // face_orientation=true, face_rotation=false, face_flip=true
      if (i == std::array<T, 3>{{j[2], j[0], j[1]}})
        return 5;
 
      // face_orientation=false, face_rotation=false, face_flip=false
      if (i == std::array<T, 3>{{j[0], j[2], j[1]}})
        return 0;
 
      // face_orientation=false, face_rotation=true, face_flip=false
      if (i == std::array<T, 3>{{j[2], j[1], j[0]}})
        return 2;
 
      // face_orientation=false, face_rotation=false, face_flip=true
      if (i == std::array<T, 3>{{j[1], j[0], j[2]}})
        return 4;
    }
  else if (*this == ReferenceCells::Quadrilateral)
    {
      const std::array<T, 4> i{
        {vertices_0[0], vertices_0[1], vertices_0[2], vertices_0[3]}};
      const std::array<T, 4> j{
        {vertices_1[0], vertices_1[1], vertices_1[2], vertices_1[3]}};
 
      // face_orientation=true, face_rotation=false, face_flip=false
      if (i == std::array<T, 4>{{j[0], j[1], j[2], j[3]}})
        return 1;
 
      // face_orientation=true, face_rotation=true, face_flip=false
      if (i == std::array<T, 4>{{j[2], j[0], j[3], j[1]}})
        return 3;
 
      // face_orientation=true, face_rotation=false, face_flip=true
      if (i == std::array<T, 4>{{j[3], j[2], j[1], j[0]}})
        return 5;
 
      // face_orientation=true, face_rotation=true, face_flip=true
      if (i == std::array<T, 4>{{j[1], j[3], j[0], j[2]}})
        return 7;
 
      // face_orientation=false, face_rotation=false, face_flip=false
      if (i == std::array<T, 4>{{j[0], j[2], j[1], j[3]}})
        return 0;
 
      // face_orientation=false, face_rotation=true, face_flip=false
      if (i == std::array<T, 4>{{j[2], j[3], j[0], j[1]}})
        return 2;
 
      // face_orientation=false, face_rotation=false, face_flip=true
      if (i == std::array<T, 4>{{j[3], j[1], j[2], j[0]}})
        return 4;
 
      // face_orientation=false, face_rotation=true, face_flip=true
      if (i == std::array<T, 4>{{j[1], j[0], j[3], j[2]}})
        return 6;
    }
 
  Assert(false, (internal::NoPermutation<T, N>(*this, vertices_0, vertices_1)));
 
  return -1;
}
 
 
 
template <typename T, std::size_t N>
inline std::array<T, N>
ReferenceCell::permute_according_orientation(
  const std::array<T, N> &vertices,
  const unsigned int      orientation) const
{
  std::array<T, 4> temp;
 
  if (*this == ReferenceCells::Line)
    {
      switch (orientation)
        {
          case 1:
            temp = {{vertices[0], vertices[1]}};
            break;
          case 0:
            temp = {{vertices[1], vertices[0]}};
            break;
          default:
            Assert(false, ExcNotImplemented());
        }
    }
  else if (*this == ReferenceCells::Triangle)
    {
      switch (orientation)
        {
          case 1:
            temp = {{vertices[0], vertices[1], vertices[2]}};
            break;
          case 3:
            temp = {{vertices[1], vertices[2], vertices[0]}};
            break;
          case 5:
            temp = {{vertices[2], vertices[0], vertices[1]}};
            break;
          case 0:
            temp = {{vertices[0], vertices[2], vertices[1]}};
            break;
          case 2:
            temp = {{vertices[2], vertices[1], vertices[0]}};
            break;
          case 4:
            temp = {{vertices[1], vertices[0], vertices[2]}};
            break;
          default:
            Assert(false, ExcNotImplemented());
        }
    }
  else if (*this == ReferenceCells::Quadrilateral)
    {
      switch (orientation)
        {
          case 1:
            temp = {{vertices[0], vertices[1], vertices[2], vertices[3]}};
            break;
          case 3:
            temp = {{vertices[2], vertices[0], vertices[3], vertices[1]}};
            break;
          case 5:
            temp = {{vertices[3], vertices[2], vertices[1], vertices[0]}};
            break;
          case 7:
            temp = {{vertices[1], vertices[3], vertices[0], vertices[2]}};
            break;
          case 0:
            temp = {{vertices[0], vertices[2], vertices[1], vertices[3]}};
            break;
          case 2:
            temp = {{vertices[2], vertices[3], vertices[0], vertices[1]}};
            break;
          case 4:
            temp = {{vertices[3], vertices[1], vertices[2], vertices[0]}};
            break;
          case 6:
            temp = {{vertices[1], vertices[0], vertices[3], vertices[2]}};
            break;
          default:
            Assert(false, ExcNotImplemented());
        }
    }
  else
    {
      AssertThrow(false, ExcNotImplemented());
    }
 
  std::array<T, N> temp_;
  std::copy_n(temp.begin(), N, temp_.begin());
 
  return temp_;
}
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif