tria_accessor.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__tria_accessor_h
#define dealii__tria_accessor_h


#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/point.h>
#include <deal.II/grid/tria_iterator_base.h>
#include <deal.II/grid/tria_iterator_selector.h>
#include <deal.II/grid/cell_id.h>

#include <utility>


DEAL_II_NAMESPACE_OPEN

template <int dim, int spacedim> class Triangulation;
template <typename Accessor> class TriaRawIterator;
template <typename Accessor> class TriaIterator;
template <typename Accessor> class TriaActiveIterator;

template <int dim, int spacedim> class Boundary;
template <int dim, int spacedim> class Manifold;


namespace internal
{
  namespace Triangulation
  {
    template <int dim> class TriaObject;
    template <typename G> class TriaObjects;
    struct Implementation;
  }

  namespace TriaAccessor
  {
    struct Implementation;

    template <int structdim, int dim> struct PresentLevelType
    {
      struct type
      {
        type ()
        {}

        type (const int level)
        {
          Assert (level == 0, ExcInternalError());
          (void)level; // removes -Wunused-parameter warning in optimized mode
        }

        operator int () const
        {
          return 0;
        }

        void operator ++ () const
        {
          Assert (false, ExcInternalError());
        }

        void operator -- () const
        {
          Assert (false, ExcInternalError());
        }
      };
    };


    template <int dim> struct PresentLevelType<dim,dim>
    {
      typedef int type;
    };

  }
}
template <int structdim, int dim, int spacedim> class TriaAccessor;
template <int dim, int spacedim>                class TriaAccessor<0, dim, spacedim>;
template <int spacedim>                         class TriaAccessor<0, 1, spacedim>;

// note: the file tria_accessor.templates.h is included at the end of
// this file.  this includes a lot of templates. originally, this was
// only done in debug mode, but led to cyclic reduction problems and
// so is now on by default.


namespace TriaAccessorExceptions
{
//TODO: Write documentation!
  DeclException0 (ExcCellNotUsed);
  DeclException0 (ExcCellNotActive);
  DeclException0 (ExcCellHasNoChildren);
  DeclException0 (ExcCellHasNoParent);
//TODO: Write documentation!
  DeclException0 (ExcUnusedCellAsChild);
//TODO: Write documentation!
  DeclException1 (ExcCantSetChildren,
                  int,
                  << "You can only set the child index if the cell has no "
                  << "children, or clear it. The given "
                  << "index was " << arg1 << " (-1 means: clear children)");
//TODO: Write documentation!
  DeclException0 (ExcUnusedCellAsNeighbor);
//TODO: Write documentation!
  DeclException0 (ExcUncaughtCase);
//TODO: Write documentation!
  DeclException0 (ExcDereferenceInvalidObject);
//TODO: Write documentation!
  DeclException0 (ExcCantCompareIterators);
//TODO: Write documentation!
  DeclException0 (ExcNeighborIsCoarser);
//TODO: Write documentation!
  DeclException0 (ExcNeighborIsNotCoarser);
  DeclException0 (ExcFacesHaveNoLevel);
  DeclException0 (ExcNoPeriodicNeighbor);
//TODO: Write documentation!
  DeclException1 (ExcSetOnlyEvenChildren,
                  int,
                  << "You can only set the child index of an even numbered child."
                  << "The number of the child given was " << arg1 << ".");
}


template <int structdim, int dim, int spacedim=dim>
class TriaAccessorBase
{
public:
  static const unsigned int space_dimension = spacedim;

  static const unsigned int dimension = dim;

  static const unsigned int structure_dimension = structdim;

protected:
  typedef void AccessorData;

  TriaAccessorBase (const Triangulation<dim,spacedim> *parent =  0,
                    const int                 level  = -1,
                    const int                 index  = -1,
                    const AccessorData             * =  0);

  TriaAccessorBase (const TriaAccessorBase &);

  void copy_from (const TriaAccessorBase &);

  TriaAccessorBase &operator = (const TriaAccessorBase &);

  bool operator < (const TriaAccessorBase &other) const;

protected:
  void operator = (const TriaAccessorBase *);

  bool operator == (const TriaAccessorBase &) const;

  bool operator != (const TriaAccessorBase &) const;

  void operator ++ ();

  void operator -- ();
  ::internal::Triangulation::TriaObjects<::internal::Triangulation::TriaObject<structdim> > &
  objects () const;

public:
  typedef void *LocalData;

  int level () const;

  int index () const;

  IteratorState::IteratorStates state () const;

  const Triangulation<dim,spacedim> &get_triangulation () const;

protected:
  typename ::internal::TriaAccessor::PresentLevelType<structdim,dim>::type present_level;

  int present_index;

  const Triangulation<dim,spacedim> *tria;

private:

  template <typename Accessor> friend class TriaRawIterator;
  template <typename Accessor> friend class TriaIterator;
  template <typename Accessor> friend class TriaActiveIterator;
};



template <int structdim, int dim, int spacedim=dim>
class InvalidAccessor :  public TriaAccessorBase<structdim,dim,spacedim>
{
public:
  typedef typename TriaAccessorBase<structdim,dim,spacedim>::AccessorData AccessorData;

  InvalidAccessor (const Triangulation<dim,spacedim> *parent     =  0,
                   const int                 level      = -1,
                   const int                 index      = -1,
                   const AccessorData       *local_data =  0);

  InvalidAccessor (const InvalidAccessor &);

  template <typename OtherAccessor>
  InvalidAccessor (const OtherAccessor &);

  void copy_from (const InvalidAccessor &);

  bool operator == (const InvalidAccessor &) const;
  bool operator != (const InvalidAccessor &) const;

  void operator ++ () const;
  void operator -- () const;

  bool used () const;

  bool has_children () const;
};



template <int structdim, int dim, int spacedim>
class TriaAccessor : public TriaAccessorBase<structdim, dim, spacedim>
{
public:
  typedef
  typename TriaAccessorBase<structdim,dim,spacedim>::AccessorData
  AccessorData;

  TriaAccessor (const Triangulation<dim,spacedim> *parent     =  0,
                const int                 level      = -1,
                const int                 index      = -1,
                const AccessorData       *local_data =  0);

  template <int structdim2, int dim2, int spacedim2>
  TriaAccessor (const InvalidAccessor<structdim2,dim2,spacedim2> &);

  template <int structdim2, int dim2, int spacedim2>
  TriaAccessor (const TriaAccessor<structdim2,dim2,spacedim2> &);

  bool used () const;

  typename ::internal::Triangulation::Iterators<dim,spacedim>::vertex_iterator
  vertex_iterator (const unsigned int i) const;

  unsigned int vertex_index (const unsigned int i) const;

  Point<spacedim> &vertex (const unsigned int i) const;

  typename ::internal::Triangulation::Iterators<dim,spacedim>::line_iterator
  line (const unsigned int i) const;

  unsigned int line_index (const unsigned int i) const;

  typename ::internal::Triangulation::Iterators<dim,spacedim>::quad_iterator
  quad (const unsigned int i) const;

  unsigned int quad_index (const unsigned int i) const;
  bool face_orientation (const unsigned int face) const;

  bool face_flip (const unsigned int face) const;

  bool face_rotation (const unsigned int face) const;

  bool line_orientation (const unsigned int line) const;
  bool has_children () const;

  unsigned int n_children() const;

  unsigned int number_of_children () const;

  unsigned int max_refinement_depth () const;

  TriaIterator<TriaAccessor<structdim,dim,spacedim> >
  child (const unsigned int i) const;

  TriaIterator<TriaAccessor<structdim,dim,spacedim> >
  isotropic_child (const unsigned int i) const;

  RefinementCase<structdim> refinement_case () const;

  int child_index (const unsigned int i) const;

  int isotropic_child_index (const unsigned int i) const;
  types::boundary_id boundary_id () const;

  void set_boundary_id (const types::boundary_id) const;

  void set_all_boundary_ids (const types::boundary_id) const;

  bool at_boundary () const;

  const Boundary<dim,spacedim> &get_boundary () const;

  const Manifold<dim,spacedim> &get_manifold () const;

  types::manifold_id manifold_id () const;

  void set_manifold_id (const types::manifold_id) const;

  void set_all_manifold_ids (const types::manifold_id) const;

  bool user_flag_set () const;

  void set_user_flag () const;

  void clear_user_flag () const;

  void recursively_set_user_flag () const;

  void recursively_clear_user_flag () const;

  void clear_user_data () const;

  void set_user_pointer (void *p) const;

  void clear_user_pointer () const;

  void *user_pointer () const;

  void recursively_set_user_pointer (void *p) const;

  void recursively_clear_user_pointer () const;

  void set_user_index (const unsigned int p) const;

  void clear_user_index () const;

  unsigned int user_index () const;

  void recursively_set_user_index (const unsigned int p) const;

  void recursively_clear_user_index () const;
  double diameter () const;

  double extent_in_direction (const unsigned int axis) const;

  double minimum_vertex_distance () const;

  Point<spacedim> intermediate_point(const Point<structdim> &coordinates) const;

  Point<spacedim> center (const bool respect_manifold=false,
                          const bool use_laplace_transformation=false) const;

  Point<spacedim> barycenter () const;

  double measure () const;

  bool
  is_translation_of (const TriaIterator<TriaAccessor<structdim,dim,spacedim> > &o) const;

private:
  void set_boundary_id_internal(const types::boundary_id id) const;

  void set (const ::internal::Triangulation::TriaObject<structdim> &o) const;

  void set_line_orientation (const unsigned int line,
                             const bool         orientation) const;

  void set_face_orientation (const unsigned int face,
                             const bool         orientation) const;

  void set_face_flip (const unsigned int face,
                      const bool         flip) const;

  void set_face_rotation (const unsigned int face,
                          const bool         rotation) const;

  void set_used_flag () const;

  void clear_used_flag () const;

  void set_refinement_case (const RefinementCase<structdim> &ref_case) const;

  void clear_refinement_case () const;

  void set_children (const unsigned int i, const int index) const;

  void clear_children () const;

private:
  void operator = (const TriaAccessor &);

  template <int, int> friend class Triangulation;

  friend struct ::internal::Triangulation::Implementation;
  friend struct ::internal::TriaAccessor::Implementation;
};






template<int dim, int spacedim>
class TriaAccessor<0, dim, spacedim>
{
public:
  static const unsigned int space_dimension = spacedim;

  static const unsigned int dimension = dim;

  static const unsigned int structure_dimension = 0;

  typedef void AccessorData;

  TriaAccessor (const Triangulation<dim,spacedim> *tria,
                const unsigned int    vertex_index);

  TriaAccessor (const Triangulation<dim,spacedim> *tria  = NULL,
                const int                          level = 0,
                const int                          index = 0,
                const AccessorData                     * = 0);

  template <int structdim2, int dim2, int spacedim2>
  TriaAccessor (const TriaAccessor<structdim2,dim2,spacedim2> &);

  template <int structdim2, int dim2, int spacedim2>
  TriaAccessor (const InvalidAccessor<structdim2,dim2,spacedim2> &);

  IteratorState::IteratorStates state () const;

  static int level ();

  int index () const;

  void operator ++ ();

  void operator -- ();
  bool operator == (const TriaAccessor &) const;

  bool operator != (const TriaAccessor &) const;

  unsigned int vertex_index (const unsigned int i = 0) const;

  Point<spacedim> &vertex (const unsigned int i = 0) const;

  typename ::internal::Triangulation::Iterators<dim,spacedim>::line_iterator
  static line (const unsigned int);

  static unsigned int line_index (const unsigned int i);

  static
  typename ::internal::Triangulation::Iterators<dim,spacedim>::quad_iterator
  quad (const unsigned int i);

  static unsigned int quad_index (const unsigned int i);

  double diameter () const;

  double extent_in_direction (const unsigned int axis) const;

  Point<spacedim> center (const bool respect_manifold=false,
                          const bool use_laplace_transformation=false) const;

  double measure () const;
  static bool face_orientation (const unsigned int face);

  static bool face_flip (const unsigned int face);

  static bool face_rotation (const unsigned int face);

  static bool line_orientation (const unsigned int line);

  static bool has_children ();

  static unsigned int n_children();

  static unsigned int number_of_children ();

  static unsigned int max_refinement_depth ();

  static
  TriaIterator<TriaAccessor<0,dim,spacedim> >
  child (const unsigned int);

  static
  TriaIterator<TriaAccessor<0,dim,spacedim> >
  isotropic_child (const unsigned int);

  static
  RefinementCase<0> refinement_case ();

  static
  int child_index (const unsigned int i);

  static
  int isotropic_child_index (const unsigned int i);
  bool used () const;

protected:
  void copy_from (const TriaAccessor &);

  const Triangulation<dim,spacedim> *tria;

  unsigned int  global_vertex_index;

private:

  template <typename Accessor> friend class TriaRawIterator;
  template <typename Accessor> friend class TriaIterator;
  template <typename Accessor> friend class TriaActiveIterator;
};



template <int spacedim>
class TriaAccessor<0, 1, spacedim>
{
public:
  static const unsigned int space_dimension = spacedim;

  static const unsigned int dimension = 1;

  static const unsigned int structure_dimension = 0;

  typedef void AccessorData;

  enum VertexKind
  {
    left_vertex,
    interior_vertex,
    right_vertex
  };

  TriaAccessor (const Triangulation<1,spacedim> *tria,
                const VertexKind      vertex_kind,
                const unsigned int    vertex_index);

  TriaAccessor (const Triangulation<1,spacedim> *tria = 0,
                const int = 0,
                const int = 0,
                const AccessorData * = 0);

  template <int structdim2, int dim2, int spacedim2>
  TriaAccessor (const TriaAccessor<structdim2,dim2,spacedim2> &);

  template <int structdim2, int dim2, int spacedim2>
  TriaAccessor (const InvalidAccessor<structdim2,dim2,spacedim2> &);

  void copy_from (const TriaAccessor &);

  static IteratorState::IteratorStates state ();

  static int level ();

  int index () const;

  void operator ++ () const;

  void operator -- () const;
  bool operator == (const TriaAccessor &) const;

  bool operator != (const TriaAccessor &) const;

  unsigned int vertex_index (const unsigned int i = 0) const;

  Point<spacedim> &vertex (const unsigned int i = 0) const;

  Point<spacedim> center () const;

  typename ::internal::Triangulation::Iterators<1,spacedim>::line_iterator
  static line (const unsigned int);

  static unsigned int line_index (const unsigned int i);

  static
  typename ::internal::Triangulation::Iterators<1,spacedim>::quad_iterator
  quad (const unsigned int i);

  static unsigned int quad_index (const unsigned int i);

  bool at_boundary () const;

  types::boundary_id boundary_id () const;

  types::manifold_id manifold_id () const;


  static bool face_orientation (const unsigned int face);

  static bool face_flip (const unsigned int face);

  static bool face_rotation (const unsigned int face);

  static bool line_orientation (const unsigned int line);

  static bool has_children ();

  static unsigned int n_children();

  static unsigned int number_of_children ();

  static unsigned int max_refinement_depth ();

  static
  TriaIterator<TriaAccessor<0,1,spacedim> >
  child (const unsigned int);

  static
  TriaIterator<TriaAccessor<0,1,spacedim> >
  isotropic_child (const unsigned int);

  static
  RefinementCase<0> refinement_case ();

  static
  int child_index (const unsigned int i);

  static
  int isotropic_child_index (const unsigned int i);
  void
  set_boundary_id (const types::boundary_id);

  void
  set_manifold_id (const types::manifold_id);

  void
  set_all_boundary_ids (const types::boundary_id);

  void
  set_all_manifold_ids (const types::manifold_id);
  bool used () const;

protected:
  const Triangulation<1,spacedim> *tria;

  VertexKind      vertex_kind;

  unsigned int    global_vertex_index;
};




template <int dim, int spacedim=dim>
class CellAccessor :  public TriaAccessor<dim,dim,spacedim>
{
public:
  typedef typename TriaAccessor<dim,dim,spacedim>::AccessorData AccessorData;

  typedef Triangulation<dim, spacedim> Container;

  CellAccessor (const Triangulation<dim,spacedim> *parent     =  0,
                const int                 level      = -1,
                const int                 index      = -1,
                const AccessorData       *local_data =  0);

  CellAccessor (const TriaAccessor<dim,dim,spacedim> &cell_accessor);

  template <int structdim2, int dim2, int spacedim2>
  CellAccessor (const InvalidAccessor<structdim2,dim2,spacedim2> &);

  template <int structdim2, int dim2, int spacedim2>
  CellAccessor (const TriaAccessor<structdim2,dim2,spacedim2> &);

  TriaIterator<CellAccessor<dim, spacedim> >
  child (const unsigned int i) const;

  TriaIterator<TriaAccessor<dim-1,dim,spacedim> >
  face (const unsigned int i) const;

  unsigned int
  face_index (const unsigned int i) const;

  TriaIterator<CellAccessor<dim, spacedim> >
  neighbor_child_on_subface (const unsigned int face_no,
                             const unsigned int subface_no) const;

  TriaIterator<CellAccessor<dim, spacedim>  >
  neighbor (const unsigned int i) const;

  int neighbor_index (const unsigned int i) const;

  int neighbor_level (const unsigned int i) const;

  unsigned int neighbor_of_neighbor (const unsigned int neighbor) const;

  bool neighbor_is_coarser (const unsigned int neighbor) const;

  std::pair<unsigned int, unsigned int>
  neighbor_of_coarser_neighbor (const unsigned int neighbor) const;

  unsigned int neighbor_face_no (const unsigned int neighbor) const;

  bool has_periodic_neighbor(const unsigned int i) const;

  TriaIterator<CellAccessor<dim, spacedim> >
  periodic_neighbor (const unsigned int i) const;

  TriaIterator<CellAccessor<dim, spacedim> >
  neighbor_or_periodic_neighbor (const unsigned int i) const;

  TriaIterator<CellAccessor<dim, spacedim> >
  periodic_neighbor_child_on_subface (const unsigned int face_no,
                                      const unsigned int subface_no) const;

  std::pair<unsigned int, unsigned int>
  periodic_neighbor_of_coarser_periodic_neighbor (const unsigned face_no) const;

  int
  periodic_neighbor_index (const unsigned int i) const;

  int
  periodic_neighbor_level (const unsigned int i) const;

  unsigned int
  periodic_neighbor_of_periodic_neighbor (const unsigned int i) const;

  unsigned int
  periodic_neighbor_face_no (const unsigned int i) const;

  bool
  periodic_neighbor_is_coarser (const unsigned int i) const;

  bool at_boundary (const unsigned int i) const;

  bool at_boundary () const;

  bool has_boundary_lines () const;
  RefinementCase<dim> refine_flag_set () const;

  void set_refine_flag (const RefinementCase<dim> ref_case = RefinementCase<dim>::isotropic_refinement) const;

  void clear_refine_flag () const;

  bool flag_for_face_refinement (const unsigned int face_no,
                                 const RefinementCase<dim-1> &face_refinement_case=RefinementCase<dim-1>::isotropic_refinement) const;

  bool flag_for_line_refinement (const unsigned int line_no) const;

  ::internal::SubfaceCase<dim> subface_case(const unsigned int face_no) const;

  bool coarsen_flag_set () const;

  void set_coarsen_flag () const;

  void clear_coarsen_flag () const;
  types::material_id material_id () const;

  void set_material_id (const types::material_id new_material_id) const;

  void recursively_set_material_id (const types::material_id new_material_id) const;
  types::subdomain_id subdomain_id () const;

  void set_subdomain_id (const types::subdomain_id new_subdomain_id) const;

  types::subdomain_id level_subdomain_id () const;

  void set_level_subdomain_id (const types::subdomain_id new_level_subdomain_id) const;


  void recursively_set_subdomain_id (const types::subdomain_id new_subdomain_id) const;
  bool direction_flag () const;

  unsigned int active_cell_index () const;

  int parent_index () const;

  TriaIterator<CellAccessor<dim,spacedim> >
  parent () const;

  bool active () const;

  bool is_locally_owned () const;

  bool is_locally_owned_on_level () const;

  bool is_ghost () const;

  bool is_artificial () const;

  bool point_inside (const Point<spacedim> &p) const;

  void set_neighbor (const unsigned int i,
                     const TriaIterator<CellAccessor<dim, spacedim> > &pointer) const;

  CellId id() const;

  DeclException0 (ExcRefineCellNotActive);
  DeclException0 (ExcCellFlaggedForRefinement);
  DeclException0 (ExcCellFlaggedForCoarsening);

protected:
  unsigned int neighbor_of_neighbor_internal (const unsigned int neighbor) const;

  template<int dim_,int spacedim_ >
  bool point_inside_codim(const Point<spacedim_> &p) const;



private:
  void set_active_cell_index (const unsigned int active_cell_index);

  void set_parent (const unsigned int parent_index);

  void set_direction_flag (const bool new_direction_flag) const;
  void operator = (const CellAccessor<dim, spacedim> &);

  template <int, int> friend class Triangulation;

  friend struct ::internal::Triangulation::Implementation;
};



/* -------------- declaration of explicit
   specializations and general templates ------------- */


template <int structdim, int dim, int spacedim>
template <typename OtherAccessor>
InvalidAccessor<structdim, dim, spacedim>::
InvalidAccessor (const OtherAccessor &)
{
  Assert (false,
          ExcMessage ("You are attempting an illegal conversion between "
                      "iterator/accessor types. The constructor you call "
                      "only exists to make certain template constructs "
                      "easier to write as dimension independent code but "
                      "the conversion is not valid in the current context."));
}



template <int structdim, int dim, int spacedim>
template <int structdim2, int dim2, int spacedim2>
TriaAccessor<structdim,dim,spacedim>::
TriaAccessor (const InvalidAccessor<structdim2,dim2,spacedim2> &)
{
  Assert (false,
          ExcMessage ("You are attempting an illegal conversion between "
                      "iterator/accessor types. The constructor you call "
                      "only exists to make certain template constructs "
                      "easier to write as dimension independent code but "
                      "the conversion is not valid in the current context."));
}



template <int dim, int spacedim>
template <int structdim2, int dim2, int spacedim2>
CellAccessor<dim,spacedim>::
CellAccessor (const InvalidAccessor<structdim2,dim2,spacedim2> &)
{
  Assert (false,
          ExcMessage ("You are attempting an illegal conversion between "
                      "iterator/accessor types. The constructor you call "
                      "only exists to make certain template constructs "
                      "easier to write as dimension independent code but "
                      "the conversion is not valid in the current context."));
}



template <int structdim, int dim, int spacedim>
template <int structdim2, int dim2, int spacedim2>
TriaAccessor<structdim,dim,spacedim>::
TriaAccessor (const TriaAccessor<structdim2,dim2,spacedim2> &)
{
  Assert (false,
          ExcMessage ("You are attempting an illegal conversion between "
                      "iterator/accessor types. The constructor you call "
                      "only exists to make certain template constructs "
                      "easier to write as dimension independent code but "
                      "the conversion is not valid in the current context."));
}



template <int dim, int spacedim>
template <int structdim2, int dim2, int spacedim2>
CellAccessor<dim,spacedim>::
CellAccessor (const TriaAccessor<structdim2,dim2,spacedim2> &)
{
  Assert (false,
          ExcMessage ("You are attempting an illegal conversion between "
                      "iterator/accessor types. The constructor you call "
                      "only exists to make certain template constructs "
                      "easier to write as dimension independent code but "
                      "the conversion is not valid in the current context."));
}


#ifndef DOXYGEN

template <> bool CellAccessor<1,1>::point_inside (const Point<1> &) const;
template <> bool CellAccessor<2,2>::point_inside (const Point<2> &) const;
template <> bool CellAccessor<3,3>::point_inside (const Point<3> &) const;
template <> bool CellAccessor<1,2>::point_inside (const Point<2> &) const;
template <> bool CellAccessor<1,3>::point_inside (const Point<3> &) const;
template <> bool CellAccessor<2,3>::point_inside (const Point<3> &) const;
// -------------------------------------------------------------------

template <> void TriaAccessor<3,3,3>::set_all_manifold_ids (const types::manifold_id) const;

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

// include more templates in debug and optimized mode
#  include "tria_accessor.templates.h"

#endif