reference_cell.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2020 - 2023 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 <boost/container/small_vector.hpp>
 
#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
{
  constexpr ReferenceCell
  make_reference_cell_from_int(const std::uint8_t kind);
} // namespace internal
 
 
 
class ReferenceCell
{
public:
  static ReferenceCell
  n_vertices_to_type(const int dim, const unsigned int n_vertices);
 
  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;
 
  static constexpr unsigned char
  default_combined_face_orientation();
 
  static constexpr unsigned char
  reversed_combined_line_orientation();
 
  unsigned int
  child_cell_on_face(const unsigned int  face,
                     const unsigned int  subface,
                     const unsigned char face_orientation =
                       default_combined_face_orientation()) 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
  line_to_cell_vertices(const unsigned int line,
                        const unsigned int vertex) 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;
 
  template <int dim>
  Point<dim>
  face_vertex_location(const unsigned int face,
                       const unsigned int vertex) 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;
 
  template <int dim>
  Point<dim>
  closest_point(const Point<dim> &p) const;
 
  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;
 
  unsigned int
  n_face_orientations(const unsigned int face_no) const;
 
  template <typename T, std::size_t N>
  DEAL_II_DEPRECATED_EARLY unsigned char
  compute_orientation(const std::array<T, N> &vertices_0,
                      const std::array<T, N> &vertices_1) const;
 
  template <typename T>
  unsigned char
  get_combined_orientation(const ArrayView<const T> &vertices_0,
                           const ArrayView<const T> &vertices_1) const;
 
 
  template <typename T, std::size_t N>
  DEAL_II_DEPRECATED_EARLY std::array<T, N>
  permute_according_orientation(const std::array<T, N> &vertices,
                                const unsigned int      orientation) const;
 
  template <typename T>
  boost::container::small_vector<T, 8>
  permute_by_combined_orientation(const ArrayView<const T> &vertices,
                                  const unsigned char       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;
 
  template <int dim>
  unsigned int
  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;
 
  unsigned int
  vtk_vertex_to_deal_vertex(const unsigned int vertex_index) 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*/);
 
  static constexpr std::size_t
  memory_consumption();
 
private:
  std::uint8_t kind;
 
  constexpr ReferenceCell(const std::uint8_t kind);
 
  static constexpr ndarray<unsigned int, 2, 2> line_vertex_permutations = {
    {{{1, 0}}, {{0, 1}}}};
 
 
  static constexpr ndarray<unsigned int, 6, 3> triangle_vertex_permutations = {
    {{{0, 2, 1}},
     {{0, 1, 2}},
     {{2, 1, 0}},
     {{1, 2, 0}},
     {{1, 0, 2}},
     {{2, 0, 1}}}};
 
  static constexpr ndarray<unsigned int, 8, 4>
    quadrilateral_vertex_permutations = {{{{0, 2, 1, 3}},
                                          {{0, 1, 2, 3}},
                                          {{2, 3, 0, 1}},
                                          {{2, 0, 3, 1}},
                                          {{3, 1, 2, 0}},
                                          {{3, 2, 1, 0}},
                                          {{1, 0, 3, 2}},
                                          {{1, 3, 0, 2}}}};
 
  friend constexpr ReferenceCell
  internal::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
{
  inline constexpr ReferenceCell
  make_reference_cell_from_int(const std::uint8_t kind)
  {
#ifndef DEAL_II_CXX14_CONSTEXPR_BUG
    // Make sure these are the only indices from which objects can be
    // created.
    Assert((kind == std::numeric_limits<std::uint8_t>::max()) || (kind < 8),
           ExcInternalError());
#endif
 
    // Call the private constructor, which we can from here because this
    // function is a 'friend'.
    return {kind};
  }
} // namespace internal
 
 
 
namespace ReferenceCells
{
  constexpr const ReferenceCell Vertex =
    internal::make_reference_cell_from_int(0);
  constexpr const ReferenceCell Line =
    internal::make_reference_cell_from_int(1);
  constexpr const ReferenceCell Triangle =
    internal::make_reference_cell_from_int(2);
  constexpr const ReferenceCell Quadrilateral =
    internal::make_reference_cell_from_int(3);
  constexpr const ReferenceCell Tetrahedron =
    internal::make_reference_cell_from_int(4);
  constexpr const ReferenceCell Pyramid =
    internal::make_reference_cell_from_int(5);
  constexpr const ReferenceCell Wedge =
    internal::make_reference_cell_from_int(6);
  constexpr const ReferenceCell Hexahedron =
    internal::make_reference_cell_from_int(7);
  constexpr const ReferenceCell Invalid =
    internal::make_reference_cell_from_int(
      std::numeric_limits<std::uint8_t>::max());
 
  template <int dim>
  constexpr const ReferenceCell &
  get_simplex();
 
  template <int dim>
  constexpr const ReferenceCell &
  get_hypercube();
} // namespace ReferenceCells
 
 
 
inline constexpr ReferenceCell::ReferenceCell()
  : ReferenceCell(ReferenceCells::Invalid)
{}
 
 
 
template <class Archive>
inline void
ReferenceCell::serialize(Archive &archive, const unsigned int /*version*/)
{
  archive &kind;
}
 
 
 
inline constexpr std::size_t
ReferenceCell::memory_consumption()
{
  return sizeof(ReferenceCells::Invalid);
}
 
 
 
inline ArrayView<const unsigned int>
ReferenceCell::faces_for_given_vertex(const unsigned int vertex) const
{
  AssertIndexRange(vertex, n_vertices());
  switch (this->kind)
    {
      case ReferenceCells::Line:
        return {&GeometryInfo<2>::vertex_to_face[vertex][0], 1};
      case ReferenceCells::Quadrilateral:
        return {&GeometryInfo<2>::vertex_to_face[vertex][0], 2};
      case ReferenceCells::Triangle:
        {
          static constexpr ndarray<unsigned int, 3, 2> table = {
            {{{0, 2}}, {{0, 1}}, {{1, 2}}}};
          return table[vertex];
        }
      case ReferenceCells::Tetrahedron:
        {
          static constexpr ndarray<unsigned int, 4, 3> table = {
            {{{0, 1, 2}}, {{0, 1, 3}}, {{0, 2, 3}}, {{1, 2, 3}}}};
 
          return table[vertex];
        }
      case ReferenceCells::Pyramid:
        {
          static constexpr unsigned int X = numbers::invalid_unsigned_int;
          static constexpr 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};
        }
      case ReferenceCells::Wedge:
        {
          AssertIndexRange(vertex, 6);
          static constexpr 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];
        }
      case ReferenceCells::Hexahedron:
        return {&GeometryInfo<3>::vertex_to_face[vertex][0], 3};
      default:
        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
{
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        return 0;
      case ReferenceCells::Line:
        return 1;
      case ReferenceCells::Triangle:
      case ReferenceCells::Quadrilateral:
        return 2;
      case ReferenceCells::Tetrahedron:
      case ReferenceCells::Pyramid:
      case ReferenceCells::Wedge:
      case ReferenceCells::Hexahedron:
        return 3;
      default:
        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
{
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        return 1;
      case ReferenceCells::Line:
        return 2;
      case ReferenceCells::Triangle:
        return 3;
      case ReferenceCells::Quadrilateral:
        return 4;
      case ReferenceCells::Tetrahedron:
        return 4;
      case ReferenceCells::Pyramid:
        return 5;
      case ReferenceCells::Wedge:
        return 6;
      case ReferenceCells::Hexahedron:
        return 8;
      default:
        Assert(false, ExcNotImplemented());
    }
 
  return numbers::invalid_unsigned_int;
}
 
 
 
inline unsigned int
ReferenceCell::n_lines() const
{
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        return 0;
      case ReferenceCells::Line:
        return 1;
      case ReferenceCells::Triangle:
        return 3;
      case ReferenceCells::Quadrilateral:
        return 4;
      case ReferenceCells::Tetrahedron:
        return 6;
      case ReferenceCells::Pyramid:
        return 8;
      case ReferenceCells::Wedge:
        return 9;
      case ReferenceCells::Hexahedron:
        return 12;
      default:
        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());
 
  switch (dim)
    {
      case 0:
        {
          if (*this == ReferenceCells::Vertex)
            return Point<dim>();
          break;
        }
      case 1:
        {
          static const Point<dim> vertices[2] = {
            Point<dim>(),              // the origin
            Point<dim>::unit_vector(0) // unit point along x-axis
          };
          if (*this == ReferenceCells::Line)
            return vertices[v];
          break;
        }
      case 2:
        {
          switch (this->kind)
            {
              case 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];
                }
              case 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];
                }
            }
          break;
        }
      case 3:
        {
          switch (this->kind)
            {
              case 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];
                }
              case 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];
                }
              case 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];
                }
              case 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];
                }
            }
          break;
        }
      default:
        Assert(false, ExcNotImplemented());
    }
 
  Assert(false, ExcNotImplemented());
  return Point<dim>();
}
 
 
inline unsigned int
ReferenceCell::n_faces() const
{
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        return 0;
      case ReferenceCells::Line:
        return 2;
      case ReferenceCells::Triangle:
        return 3;
      case ReferenceCells::Quadrilateral:
        return 4;
      case ReferenceCells::Tetrahedron:
        return 4;
      case ReferenceCells::Pyramid:
        return 5;
      case ReferenceCells::Wedge:
        return 5;
      case ReferenceCells::Hexahedron:
        return 6;
      default:
        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
{
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        return 0;
      case ReferenceCells::Line:
        return 2;
      case ReferenceCells::Triangle:
        return 4;
      case ReferenceCells::Quadrilateral:
        return 4;
      case ReferenceCells::Tetrahedron:
        return 8;
      case ReferenceCells::Pyramid:
        // We haven't yet decided how to refine pyramids. Update this when we
        // have
        Assert(false, ExcNotImplemented());
        return numbers::invalid_unsigned_int;
      case ReferenceCells::Wedge:
        return 8;
      case ReferenceCells::Hexahedron:
        return 8;
      default:
        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());
 
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        return ReferenceCells::Invalid;
      case ReferenceCells::Line:
        return ReferenceCells::Vertex;
      case ReferenceCells::Triangle:
      case ReferenceCells::Quadrilateral:
        return ReferenceCells::Line;
      case ReferenceCells::Tetrahedron:
        return ReferenceCells::Triangle;
      case ReferenceCells::Pyramid:
        if (face_no == 0)
          return ReferenceCells::Quadrilateral;
        else
          return ReferenceCells::Triangle;
      case ReferenceCells::Wedge:
        if (face_no > 1)
          return ReferenceCells::Quadrilateral;
        else
          return ReferenceCells::Triangle;
      case ReferenceCells::Hexahedron:
        return ReferenceCells::Quadrilateral;
      default:
        Assert(false, ExcNotImplemented());
    }
 
  return ReferenceCells::Invalid;
}
 
 
 
inline constexpr unsigned char
ReferenceCell::default_combined_face_orientation()
{
  // Our convention is that 'orientation' has a default value of true and
  // occupies the least-significant bit while rotate and flip have default
  // values of 'false' and occupy the second and third bits.
  return 0b001;
}
 
 
 
inline constexpr unsigned char
ReferenceCell::reversed_combined_line_orientation()
{
  // For a reversed line 'orientation' is false and neither flip nor rotate are
  // defined.
  return 0b000;
}
 
 
 
inline unsigned int
ReferenceCell::child_cell_on_face(
  const unsigned int  face,
  const unsigned int  subface,
  const unsigned char combined_face_orientation) const
{
  AssertIndexRange(face, n_faces());
  AssertIndexRange(subface, face_reference_cell(face).n_isotropic_children());
 
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
      case ReferenceCells::Line:
        {
          Assert(false, ExcNotImplemented());
          break;
        }
      case ReferenceCells::Triangle:
        {
          static constexpr ndarray<unsigned int, 3, 2> subcells = {
            {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
 
          return subcells[face][subface];
        }
      case ReferenceCells::Quadrilateral:
        {
          const bool face_orientation =
            Utilities::get_bit(combined_face_orientation, 0);
          const bool face_flip =
            Utilities::get_bit(combined_face_orientation, 2);
          const bool face_rotation =
            Utilities::get_bit(combined_face_orientation, 1);
 
          return GeometryInfo<2>::child_cell_on_face(
            RefinementCase<2>(RefinementPossibilities<2>::isotropic_refinement),
            face,
            subface,
            face_orientation,
            face_flip,
            face_rotation);
        }
      case ReferenceCells::Tetrahedron:
      case ReferenceCells::Pyramid:
      case ReferenceCells::Wedge:
        {
          Assert(false, ExcNotImplemented());
          break;
        }
      case ReferenceCells::Hexahedron:
        {
          const bool face_orientation =
            Utilities::get_bit(combined_face_orientation, 0);
          const bool face_flip =
            Utilities::get_bit(combined_face_orientation, 2);
          const bool face_rotation =
            Utilities::get_bit(combined_face_orientation, 1);
 
          return GeometryInfo<3>::child_cell_on_face(
            RefinementCase<3>(RefinementPossibilities<3>::isotropic_refinement),
            face,
            subface,
            face_orientation,
            face_flip,
            face_rotation);
        }
      default:
        Assert(false, ExcNotImplemented());
    }
 
  return numbers::invalid_unsigned_int;
}
 
 
 
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;
 
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
      case ReferenceCells::Line:
        Assert(false, ExcNotImplemented());
        break;
      case ReferenceCells::Triangle:
        {
          static constexpr ndarray<unsigned int, 6, 2> table = {
            {{{0, 0}}, {{0, 1}}, {{1, 1}}, {{X, X}}, {{X, X}}, {{X, X}}}};
 
          return table[vertex];
        }
      case ReferenceCells::Quadrilateral:
        {
          return GeometryInfo<2>::standard_quad_vertex_to_line_vertex_index(
            vertex);
        }
      case ReferenceCells::Tetrahedron:
        {
          static constexpr ndarray<unsigned int, 6, 2> table = {
            {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 2}}, {{X, X}}, {{X, X}}}};
 
          return table[vertex];
        }
      case ReferenceCells::Pyramid:
        {
          static constexpr ndarray<unsigned int, 6, 2> table = {
            {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{0, 3}}, {{1, 2}}, {{X, X}}}};
 
          return table[vertex];
        }
      case ReferenceCells::Wedge:
        {
          static constexpr ndarray<unsigned int, 6, 2> table = {
            {{{0, 1}}, {{0, 0}}, {{0, 2}}, {{1, 0}}, {{1, 1}}, {{1, 2}}}};
 
          return table[vertex];
        }
      case ReferenceCells::Hexahedron:
        {
          return GeometryInfo<3>::standard_hex_vertex_to_quad_vertex_index(
            vertex);
        }
      default:
        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());
 
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
      case ReferenceCells::Line:
      case ReferenceCells::Triangle:
      case ReferenceCells::Quadrilateral:
        {
          Assert(false, ExcNotImplemented());
          break;
        }
      case 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];
        }
      case 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];
        }
      case 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];
        }
      case ReferenceCells::Hexahedron:
        {
          return GeometryInfo<3>::standard_hex_line_to_quad_line_index(line);
        }
      default:
        Assert(false, ExcNotImplemented());
    }
 
  return {};
}
 
 
inline unsigned int
ReferenceCell::line_to_cell_vertices(const unsigned int line,
                                     const unsigned int vertex) const
{
  AssertIndexRange(vertex, 2);
  AssertIndexRange(line, n_lines());
 
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
      case ReferenceCells::Line:
        return vertex;
      case ReferenceCells::Triangle:
        {
          static constexpr ndarray<unsigned int, 3, 2> table = {
            {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
          return table[line][vertex];
        }
      case ReferenceCells::Quadrilateral:
        {
          static constexpr ndarray<unsigned int, 4, 2> table = {
            {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
          return table[line][vertex];
        }
      case ReferenceCells::Tetrahedron:
        {
          static constexpr ndarray<unsigned int, 6, 2> table = {
            {{{0, 1}}, {{1, 2}}, {{2, 0}}, {{0, 3}}, {{1, 3}}, {{2, 3}}}};
          return table[line][vertex];
        }
      case ReferenceCells::Pyramid:
        {
          static constexpr ndarray<unsigned int, 8, 2> table = {{{{0, 2}},
                                                                 {{1, 3}},
                                                                 {{0, 1}},
                                                                 {{2, 3}},
                                                                 {{4, 0}},
                                                                 {{1, 4}},
                                                                 {{2, 4}},
                                                                 {{4, 3}}}};
          return table[line][vertex];
        }
      case ReferenceCells::Wedge:
        {
          static constexpr ndarray<unsigned int, 9, 2> table = {{{{1, 0}},
                                                                 {{2, 1}},
                                                                 {{0, 2}},
                                                                 {{3, 4}},
                                                                 {{4, 5}},
                                                                 {{5, 3}},
                                                                 {{0, 3}},
                                                                 {{1, 4}},
                                                                 {{2, 5}}}};
          return table[line][vertex];
        }
      case ReferenceCells::Hexahedron:
        {
          // first four lines comprise the bottom face, next four are the top,
          // and the last four are 'bottom to top'
          static constexpr ndarray<unsigned int, 12, 2> table = {{{{0, 2}},
                                                                  {{1, 3}},
                                                                  {{0, 1}},
                                                                  {{2, 3}},
                                                                  {{4, 6}},
                                                                  {{5, 7}},
                                                                  {{4, 5}},
                                                                  {{6, 7}},
                                                                  {{0, 4}},
                                                                  {{1, 5}},
                                                                  {{2, 6}},
                                                                  {{3, 7}}}};
          return table[line][vertex];
        }
 
      default:
        Assert(false, ExcNotImplemented());
    }
 
  return numbers::invalid_unsigned_int;
}
 
 
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;
 
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        {
          Assert(false, ExcNotImplemented());
          break;
        }
      case 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));
        }
      case ReferenceCells::Triangle:
        {
          return face;
        }
      case 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));
        }
      case ReferenceCells::Tetrahedron:
        {
          static constexpr 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)];
        }
      case ReferenceCells::Pyramid:
        {
          static constexpr 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, X}}}};
 
          return table[face][standard_to_real_face_line(
            line, face, face_orientation)];
        }
      case ReferenceCells::Wedge:
        {
          static constexpr 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)];
        }
      case 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));
        }
      default:
        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());
 
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        {
          Assert(false, ExcNotImplemented());
          break;
        }
      case 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));
        }
      case ReferenceCells::Triangle:
        {
          static constexpr ndarray<unsigned int, 3, 2> table = {
            {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
 
          return table[face][face_orientation != 0u ? vertex : (1 - vertex)];
        }
      case 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));
        }
      case ReferenceCells::Tetrahedron:
        {
          static constexpr 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)];
        }
      case ReferenceCells::Pyramid:
        {
          constexpr auto X = numbers::invalid_unsigned_int;
          static constexpr 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)];
        }
      case ReferenceCells::Wedge:
        {
          constexpr auto X = numbers::invalid_unsigned_int;
          static constexpr 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)];
        }
      case 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));
        }
      default:
        Assert(false, ExcNotImplemented());
    }
 
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int dim>
Point<dim>
ReferenceCell::face_vertex_location(const unsigned int face,
                                    const unsigned int vertex) const
{
  AssertIndexRange(face, n_faces());
  AssertIndexRange(vertex, face_reference_cell(face).n_vertices());
  AssertDimension(dim, this->get_dimension());
 
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        {
          Assert(false, ExcNotImplemented());
          break;
        }
      case ReferenceCells::Line:
        {
          if (face == 0)
            return Point<dim>(0);
          else
            return Point<dim>(1);
        }
      case ReferenceCells::Triangle:
        {
          static const ndarray<Point<dim>, 3, 2> table = {
            {{{Point<dim>(0, 0), Point<dim>(1, 0)}},
             {{Point<dim>(1, 0), Point<dim>(0, 1)}},
             {{Point<dim>(0, 1), Point<dim>(0, 0)}}}};
 
          return table[face][vertex];
        }
      case ReferenceCells::Quadrilateral:
        {
          static const ndarray<Point<dim>, 4, 2> table = {
            {{{Point<dim>(0, 0), Point<dim>(0, 1)}},
             {{Point<dim>(1, 0), Point<dim>(1, 1)}},
             {{Point<dim>(0, 0), Point<dim>(1, 0)}},
             {{Point<dim>(0, 1), Point<dim>(1, 1)}}}};
 
          return table[face][vertex];
        }
      case ReferenceCells::Tetrahedron:
        {
          static const ndarray<Point<dim>, 4, 3> table = {
            {{{Point<dim>(0.0, 0.0, 0.0),
               Point<dim>(1.0, 0.0, 0.0),
               Point<dim>(0.0, 1.0, 0.0)}},
             {{Point<dim>(1.0, 0.0, 0.0),
               Point<dim>(0.0, 0.0, 0.0),
               Point<dim>(0.0, 0.0, 1.0)}},
             {{Point<dim>(0.0, 0.0, 0.0),
               Point<dim>(0.0, 1.0, 0.0),
               Point<dim>(0.0, 0.0, 1.0)}},
             {{Point<dim>(0.0, 1.0, 0.0),
               Point<dim>(1.0, 0.0, 0.0),
               Point<dim>(0.0, 0.0, 1.0)}}}};
 
          return table[face][vertex];
        }
      case ReferenceCells::Pyramid:
        {
          // We would like to use ndarray here as we do above, but
          // not every face has the same number of vertices and so
          // the inner array needs to be a dynamic vector:
          static const std::array<std::vector<Point<dim>>, 5> table = {
            {{{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>(-1.0, -1.0, 0.0),
               Point<dim>(-1.0, +1.0, 0.0),
               Point<dim>(+0.0, +0.0, 1.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)}},
             {{Point<dim>(+1.0, -1.0, 0.0),
               Point<dim>(-1.0, -1.0, 0.0),
               Point<dim>(+0.0, +0.0, 1.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 table[face][vertex];
        }
      case ReferenceCells::Wedge:
        {
          // We would like to use ndarray here as we do above, but
          // not every face has the same number of vertices and so
          // the inner array needs to be a dynamic vector:
          static const std::array<std::vector<Point<dim>>, 5> table = {
            {{{Point<dim>(1.0, 0.0, 0.0),
               Point<dim>(0.0, 0.0, 0.0),
               Point<dim>(0.0, 1.0, 0.0)}},
             {{Point<dim>(0.0, 0.0, 1.0),
               Point<dim>(1.0, 0.0, 1.0),
               Point<dim>(0.0, 1.0, 1.0)}},
             {{Point<dim>(0.0, 0.0, 0.0),
               Point<dim>(1.0, 0.0, 0.0),
               Point<dim>(0.0, 0.0, 1.0),
               Point<dim>(1.0, 0.0, 1.0)}},
             {{Point<dim>(1.0, 0.0, 0.0),
               Point<dim>(0.0, 1.0, 0.0),
               Point<dim>(1.0, 0.0, 1.0),
               Point<dim>(0.0, 1.0, 1.0)}},
             {{Point<dim>(0.0, 1.0, 0.0),
               Point<dim>(0.0, 0.0, 0.0),
               Point<dim>(0.0, 1.0, 1.0),
               Point<dim>(0.0, 0.0, 1.0)}}}};
 
          return table[face][vertex];
        }
      case ReferenceCells::Hexahedron:
        {
          /*
           * Follow these figures:
           *
           *  * <li> Faces 0 and 1:
           *  @verbatim
           *          Face 0           Face 1
           *        *-------*        *-------*
           *       /|       |       /       /|
           *      3 1       |      /       3 1
           *    y/  |       |     /      y/  |
           *    *   |x      |    *-------*   |x
           *    |   *-------*    |       |   *
           *    0  /       /     |       0  /
           *    | 2       /      |       | 2
           *    |/       /       |       |/
           *    *-------*        *-------*
           *  @endverbatim
           *
           * <li> Faces 2 and 3:
           *  @verbatim
           *        x Face 3           Face 2
           *        *---1---*        *-------*
           *       /|       |       /       /|
           *      / |       3      /       / |
           *     /  2       |    x/       /  |
           *    *   |       |    *---1---*   |
           *    |   *---0---*y   |       |   *
           *    |  /       /     |       3  /
           *    | /       /      2       | /
           *    |/       /       |       |/
           *    *-------*        *---0---*y
           *  @endverbatim
           *
           * <li> Faces 4 and 5:
           *  @verbatim
           *          Face 4         y Face 5
           *        *-------*        *---3---*
           *       /|       |       /       /|
           *      / |       |      0       1 |
           *     /  |       |     /       /  |
           *    *   |y      |    *---2---* x |
           *    |   *---3---*    |       |   *
           *    |  /       /     |       |  /
           *    | 0       1      |       | /
           *    |/       /       |       |/
           *    *---2---* x      *-------*
           *  @endverbatim
           * </ul>
           */
          static const ndarray<Point<dim>, 6, 4> table = {
            {// Face 0:
             {{Point<dim>(0, 0, 0),
               Point<dim>(0, 1, 0),
               Point<dim>(0, 0, 1),
               Point<dim>(0, 1, 1)}},
             // Face 1:
             {{Point<dim>(1, 0, 0),
               Point<dim>(1, 1, 0),
               Point<dim>(1, 0, 1),
               Point<dim>(1, 1, 1)}},
             // Face 2:
             {{Point<dim>(0, 0, 0),
               Point<dim>(0, 0, 1),
               Point<dim>(1, 0, 0),
               Point<dim>(1, 0, 1)}},
             // Face 3:
             {{Point<dim>(0, 1, 0),
               Point<dim>(0, 1, 1),
               Point<dim>(1, 1, 0),
               Point<dim>(1, 1, 1)}},
             // Face 4:
             {{Point<dim>(0, 0, 0),
               Point<dim>(1, 0, 0),
               Point<dim>(0, 1, 0),
               Point<dim>(1, 1, 0)}},
             // Face 5:
             {{Point<dim>(0, 0, 1),
               Point<dim>(1, 0, 1),
               Point<dim>(0, 1, 1),
               Point<dim>(1, 1, 1)}}}};
 
          return table[face][vertex];
        }
      default:
        Assert(false, ExcNotImplemented());
    }
 
  return {};
}
 
 
 
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());
 
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
      case ReferenceCells::Line:
        Assert(false, ExcNotImplemented());
        break;
      case ReferenceCells::Triangle:
      case ReferenceCells::Quadrilateral:
        return line_vertex_permutations[face_orientation][vertex];
      case ReferenceCells::Tetrahedron:
        return triangle_vertex_permutations[face_orientation][vertex];
      case ReferenceCells::Pyramid:
        // face 0 is a quadrilateral
        if (face == 0)
          return quadrilateral_vertex_permutations[face_orientation][vertex];
        else
          return triangle_vertex_permutations[face_orientation][vertex];
      case ReferenceCells::Wedge:
        // faces 0 and 1 are triangles
        if (face > 1)
          return quadrilateral_vertex_permutations[face_orientation][vertex];
        else
          return triangle_vertex_permutations[face_orientation][vertex];
      case ReferenceCells::Hexahedron:
        return quadrilateral_vertex_permutations[face_orientation][vertex];
      default:
        Assert(false, ExcNotImplemented());
    }
 
  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());
 
  static constexpr ndarray<unsigned int, 6, 3> triangle_table = {{{{2, 1, 0}},
                                                                  {{0, 1, 2}},
                                                                  {{1, 0, 2}},
                                                                  {{1, 2, 0}},
                                                                  {{0, 2, 1}},
                                                                  {{2, 0, 1}}}};
 
 
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
      case ReferenceCells::Line:
      case ReferenceCells::Triangle:
      case ReferenceCells::Quadrilateral:
        Assert(false, ExcNotImplemented());
        break;
      case ReferenceCells::Tetrahedron:
        return triangle_table[face_orientation][line];
      case 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
          {
            return triangle_table[face_orientation][line];
          }
      case 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
          return triangle_table[face_orientation][line];
      case ReferenceCells::Hexahedron:
        {
          static constexpr ndarray<unsigned int, 8, 4> table = {
            {{{2, 3, 0, 1}},
             {{0, 1, 2, 3}},
             {{0, 1, 3, 2}},
             {{3, 2, 0, 1}},
             {{3, 2, 1, 0}},
             {{1, 0, 3, 2}},
             {{1, 0, 2, 3}},
             {{2, 3, 1, 0}}}};
          return table[face_orientation][line];
        }
      default:
        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());
  AssertIndexRange(i, n_vertices());
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
      case ReferenceCells::Line:
      case ReferenceCells::Quadrilateral:
      case ReferenceCells::Hexahedron:
        return GeometryInfo<dim>::d_linear_shape_function(xi, i);
      // see also BarycentricPolynomials<2>::compute_value
      case ReferenceCells::Triangle:
        {
          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)];
              default:
                Assert(false, ExcInternalError());
            }
        }
      // see also BarycentricPolynomials<3>::compute_value
      case ReferenceCells::Tetrahedron:
        {
          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)];
              default:
                Assert(false, ExcInternalError());
            }
        }
      // see also ScalarLagrangePolynomialPyramid::compute_value()
      case ReferenceCells::Pyramid:
        {
          const double Q14 = 0.25;
 
          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)];
 
          const double ratio =
            (std::fabs(t - 1.0) > 1.0e-14 ? (r * s * t) / (1.0 - t) : 0.0);
 
          if (i == 0)
            return Q14 * ((1.0 - r) * (1.0 - s) - t + ratio);
          if (i == 1)
            return Q14 * ((1.0 + r) * (1.0 - s) - t - ratio);
          if (i == 2)
            return Q14 * ((1.0 - r) * (1.0 + s) - t - ratio);
          if (i == 3)
            return Q14 * ((1.0 + r) * (1.0 + s) - t + ratio);
          else
            return t;
        }
      // see also ScalarLagrangePolynomialWedge::compute_value()
      case ReferenceCells::Wedge:
        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);
      default:
        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());
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
      case ReferenceCells::Line:
      case ReferenceCells::Quadrilateral:
      case ReferenceCells::Hexahedron:
        return GeometryInfo<dim>::d_linear_shape_function_gradient(xi, i);
        // see also BarycentricPolynomials<2>::compute_grad()
      case ReferenceCells::Triangle:
        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);
            default:
              Assert(false, ExcInternalError());
          }
      default:
        Assert(false, ExcNotImplemented());
    }
 
  return Point<dim>(+0.0, +0.0, +0.0);
}
 
 
 
inline double
ReferenceCell::volume() const
{
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        return 0;
      case ReferenceCells::Line:
        return 1;
      case ReferenceCells::Triangle:
        return 1. / 2.;
      case ReferenceCells::Quadrilateral:
        return 1;
      case ReferenceCells::Tetrahedron:
        return 1. / 6.;
      case ReferenceCells::Pyramid:
        return 4. / 3.;
      case ReferenceCells::Wedge:
        return 1. / 2.;
      case ReferenceCells::Hexahedron:
        return 1;
      default:
        Assert(false, ExcNotImplemented());
    }
 
  return 0.0;
}
 
 
 
template <int dim>
inline Point<dim>
ReferenceCell::barycenter() const
{
  AssertDimension(dim, get_dimension());
 
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
        return Point<dim>();
      case ReferenceCells::Line:
        return Point<dim>(1. / 2.);
      case ReferenceCells::Triangle:
        return Point<dim>(1. / 3., 1. / 3.);
      case ReferenceCells::Quadrilateral:
        return Point<dim>(1. / 2., 1. / 2.);
      case ReferenceCells::Tetrahedron:
        return Point<dim>(1. / 4., 1. / 4., 1. / 4.);
      case ReferenceCells::Pyramid:
        return Point<dim>(0, 0, 1. / 4.);
      case ReferenceCells::Wedge:
        return Point<dim>(1. / 3, 1. / 3, 1. / 2.);
      case ReferenceCells::Hexahedron:
        return Point<dim>(1. / 2., 1. / 2., 1. / 2.);
      default:
        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());
 
  // Introduce abbreviations to silence compiler warnings about invalid
  // array accesses (that can't happen, of course, but the compiler
  // doesn't know that).
  constexpr unsigned int x_coordinate = 0;
  constexpr unsigned int y_coordinate = (dim >= 2 ? 1 : 0);
  constexpr unsigned int z_coordinate = (dim >= 3 ? 2 : 0);
 
  switch (this->kind)
    {
      case 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;
        }
      case ReferenceCells::Line:
      case ReferenceCells::Quadrilateral:
      case ReferenceCells::Hexahedron:
        {
          for (unsigned int d = 0; d < dim; ++d)
            if ((p[d] < -tolerance) || (p[d] > 1 + tolerance))
              return false;
          return true;
        }
      case ReferenceCells::Triangle:
      case ReferenceCells::Tetrahedron:
        {
          // 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));
        }
      case ReferenceCells::Pyramid:
        {
          // A pyramid only lives in the upper half-space:
          if (p[z_coordinate] < -tolerance)
            return false;
 
          // It also only lives in the space below z=1:
          if (p[z_coordinate] > 1 + tolerance)
            return false;
 
          // Within what's left of the space, a pyramid is a cone that tapers
          // towards the top. First compute the distance of the point to the
          // axis in the max norm (this is the right norm because the vertices
          // of the pyramid are at points +/-1, +/-1):
          const double distance_from_axis =
            std::max(std::fabs(p[x_coordinate]), std::fabs(p[y_coordinate]));
 
          // We are inside the pyramid if the distance from the axis is less
          // than (1-z)
          return (distance_from_axis <= 1 + tolerance - p[z_coordinate]);
        }
      case ReferenceCells::Wedge:
        {
          // The wedge we use is a triangle extruded into the third
          // dimension by one unit. So we can use the same logic as for
          // triangles above (i.e., for the simplex above, using dim==2)
          // and then check the third dimension separately.
 
          if ((p[x_coordinate] < -tolerance) || (p[y_coordinate] < -tolerance))
            return false;
 
          const double sum = p[x_coordinate] + p[y_coordinate];
          if (sum > 1 + tolerance * std::sqrt(2.0))
            return false;
 
          if (p[z_coordinate] < -tolerance)
            return false;
          if (p[z_coordinate] > 1 + tolerance)
            return false;
 
          return true;
        }
      default:
        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);
 
  switch (this->kind)
    {
      case ReferenceCells::Vertex:
      case ReferenceCells::Line:
      case ReferenceCells::Quadrilateral:
      case ReferenceCells::Hexahedron:
        AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
        return GeometryInfo<dim>::unit_tangential_vectors[face_no][i];
      case 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];
        }
      case 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];
        }
      case 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];
        }
      case 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];
        }
      default:
        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 unsigned int
ReferenceCell::n_face_orientations(const unsigned int face_no) const
{
  AssertIndexRange(face_no, n_faces());
  if (get_dimension() == 1)
    return 1;
  if (get_dimension() == 2)
    return 2;
  else if (face_reference_cell(face_no) == ReferenceCells::Quadrilateral)
    return 8;
  else if (face_reference_cell(face_no) == ReferenceCells::Triangle)
    return 6;
 
  Assert(false, ExcInternalError());
  return numbers::invalid_unsigned_int;
}
 
 
 
inline bool
ReferenceCell::standard_vs_true_line_orientation(
  const unsigned int  line,
  const unsigned char combined_face_orientation,
  const bool          line_orientation) const
{
  if (*this == ReferenceCells::Hexahedron)
    {
      static constexpr ::ndarray<bool, 2, 8> bool_table{
        {{{true, true, false, true, false, false, true, false}},
         {{true, true, true, false, false, false, false, true}}}};
 
      return (line_orientation ==
              bool_table[line / 2][combined_face_orientation]);
    }
  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>
  class NoPermutation : public ::ExceptionBase
  {
  public:
    NoPermutation(const ReferenceCell &     entity_type,
                  const ArrayView<const T> &vertices_0,
                  const ArrayView<const T> &vertices_1)
      : entity_type(entity_type)
      , vertices_0(vertices_0)
      , vertices_1(vertices_1)
    {
      Assert(vertices_0.size() >= entity_type.n_vertices(), ExcInternalError());
      Assert(vertices_1.size() >= entity_type.n_vertices(), ExcInternalError());
    }
 
    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 valid 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 ArrayView<const T> vertices_0;
 
    const ArrayView<const T> 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
{
  Assert(N >= n_vertices(),
         ExcMessage("The number of array elements must be equal to or "
                    "greater than the number of vertices of the cell "
                    "referenced by this object."));
 
  // Call the non-deprecated function, taking care of calling it only with
  // those array elements that we actually care about (see the note
  // in the documentation about the arguments potentially being
  // larger arrays than necessary).
  return get_combined_orientation(
    make_array_view(vertices_0.begin(), vertices_0.begin() + n_vertices()),
    make_array_view(vertices_1.begin(), vertices_1.begin() + n_vertices()));
}
 
 
 
template <typename T>
unsigned char
ReferenceCell::get_combined_orientation(
  const ArrayView<const T> &vertices_0,
  const ArrayView<const T> &vertices_1) const
{
  Assert(vertices_0.size() == n_vertices(),
         ExcMessage("The number of array elements must be equal to "
                    "the number of vertices of the cell "
                    "referenced by this object."));
  Assert(vertices_1.size() == n_vertices(),
         ExcMessage("The number of array elements must be equal to "
                    "the number of vertices of the cell "
                    "referenced by this object."));
 
  const auto v0_equals = [&](const std::initializer_list<const T> &list) {
    Assert(list.size() == n_vertices(), ExcInternalError());
    return std::equal(vertices_0.begin(), vertices_0.end(), std::begin(list));
  };
 
  switch (this->kind)
    {
      case ReferenceCells::Line:
        // line_orientation=true
        if (v0_equals({vertices_1[0], vertices_1[1]}))
          return 1;
 
        // line_orientation=false
        if (v0_equals({vertices_1[1], vertices_1[0]}))
          return 0;
        break;
      case ReferenceCells::Triangle:
        // face_orientation=true, face_rotation=false, face_flip=false
        if (v0_equals({vertices_1[0], vertices_1[1], vertices_1[2]}))
          return 1;
 
        // face_orientation=true, face_rotation=true, face_flip=false
        if (v0_equals({vertices_1[1], vertices_1[2], vertices_1[0]}))
          return 3;
 
        // face_orientation=true, face_rotation=false, face_flip=true
        if (v0_equals({vertices_1[2], vertices_1[0], vertices_1[1]}))
          return 5;
 
        // face_orientation=false, face_rotation=false, face_flip=false
        if (v0_equals({vertices_1[0], vertices_1[2], vertices_1[1]}))
          return 0;
 
        // face_orientation=false, face_rotation=true, face_flip=false
        if (v0_equals({vertices_1[2], vertices_1[1], vertices_1[0]}))
          return 2;
 
        // face_orientation=false, face_rotation=false, face_flip=true
        if (v0_equals({vertices_1[1], vertices_1[0], vertices_1[2]}))
          return 4;
        break;
      case ReferenceCells::Quadrilateral:
        // face_orientation=true, face_rotation=false, face_flip=false
        if (v0_equals(
              {vertices_1[0], vertices_1[1], vertices_1[2], vertices_1[3]}))
          return 1;
 
        // face_orientation=true, face_rotation=true, face_flip=false
        if (v0_equals(
              {vertices_1[2], vertices_1[0], vertices_1[3], vertices_1[1]}))
          return 3;
 
        // face_orientation=true, face_rotation=false, face_flip=true
        if (v0_equals(
              {vertices_1[3], vertices_1[2], vertices_1[1], vertices_1[0]}))
          return 5;
 
        // face_orientation=true, face_rotation=true, face_flip=true
        if (v0_equals(
              {vertices_1[1], vertices_1[3], vertices_1[0], vertices_1[2]}))
          return 7;
 
        // face_orientation=false, face_rotation=false, face_flip=false
        if (v0_equals(
              {vertices_1[0], vertices_1[2], vertices_1[1], vertices_1[3]}))
          return 0;
 
        // face_orientation=false, face_rotation=true, face_flip=false
        if (v0_equals(
              {vertices_1[2], vertices_1[3], vertices_1[0], vertices_1[1]}))
          return 2;
 
        // face_orientation=false, face_rotation=false, face_flip=true
        if (v0_equals(
              {vertices_1[3], vertices_1[1], vertices_1[2], vertices_1[0]}))
          return 4;
 
        // face_orientation=false, face_rotation=true, face_flip=true
        if (v0_equals(
              {vertices_1[1], vertices_1[0], vertices_1[3], vertices_1[2]}))
          return 6;
        break;
      default:
        Assert(false, ExcNotImplemented());
    }
 
  Assert(false, (internal::NoPermutation<T>(*this, vertices_0, vertices_1)));
  return std::numeric_limits<unsigned char>::max();
}
 
 
 
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
{
  Assert(N >= n_vertices(),
         ExcMessage("The number of array elements must be equal to or "
                    "greater than the number of vertices of the cell "
                    "referenced by this object."));
 
  // Call the non-deprecated function, taking care of calling it only with
  // those array elements that we actually care about (see the note
  // in the documentation about the arguments potentially being
  // larger arrays than necessary).
  const auto permutation = permute_by_combined_orientation(
    make_array_view(vertices.begin(), vertices.begin() + n_vertices()),
    orientation);
 
  std::array<T, N> temp;
  std::copy(permutation.begin(), permutation.end(), temp.begin());
 
  return temp;
}
 
 
 
template <typename T>
boost::container::small_vector<T, 8>
ReferenceCell::permute_by_combined_orientation(
  const ArrayView<const T> &vertices,
  const unsigned char       orientation) const
{
  Assert(vertices.size() == n_vertices(),
         ExcMessage("The number of array elements must be equal to "
                    "the number of vertices of the cell "
                    "referenced by this object."));
 
  switch (this->kind)
    {
      case ReferenceCells::Line:
        switch (orientation)
          {
            case 1:
              return {vertices[0], vertices[1]};
            case 0:
              return {vertices[1], vertices[0]};
            default:
              Assert(false, ExcNotImplemented());
          }
        break;
      case ReferenceCells::Triangle:
        switch (orientation)
          {
            case 1:
              return {vertices[0], vertices[1], vertices[2]};
            case 3:
              return {vertices[1], vertices[2], vertices[0]};
            case 5:
              return {vertices[2], vertices[0], vertices[1]};
            case 0:
              return {vertices[0], vertices[2], vertices[1]};
            case 2:
              return {vertices[2], vertices[1], vertices[0]};
            case 4:
              return {vertices[1], vertices[0], vertices[2]};
            default:
              Assert(false, ExcNotImplemented());
          }
        break;
      case ReferenceCells::Quadrilateral:
        switch (orientation)
          {
            case 1:
              return {vertices[0], vertices[1], vertices[2], vertices[3]};
            case 3:
              return {vertices[2], vertices[0], vertices[3], vertices[1]};
            case 5:
              return {vertices[3], vertices[2], vertices[1], vertices[0]};
            case 7:
              return {vertices[1], vertices[3], vertices[0], vertices[2]};
            case 0:
              return {vertices[0], vertices[2], vertices[1], vertices[3]};
            case 2:
              return {vertices[2], vertices[3], vertices[0], vertices[1]};
            case 4:
              return {vertices[3], vertices[1], vertices[2], vertices[0]};
            case 6:
              return {vertices[1], vertices[0], vertices[3], vertices[2]};
            default:
              Assert(false, ExcNotImplemented());
          }
        break;
      default:
        AssertThrow(false, ExcNotImplemented());
    }
 
  return {};
}
 
 
 
template <>
unsigned int
ReferenceCell::vtk_lexicographic_to_node_index<0>(
  const std::array<unsigned, 0> &node_indices,
  const std::array<unsigned, 0> &nodes_per_direction,
  const bool                     legacy_format) const;
 
template <>
unsigned int
ReferenceCell::vtk_lexicographic_to_node_index<1>(
  const std::array<unsigned, 1> &node_indices,
  const std::array<unsigned, 1> &nodes_per_direction,
  const bool                     legacy_format) const;
 
template <>
unsigned int
ReferenceCell::vtk_lexicographic_to_node_index<2>(
  const std::array<unsigned, 2> &node_indices,
  const std::array<unsigned, 2> &nodes_per_direction,
  const bool                     legacy_format) const;
 
template <>
unsigned int
ReferenceCell::vtk_lexicographic_to_node_index<3>(
  const std::array<unsigned, 3> &node_indices,
  const std::array<unsigned, 3> &nodes_per_direction,
  const bool                     legacy_format) const;
 
DEAL_II_NAMESPACE_CLOSE
 
#endif