geometry_info.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2016 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii__geometry_info_h
#define dealii__geometry_info_h


#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/point.h>

DEAL_II_NAMESPACE_OPEN


class GeometryPrimitive
{
public:
  enum Object
  {
    vertex = 0,
    line   = 1,
    quad   = 2,
    hex    = 3
  };

  GeometryPrimitive (const Object object);

  GeometryPrimitive (const unsigned int object_dimension);

  operator unsigned int () const;

private:
  Object object;
};




template <int dim>
struct RefinementPossibilities
{
  enum Possibilities
  {
    no_refinement = 0,

    isotropic_refinement = static_cast<unsigned char>(-1)
  };
};



template <>
struct RefinementPossibilities<1>
{
  enum Possibilities
  {
    no_refinement = 0,
    cut_x        = 1,
    isotropic_refinement = cut_x
  };
};



template <>
struct RefinementPossibilities<2>
{
  enum Possibilities
  {
    no_refinement = 0,
    cut_x         = 1,
    cut_y         = 2,
    cut_xy        = cut_x | cut_y,

    isotropic_refinement = cut_xy
  };
};



template <>
struct RefinementPossibilities<3>
{
  enum Possibilities
  {
    no_refinement = 0,
    cut_x         = 1,
    cut_y         = 2,
    cut_xy        = cut_x | cut_y,
    cut_z         = 4,
    cut_xz        = cut_x | cut_z,
    cut_yz        = cut_y | cut_z,
    cut_xyz       = cut_x | cut_y | cut_z,

    isotropic_refinement = cut_xyz
  };
};



template <int dim>
class RefinementCase : public RefinementPossibilities<dim>
{
public:
  RefinementCase ();

  RefinementCase (const typename RefinementPossibilities<dim>::Possibilities refinement_case);

  explicit RefinementCase (const unsigned char refinement_case);

  operator unsigned char () const;

  RefinementCase operator | (const RefinementCase &r) const;

  RefinementCase operator & (const RefinementCase &r) const;

  RefinementCase operator ~ () const;


  static
  RefinementCase cut_axis (const unsigned int i);

  static std::size_t memory_consumption ();

  template <class Archive>
  void serialize(Archive &ar,
                 const unsigned int version);

  DeclException1 (ExcInvalidRefinementCase,
                  int,
                  << "The refinement flags given (" << arg1 << ") contain set bits that do not "
                  << "make sense for the space dimension of the object to which they are applied.");

private:
unsigned char value :
  (dim > 0 ? dim : 1);
};


namespace internal
{


  template <int dim>
  struct SubfacePossibilities
  {
    enum Possibilities
    {
      case_none = 0,

      case_isotropic = static_cast<unsigned char>(-1)
    };
  };


  template <>
  struct SubfacePossibilities<0>
  {
    enum Possibilities
    {
      case_none = 0,

      case_isotropic = case_none
    };
  };



  template <>
  struct SubfacePossibilities<1>
  {
    enum Possibilities
    {
      case_none = 0,

      case_isotropic = case_none
    };
  };



  template <>
  struct SubfacePossibilities<2>
  {
    enum Possibilities
    {
      case_none = 0,
      case_x    = 1,
      case_isotropic = case_x
    };
  };



  template <>
  struct SubfacePossibilities<3>
  {
    enum Possibilities
    {
      case_none  = 0,
      case_x     = 1,
      case_x1y   = 2,
      case_x2y   = 3,
      case_x1y2y = 4,
      case_y     = 5,
      case_y1x   = 6,
      case_y2x   = 7,
      case_y1x2x = 8,
      case_xy    = 9,

      case_isotropic = case_xy
    };
  };




  template <int dim>
  class SubfaceCase : public SubfacePossibilities<dim>
  {
  public:
    SubfaceCase (const typename SubfacePossibilities<dim>::Possibilities subface_possibility);

    operator unsigned char () const;

    static std::size_t memory_consumption ();

    DeclException1 (ExcInvalidSubfaceCase,
                    int,
                    << "The subface case given (" << arg1 << ") does not make sense "
                    << "for the space dimension of the object to which they are applied.");

  private:
  unsigned char value :
    (dim == 3 ? 4 : 1);
  };

} // namespace internal



template <int dim> struct GeometryInfo;




template <>
struct GeometryInfo<0>
{

  static const unsigned int max_children_per_cell = 1;

  static const unsigned int faces_per_cell    = 0;

  static const unsigned int max_children_per_face = 0;

  static unsigned int n_children(const RefinementCase<0> &refinement_case);

  static const unsigned int vertices_per_cell = 1;

  static const unsigned int vertices_per_face = 0;

  static const unsigned int lines_per_face    = 0;

  static const unsigned int quads_per_face    = 0;

  static const unsigned int lines_per_cell    = 0;

  static const unsigned int quads_per_cell    = 0;

  static const unsigned int hexes_per_cell    = 0;
};





template <int dim>
struct GeometryInfo
{

  static const unsigned int max_children_per_cell = 1 << dim;

  static const unsigned int faces_per_cell = 2 * dim;

  static const unsigned int max_children_per_face = GeometryInfo<dim-1>::max_children_per_cell;

  static const unsigned int vertices_per_cell = 1 << dim;

  static const unsigned int vertices_per_face = GeometryInfo<dim-1>::vertices_per_cell;

  static const unsigned int lines_per_face
    = GeometryInfo<dim-1>::lines_per_cell;

  static const unsigned int quads_per_face
    = GeometryInfo<dim-1>::quads_per_cell;

  static const unsigned int lines_per_cell
    = (2*GeometryInfo<dim-1>::lines_per_cell +
       GeometryInfo<dim-1>::vertices_per_cell);

  static const unsigned int quads_per_cell
    = (2*GeometryInfo<dim-1>::quads_per_cell +
       GeometryInfo<dim-1>::lines_per_cell);

  static const unsigned int hexes_per_cell
    = (2*GeometryInfo<dim-1>::hexes_per_cell +
       GeometryInfo<dim-1>::quads_per_cell);

  static const unsigned int ucd_to_deal[vertices_per_cell];

  static const unsigned int dx_to_deal[vertices_per_cell];

  static const unsigned int vertex_to_face[vertices_per_cell][dim];

  static
  unsigned int
  n_children(const RefinementCase<dim> &refinement_case);

  static
  unsigned int
  n_subfaces(const internal::SubfaceCase<dim> &subface_case);

  static
  double
  subface_ratio(const internal::SubfaceCase<dim> &subface_case,
                const unsigned int subface_no);

  static
  RefinementCase<dim-1>
  face_refinement_case (const RefinementCase<dim> &cell_refinement_case,
                        const unsigned int face_no,
                        const bool face_orientation = true,
                        const bool face_flip        = false,
                        const bool face_rotation    = false);

  static
  RefinementCase<dim>
  min_cell_refinement_case_for_face_refinement
  (const RefinementCase<dim-1> &face_refinement_case,
   const unsigned int face_no,
   const bool face_orientation = true,
   const bool face_flip        = false,
   const bool face_rotation    = false);

  static
  RefinementCase<1>
  line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
                       const unsigned int line_no);

  static
  RefinementCase<dim>
  min_cell_refinement_case_for_line_refinement(const unsigned int line_no);

  static
  unsigned int
  child_cell_on_face (const RefinementCase<dim> &ref_case,
                      const unsigned int face,
                      const unsigned int subface,
                      const bool face_orientation = true,
                      const bool face_flip        = false,
                      const bool face_rotation    = false,
                      const RefinementCase<dim-1> &face_refinement_case
                      = RefinementCase<dim-1>::isotropic_refinement);

  static
  unsigned int
  line_to_cell_vertices (const unsigned int line,
                         const unsigned int vertex);

  static
  unsigned int
  face_to_cell_vertices (const unsigned int face,
                         const unsigned int vertex,
                         const bool face_orientation = true,
                         const bool face_flip        = false,
                         const bool face_rotation    = false);

  static
  unsigned int
  face_to_cell_lines (const unsigned int face,
                      const unsigned int line,
                      const bool face_orientation = true,
                      const bool face_flip        = false,
                      const bool face_rotation    = false);

  static
  unsigned int
  standard_to_real_face_vertex (const unsigned int vertex,
                                const bool face_orientation = true,
                                const bool face_flip        = false,
                                const bool face_rotation    = false);

  static
  unsigned int
  real_to_standard_face_vertex (const unsigned int vertex,
                                const bool face_orientation = true,
                                const bool face_flip        = false,
                                const bool face_rotation    = false);

  static
  unsigned int
  standard_to_real_face_line (const unsigned int line,
                              const bool face_orientation = true,
                              const bool face_flip        = false,
                              const bool face_rotation    = false);

  static
  unsigned int
  real_to_standard_face_line (const unsigned int line,
                              const bool face_orientation = true,
                              const bool face_flip        = false,
                              const bool face_rotation    = false);

  static
  Point<dim>
  unit_cell_vertex (const unsigned int vertex);

  static
  unsigned int
  child_cell_from_point (const Point<dim> &p);

  static
  Point<dim>
  cell_to_child_coordinates (const Point<dim>          &p,
                             const unsigned int         child_index,
                             const RefinementCase<dim>  refine_case
                             = RefinementCase<dim>::isotropic_refinement);

  static
  Point<dim>
  child_to_cell_coordinates (const Point<dim>          &p,
                             const unsigned int         child_index,
                             const RefinementCase<dim>  refine_case
                             = RefinementCase<dim>::isotropic_refinement);

  static
  bool
  is_inside_unit_cell (const Point<dim> &p);

  static
  bool
  is_inside_unit_cell (const Point<dim> &p,
                       const double eps);

  static
  Point<dim>
  project_to_unit_cell (const Point<dim> &p);

  static
  double
  distance_to_unit_cell (const Point<dim> &p);

  static
  double
  d_linear_shape_function (const Point<dim> &xi,
                           const unsigned int i);

  static
  Tensor<1,dim>
  d_linear_shape_function_gradient (const Point<dim> &xi,
                                    const unsigned int i);

  template <int spacedim>
  static void
  alternating_form_at_vertices
#ifndef DEAL_II_CONSTEXPR_BUG
  (const Point<spacedim> (&vertices)[vertices_per_cell],
   Tensor<spacedim-dim,spacedim> (&forms)[vertices_per_cell]);
#else
  (const Point<spacedim> *vertices,
   Tensor<spacedim-dim,spacedim> *forms);
#endif

  static const unsigned int unit_normal_direction[faces_per_cell];

  static const int unit_normal_orientation[faces_per_cell];

  static const unsigned int opposite_face[faces_per_cell];


  DeclException1 (ExcInvalidCoordinate,
                  double,
                  << "The coordinates must satisfy 0 <= x_i <= 1, "
                  << "but here we have x_i=" << arg1);

  DeclException3 (ExcInvalidSubface,
                  int, int, int,
                  << "RefinementCase<dim> " << arg1 << ": face " << arg2
                  << " has no subface " << arg3);
};




#ifndef DOXYGEN


/* -------------- declaration of explicit specializations ------------- */

#ifndef DEAL_II_MEMBER_ARRAY_SPECIALIZATION_BUG
template <>
const unsigned int GeometryInfo<1>::unit_normal_direction[faces_per_cell];
template <>
const unsigned int GeometryInfo<2>::unit_normal_direction[faces_per_cell];
template <>
const unsigned int GeometryInfo<3>::unit_normal_direction[faces_per_cell];
template <>
const unsigned int GeometryInfo<4>::unit_normal_direction[faces_per_cell];

template <>
const int GeometryInfo<1>::unit_normal_orientation[faces_per_cell];
template <>
const int GeometryInfo<2>::unit_normal_orientation[faces_per_cell];
template <>
const int GeometryInfo<3>::unit_normal_orientation[faces_per_cell];
template <>
const int GeometryInfo<4>::unit_normal_orientation[faces_per_cell];

template <>
const unsigned int GeometryInfo<1>::opposite_face[faces_per_cell];
template <>
const unsigned int GeometryInfo<2>::opposite_face[faces_per_cell];
template <>
const unsigned int GeometryInfo<3>::opposite_face[faces_per_cell];
template <>
const unsigned int GeometryInfo<4>::opposite_face[faces_per_cell];
#endif


template <>
Tensor<1,1>
GeometryInfo<1>::
d_linear_shape_function_gradient (const Point<1> &xi,
                                  const unsigned int i);
template <>
Tensor<1,2>
GeometryInfo<2>::
d_linear_shape_function_gradient (const Point<2> &xi,
                                  const unsigned int i);
template <>
Tensor<1,3>
GeometryInfo<3>::
d_linear_shape_function_gradient (const Point<3> &xi,
                                  const unsigned int i);




/* -------------- inline functions ------------- */


inline
GeometryPrimitive::GeometryPrimitive (const Object object)
  :
  object (object)
{}



inline
GeometryPrimitive::GeometryPrimitive (const unsigned int object_dimension)
  :
  object (static_cast<Object>(object_dimension))
{}


inline
GeometryPrimitive::operator unsigned int () const
{
  return static_cast<unsigned int>(object);
}



namespace internal
{

  template <int dim>
  inline
  SubfaceCase<dim>::SubfaceCase (const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
    :
    value (subface_possibility)
  {}


  template <int dim>
  inline
  SubfaceCase<dim>::operator unsigned char () const
  {
    return value;
  }


} // namespace internal


template <int dim>
inline
RefinementCase<dim>
RefinementCase<dim>::cut_axis (const unsigned int)
{
  Assert (false, ExcInternalError());
  return static_cast<unsigned char>(-1);
}


template <>
inline
RefinementCase<1>
RefinementCase<1>::cut_axis (const unsigned int i)
{
  const unsigned int dim = 1;
  Assert (i < dim, ExcIndexRange(i, 0, dim));

  static const RefinementCase options[dim] = { RefinementPossibilities<1>::cut_x };
  return options[i];
}



template <>
inline
RefinementCase<2>
RefinementCase<2>::cut_axis (const unsigned int i)
{
  const unsigned int dim = 2;
  Assert (i < dim, ExcIndexRange(i, 0, dim));

  static const RefinementCase options[dim] = { RefinementPossibilities<2>::cut_x,
                                               RefinementPossibilities<2>::cut_y
                                             };
  return options[i];
}



template <>
inline
RefinementCase<3>
RefinementCase<3>::cut_axis (const unsigned int i)
{
  const unsigned int dim = 3;
  Assert (i < dim, ExcIndexRange(i, 0, dim));

  static const RefinementCase options[dim] = { RefinementPossibilities<3>::cut_x,
                                               RefinementPossibilities<3>::cut_y,
                                               RefinementPossibilities<3>::cut_z
                                             };
  return options[i];
}



template <int dim>
inline
RefinementCase<dim>::RefinementCase ()
  :
  value (RefinementPossibilities<dim>::no_refinement)
{}



template <int dim>
inline
RefinementCase<dim>::
RefinementCase (const typename RefinementPossibilities<dim>::Possibilities refinement_case)
  :
  value (refinement_case)
{
  // check that only those bits of
  // the given argument are set that
  // make sense for a given space
  // dimension
  Assert ((refinement_case & RefinementPossibilities<dim>::isotropic_refinement) ==
          refinement_case,
          ExcInvalidRefinementCase (refinement_case));
}



template <int dim>
inline
RefinementCase<dim>::RefinementCase (const unsigned char refinement_case)
  :
  value (refinement_case)
{
  // check that only those bits of
  // the given argument are set that
  // make sense for a given space
  // dimension
  Assert ((refinement_case & RefinementPossibilities<dim>::isotropic_refinement) ==
          refinement_case,
          ExcInvalidRefinementCase (refinement_case));
}



template <int dim>
inline
RefinementCase<dim>::operator unsigned char () const
{
  return value;
}



template <int dim>
inline
RefinementCase<dim>
RefinementCase<dim>::operator | (const RefinementCase<dim> &r) const
{
  return RefinementCase<dim>(static_cast<unsigned char> (value | r.value));
}



template <int dim>
inline
RefinementCase<dim>
RefinementCase<dim>::operator & (const RefinementCase<dim> &r) const
{
  return RefinementCase<dim>(static_cast<unsigned char> (value & r.value));
}



template <int dim>
inline
RefinementCase<dim>
RefinementCase<dim>::operator ~ () const
{
  return RefinementCase<dim>(static_cast<unsigned char> (
                               (~value) & RefinementPossibilities<dim>::isotropic_refinement));
}




template <int dim>
inline
std::size_t
RefinementCase<dim>::memory_consumption ()
{
  return sizeof(RefinementCase<dim>);
}



template <int dim>
template <class Archive>
void RefinementCase<dim>::serialize (Archive &ar,
                                     const unsigned int)
{
  // serialization can't deal with bitfields, so copy from/to a full sized
  // unsigned char
  unsigned char uchar_value = value;
  ar &uchar_value;
  value = uchar_value;
}




template <>
inline
Point<1>
GeometryInfo<1>::unit_cell_vertex (const unsigned int vertex)
{
  Assert (vertex < vertices_per_cell,
          ExcIndexRange (vertex, 0, vertices_per_cell));

  return Point<1>(static_cast<double>(vertex));
}



template <>
inline
Point<2>
GeometryInfo<2>::unit_cell_vertex (const unsigned int vertex)
{
  Assert (vertex < vertices_per_cell,
          ExcIndexRange (vertex, 0, vertices_per_cell));

  return Point<2>(vertex%2, vertex/2);
}



template <>
inline
Point<3>
GeometryInfo<3>::unit_cell_vertex (const unsigned int vertex)
{
  Assert (vertex < vertices_per_cell,
          ExcIndexRange (vertex, 0, vertices_per_cell));

  return Point<3>(vertex%2, vertex/2%2, vertex/4);
}



template <int dim>
inline
Point<dim>
GeometryInfo<dim>::unit_cell_vertex (const unsigned int)
{
  Assert(false, ExcNotImplemented());

  return Point<dim> ();
}



template <>
inline
unsigned int
GeometryInfo<1>::child_cell_from_point (const Point<1> &p)
{
  Assert ((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));

  return (p[0] <= 0.5 ? 0 : 1);
}



template <>
inline
unsigned int
GeometryInfo<2>::child_cell_from_point (const Point<2> &p)
{
  Assert ((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
  Assert ((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));

  return (p[0] <= 0.5 ?
          (p[1] <= 0.5 ? 0 : 2) :
          (p[1] <= 0.5 ? 1 : 3));
}



template <>
inline
unsigned int
GeometryInfo<3>::child_cell_from_point (const Point<3> &p)
{
  Assert ((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
  Assert ((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
  Assert ((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));

  return (p[0] <= 0.5 ?
          (p[1] <= 0.5 ?
           (p[2] <= 0.5 ? 0 : 4) :
           (p[2] <= 0.5 ? 2 : 6)) :
          (p[1] <= 0.5 ?
           (p[2] <= 0.5 ? 1 : 5) :
           (p[2] <= 0.5 ? 3 : 7)));
}


template <int dim>
inline
unsigned int
GeometryInfo<dim>::child_cell_from_point (const Point<dim> &)
{
  Assert(false, ExcNotImplemented());

  return 0;
}



template <>
inline
Point<1>
GeometryInfo<1>::cell_to_child_coordinates (const Point<1>         &p,
                                            const unsigned int      child_index,
                                            const RefinementCase<1> refine_case)

{
  Assert (child_index < 2,
          ExcIndexRange (child_index, 0, 2));
  Assert (refine_case==RefinementCase<1>::cut_x,
          ExcInternalError());
  (void)refine_case; // removes -Wunused-parameter warning in optimized mode

  return Point<1>(p*2.0-unit_cell_vertex(child_index));
}



template <>
inline
Point<2>
GeometryInfo<2>::cell_to_child_coordinates (const Point<2>         &p,
                                            const unsigned int      child_index,
                                            const RefinementCase<2> refine_case)

{
  Assert (child_index < GeometryInfo<2>::n_children(refine_case),
          ExcIndexRange (child_index, 0, GeometryInfo<2>::n_children(refine_case)));

  Point<2> point=p;
  switch (refine_case)
    {
    case RefinementCase<2>::cut_x:
      point[0]*=2.0;
      if (child_index==1)
        point[0]-=1.0;
      break;
    case RefinementCase<2>::cut_y:
      point[1]*=2.0;
      if (child_index==1)
        point[1]-=1.0;
      break;
    case RefinementCase<2>::cut_xy:
      point*=2.0;
      point-=unit_cell_vertex(child_index);
      break;
    default:
      Assert(false, ExcInternalError());
    }

  return point;
}



template <>
inline
Point<3>
GeometryInfo<3>::cell_to_child_coordinates (const Point<3>         &p,
                                            const unsigned int      child_index,
                                            const RefinementCase<3> refine_case)

{
  Assert (child_index < GeometryInfo<3>::n_children(refine_case),
          ExcIndexRange (child_index, 0, GeometryInfo<3>::n_children(refine_case)));

  Point<3> point=p;
  // there might be a cleverer way to do
  // this, but since this function is called
  // in very few places for initialization
  // purposes only, I don't care at the
  // moment
  switch (refine_case)
    {
    case RefinementCase<3>::cut_x:
      point[0]*=2.0;
      if (child_index==1)
        point[0]-=1.0;
      break;
    case RefinementCase<3>::cut_y:
      point[1]*=2.0;
      if (child_index==1)
        point[1]-=1.0;
      break;
    case RefinementCase<3>::cut_z:
      point[2]*=2.0;
      if (child_index==1)
        point[2]-=1.0;
      break;
    case RefinementCase<3>::cut_xy:
      point[0]*=2.0;
      point[1]*=2.0;
      if (child_index%2==1)
        point[0]-=1.0;
      if (child_index/2==1)
        point[1]-=1.0;
      break;
    case RefinementCase<3>::cut_xz:
      // careful, this is slightly
      // different from xy and yz due to
      // different internal numbering of
      // children!
      point[0]*=2.0;
      point[2]*=2.0;
      if (child_index/2==1)
        point[0]-=1.0;
      if (child_index%2==1)
        point[2]-=1.0;
      break;
    case RefinementCase<3>::cut_yz:
      point[1]*=2.0;
      point[2]*=2.0;
      if (child_index%2==1)
        point[1]-=1.0;
      if (child_index/2==1)
        point[2]-=1.0;
      break;
    case RefinementCase<3>::cut_xyz:
      point*=2.0;
      point-=unit_cell_vertex(child_index);
      break;
    default:
      Assert(false, ExcInternalError());
    }

  return point;
}



template <int dim>
inline
Point<dim>
GeometryInfo<dim>::cell_to_child_coordinates (const Point<dim>         &/*p*/,
                                              const unsigned int        /*child_index*/,
                                              const RefinementCase<dim> /*refine_case*/)

{
  AssertThrow (false, ExcNotImplemented());
  return Point<dim>();
}



template <>
inline
Point<1>
GeometryInfo<1>::child_to_cell_coordinates (const Point<1>         &p,
                                            const unsigned int      child_index,
                                            const RefinementCase<1> refine_case)

{
  Assert (child_index < 2,
          ExcIndexRange (child_index, 0, 2));
  Assert (refine_case==RefinementCase<1>::cut_x,
          ExcInternalError());
  (void)refine_case; // removes -Wunused-parameter warning in optimized mode

  return (p+unit_cell_vertex(child_index))*0.5;
}



template <>
inline
Point<3>
GeometryInfo<3>::child_to_cell_coordinates (const Point<3>         &p,
                                            const unsigned int      child_index,
                                            const RefinementCase<3> refine_case)

{
  Assert (child_index < GeometryInfo<3>::n_children(refine_case),
          ExcIndexRange (child_index, 0, GeometryInfo<3>::n_children(refine_case)));

  Point<3> point=p;
  // there might be a cleverer way to do
  // this, but since this function is called
  // in very few places for initialization
  // purposes only, I don't care at the
  // moment
  switch (refine_case)
    {
    case RefinementCase<3>::cut_x:
      if (child_index==1)
        point[0]+=1.0;
      point[0]*=0.5;
      break;
    case RefinementCase<3>::cut_y:
      if (child_index==1)
        point[1]+=1.0;
      point[1]*=0.5;
      break;
    case RefinementCase<3>::cut_z:
      if (child_index==1)
        point[2]+=1.0;
      point[2]*=0.5;
      break;
    case RefinementCase<3>::cut_xy:
      if (child_index%2==1)
        point[0]+=1.0;
      if (child_index/2==1)
        point[1]+=1.0;
      point[0]*=0.5;
      point[1]*=0.5;
      break;
    case RefinementCase<3>::cut_xz:
      // careful, this is slightly
      // different from xy and yz due to
      // different internal numbering of
      // children!
      if (child_index/2==1)
        point[0]+=1.0;
      if (child_index%2==1)
        point[2]+=1.0;
      point[0]*=0.5;
      point[2]*=0.5;
      break;
    case RefinementCase<3>::cut_yz:
      if (child_index%2==1)
        point[1]+=1.0;
      if (child_index/2==1)
        point[2]+=1.0;
      point[1]*=0.5;
      point[2]*=0.5;
      break;
    case RefinementCase<3>::cut_xyz:
      point+=unit_cell_vertex(child_index);
      point*=0.5;
      break;
    default:
      Assert(false, ExcInternalError());
    }

  return point;
}



template <>
inline
Point<2>
GeometryInfo<2>::child_to_cell_coordinates (const Point<2>         &p,
                                            const unsigned int      child_index,
                                            const RefinementCase<2> refine_case)
{
  Assert (child_index < GeometryInfo<2>::n_children(refine_case),
          ExcIndexRange (child_index, 0, GeometryInfo<2>::n_children(refine_case)));

  Point<2> point=p;
  switch (refine_case)
    {
    case RefinementCase<2>::cut_x:
      if (child_index==1)
        point[0]+=1.0;
      point[0]*=0.5;
      break;
    case RefinementCase<2>::cut_y:
      if (child_index==1)
        point[1]+=1.0;
      point[1]*=0.5;
      break;
    case RefinementCase<2>::cut_xy:
      point+=unit_cell_vertex(child_index);
      point*=0.5;
      break;
    default:
      Assert(false, ExcInternalError());
    }

  return point;
}



template <int dim>
inline
Point<dim>
GeometryInfo<dim>::child_to_cell_coordinates (const Point<dim>         &/*p*/,
                                              const unsigned int        /*child_index*/,
                                              const RefinementCase<dim> /*refine_case*/)
{
  AssertThrow (false, ExcNotImplemented());
  return Point<dim>();
}



template <>
inline
bool
GeometryInfo<1>::is_inside_unit_cell (const Point<1> &p)
{
  return (p[0] >= 0.) && (p[0] <= 1.);
}



template <>
inline
bool
GeometryInfo<2>::is_inside_unit_cell (const Point<2> &p)
{
  return (p[0] >= 0.) && (p[0] <= 1.) &&
         (p[1] >= 0.) && (p[1] <= 1.);
}



template <>
inline
bool
GeometryInfo<3>::is_inside_unit_cell (const Point<3> &p)
{
  return (p[0] >= 0.) && (p[0] <= 1.) &&
         (p[1] >= 0.) && (p[1] <= 1.) &&
         (p[2] >= 0.) && (p[2] <= 1.);
}

template <>
inline
bool
GeometryInfo<1>::is_inside_unit_cell (const Point<1> &p,
                                      const double eps)
{
  return (p[0] >= -eps) && (p[0] <= 1.+eps);
}



template <>
inline
bool
GeometryInfo<2>::is_inside_unit_cell (const Point<2> &p,
                                      const double eps)
{
  const double l = -eps, u = 1+eps;
  return (p[0] >= l) && (p[0] <= u) &&
         (p[1] >= l) && (p[1] <= u);
}



template <>
inline
bool
GeometryInfo<3>::is_inside_unit_cell (const Point<3> &p,
                                      const double eps)
{
  const double l = -eps, u = 1.0+eps;
  return (p[0] >= l) && (p[0] <= u) &&
         (p[1] >= l) && (p[1] <= u) &&
         (p[2] >= l) && (p[2] <= u);
}


#endif // DOXYGEN
DEAL_II_NAMESPACE_CLOSE

#endif