dof_accessor.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2019 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_dof_accessor_h
#define dealii_dof_accessor_h


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

#include <deal.II/dofs/dof_handler.h>

#include <deal.II/grid/tria_accessor.h>

#include <deal.II/hp/dof_handler.h>

#include <vector>

DEAL_II_NAMESPACE_OPEN

template <typename number>
class FullMatrix;
template <typename number>
class Vector;
template <typename number>
class AffineConstraints;

template <typename Accessor>
class TriaRawIterator;

template <int, int>
class FiniteElement;


namespace internal
{
  namespace DoFCellAccessorImplementation
  {
    struct Implementation;
  }

  namespace DoFHandlerImplementation
  {
    struct Implementation;
    namespace Policy
    {
      struct Implementation;
    }
  } // namespace DoFHandlerImplementation

  namespace hp
  {
    namespace DoFHandlerImplementation
    {
      struct Implementation;
    }
  } // namespace hp
} // namespace internal

// note: the file dof_accessor.templates.h is included at the end of
// this file.  this includes a lot of templates and thus makes
// compilation slower, but at the same time allows for more aggressive
// inlining and thus faster code.


namespace internal
{
  namespace DoFAccessorImplementation
  {
    template <int structdim, int dim, int spacedim>
    struct Inheritance
    {
      using BaseClass = ::TriaAccessor<structdim, dim, spacedim>;
    };


    template <int dim, int spacedim>
    struct Inheritance<dim, dim, spacedim>
    {
      using BaseClass = ::CellAccessor<dim, spacedim>;
    };
  } // namespace DoFAccessorImplementation
} // namespace internal


/* -------------------------------------------------------------------------- */



template <int structdim, typename DoFHandlerType, bool level_dof_access>
class DoFAccessor
  : public ::internal::DoFAccessorImplementation::Inheritance<
      structdim,
      DoFHandlerType::dimension,
      DoFHandlerType::space_dimension>::BaseClass
{
public:
  static const unsigned int dimension = DoFHandlerType::dimension;

  static const unsigned int space_dimension = DoFHandlerType::space_dimension;

  using BaseClass = typename ::internal::DoFAccessorImplementation::
    Inheritance<structdim, dimension, space_dimension>::BaseClass;

  using AccessorData = DoFHandlerType;

  DoFAccessor();

  DoFAccessor(const Triangulation<DoFHandlerType::dimension,
                                  DoFHandlerType::space_dimension> *tria,
              const int                                             level,
              const int                                             index,
              const DoFHandlerType *dof_handler);

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

  template <int dim2, class DoFHandlerType2, bool level_dof_access2>
  DoFAccessor(const DoFAccessor<dim2, DoFHandlerType2, level_dof_access2> &);

  template <bool level_dof_access2>
  DoFAccessor(
    const DoFAccessor<structdim, DoFHandlerType, level_dof_access2> &);

  DoFAccessor<structdim, DoFHandlerType, level_dof_access> &
  operator=(const DoFAccessor<structdim, DoFHandlerType, level_dof_access>
              &da) = delete;

  const DoFHandlerType &
  get_dof_handler() const;

  template <bool level_dof_access2>
  void
  copy_from(const DoFAccessor<structdim, DoFHandlerType, level_dof_access2> &a);

  void
  copy_from(const TriaAccessorBase<structdim,
                                   DoFHandlerType::dimension,
                                   DoFHandlerType::space_dimension> &da);

  static bool
  is_level_cell();

  TriaIterator<DoFAccessor<structdim, DoFHandlerType, level_dof_access>>
  child(const unsigned int c) const;

  typename ::internal::DoFHandlerImplementation::
    Iterators<DoFHandlerType, level_dof_access>::line_iterator
    line(const unsigned int i) const;

  typename ::internal::DoFHandlerImplementation::
    Iterators<DoFHandlerType, level_dof_access>::quad_iterator
    quad(const unsigned int i) const;

  void
  get_dof_indices(
    std::vector<types::global_dof_index> &dof_indices,
    const unsigned int fe_index = DoFHandlerType::default_fe_index) const;

  void
  get_mg_dof_indices(
    const int                             level,
    std::vector<types::global_dof_index> &dof_indices,
    const unsigned int fe_index = DoFHandlerType::default_fe_index) const;

  void
  set_mg_dof_indices(
    const int                                   level,
    const std::vector<types::global_dof_index> &dof_indices,
    const unsigned int fe_index = DoFHandlerType::default_fe_index);

  types::global_dof_index
  vertex_dof_index(
    const unsigned int vertex,
    const unsigned int i,
    const unsigned int fe_index = DoFHandlerType::default_fe_index) const;

  types::global_dof_index
  mg_vertex_dof_index(
    const int          level,
    const unsigned int vertex,
    const unsigned int i,
    const unsigned int fe_index = DoFHandlerType::default_fe_index) const;

  types::global_dof_index
  dof_index(
    const unsigned int i,
    const unsigned int fe_index = DoFHandlerType::default_fe_index) const;

  types::global_dof_index
  mg_dof_index(const int level, const unsigned int i) const;

  unsigned int
  n_active_fe_indices() const;

  unsigned int
  nth_active_fe_index(const unsigned int n) const;

  std::set<unsigned int>
  get_active_fe_indices() const;

  bool
  fe_index_is_active(const unsigned int fe_index) const;

  const FiniteElement<DoFHandlerType::dimension,
                      DoFHandlerType::space_dimension> &
  get_fe(const unsigned int fe_index) const;

  DeclExceptionMsg(ExcInvalidObject,
                   "This accessor object has not been "
                   "associated with any DoFHandler object.");
  DeclException0(ExcVectorNotEmpty);
  DeclException0(ExcVectorDoesNotMatch);
  DeclException0(ExcMatrixDoesNotMatch);
  DeclException0(ExcNotActive);
  DeclException0(ExcCantCompareIterators);

protected:
  DoFHandlerType *dof_handler;

public:
  template <int dim2, class DoFHandlerType2, bool level_dof_access2>
  bool
  operator==(
    const DoFAccessor<dim2, DoFHandlerType2, level_dof_access2> &) const;

  template <int dim2, class DoFHandlerType2, bool level_dof_access2>
  bool
  operator!=(
    const DoFAccessor<dim2, DoFHandlerType2, level_dof_access2> &) const;

protected:
  void
  set_dof_handler(DoFHandlerType *dh);

  void
  set_dof_index(
    const unsigned int            i,
    const types::global_dof_index index,
    const unsigned int fe_index = DoFHandlerType::default_fe_index) const;

  void
  set_mg_dof_index(const int                     level,
                   const unsigned int            i,
                   const types::global_dof_index index) const;

  void
  set_vertex_dof_index(
    const unsigned int            vertex,
    const unsigned int            i,
    const types::global_dof_index index,
    const unsigned int fe_index = DoFHandlerType::default_fe_index) const;

  void
  set_mg_vertex_dof_index(
    const int                     level,
    const unsigned int            vertex,
    const unsigned int            i,
    const types::global_dof_index index,
    const unsigned int fe_index = DoFHandlerType::default_fe_index) const;

  template <typename>
  friend class TriaRawIterator;
  template <int, class, bool>
  friend class DoFAccessor;

private:
  template <int dim, int spacedim>
  friend class DoFHandler;
  template <int dim, int spacedim>
  friend class hp::DoFHandler;

  friend struct ::internal::DoFHandlerImplementation::Policy::
    Implementation;
  friend struct ::internal::DoFHandlerImplementation::Implementation;
  friend struct ::internal::hp::DoFHandlerImplementation::Implementation;
  friend struct ::internal::DoFCellAccessorImplementation::Implementation;
  friend struct ::internal::DoFAccessorImplementation::Implementation;
};



template <template <int, int> class DoFHandlerType,
          int  spacedim,
          bool level_dof_access>
class DoFAccessor<0, DoFHandlerType<1, spacedim>, level_dof_access>
  : public TriaAccessor<0, 1, spacedim>
{
public:
  static const unsigned int dimension = 1;

  static const unsigned int space_dimension = spacedim;

  using BaseClass = TriaAccessor<0, 1, spacedim>;

  using AccessorData = DoFHandlerType<1, spacedim>;

  DoFAccessor();

  DoFAccessor(
    const Triangulation<1, spacedim> *                      tria,
    const typename TriaAccessor<0, 1, spacedim>::VertexKind vertex_kind,
    const unsigned int                                      vertex_index,
    const DoFHandlerType<1, spacedim> *                     dof_handler);

  DoFAccessor(const Triangulation<1, spacedim> *,
              const int                                      = 0,
              const int                                      = 0,
              const DoFHandlerType<1, spacedim> *dof_handler = 0);

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

  template <int dim2, class DoFHandlerType2, bool level_dof_access2>
  DoFAccessor(const DoFAccessor<dim2, DoFHandlerType2, level_dof_access2> &);

  DoFAccessor<0, DoFHandlerType<1, spacedim>, level_dof_access> &
  operator=(const DoFAccessor<0, DoFHandlerType<1, spacedim>, level_dof_access>
              &da) = delete;

  const DoFHandlerType<1, spacedim> &
  get_dof_handler() const;

  template <bool level_dof_access2>
  void
  copy_from(
    const DoFAccessor<0, DoFHandlerType<1, spacedim>, level_dof_access2> &a);

  void
  copy_from(const TriaAccessorBase<0, 1, spacedim> &da);

  TriaIterator<DoFAccessor<0, DoFHandlerType<1, spacedim>, level_dof_access>>
  child(const unsigned int c) const;

  typename ::internal::DoFHandlerImplementation::
    Iterators<DoFHandlerType<1, spacedim>, level_dof_access>::line_iterator
    line(const unsigned int i) const;

  typename ::internal::DoFHandlerImplementation::
    Iterators<DoFHandlerType<1, spacedim>, level_dof_access>::quad_iterator
    quad(const unsigned int i) const;

  void
  get_dof_indices(
    std::vector<types::global_dof_index> &dof_indices,
    const unsigned int fe_index = AccessorData::default_fe_index) const;

  void
  get_mg_dof_indices(
    const int                             level,
    std::vector<types::global_dof_index> &dof_indices,
    const unsigned int fe_index = AccessorData::default_fe_index) const;

  types::global_dof_index
  vertex_dof_index(
    const unsigned int vertex,
    const unsigned int i,
    const unsigned int fe_index = AccessorData::default_fe_index) const;

  types::global_dof_index
  dof_index(const unsigned int i,
            const unsigned int fe_index = AccessorData::default_fe_index) const;

  unsigned int
  n_active_fe_indices() const;

  unsigned int
  nth_active_fe_index(const unsigned int n) const;

  bool
  fe_index_is_active(const unsigned int fe_index) const;

  const FiniteElement<DoFHandlerType<1, spacedim>::dimension,
                      DoFHandlerType<1, spacedim>::space_dimension> &
  get_fe(const unsigned int fe_index) const;

  DeclExceptionMsg(ExcInvalidObject,
                   "This accessor object has not been "
                   "associated with any DoFHandler object.");
  DeclException0(ExcVectorNotEmpty);
  DeclException0(ExcVectorDoesNotMatch);
  DeclException0(ExcMatrixDoesNotMatch);
  DeclException0(ExcNotActive);
  DeclException0(ExcCantCompareIterators);

protected:
  DoFHandlerType<1, spacedim> *dof_handler;

  template <int dim2, class DoFHandlerType2, bool level_dof_access2>
  bool
  operator==(
    const DoFAccessor<dim2, DoFHandlerType2, level_dof_access2> &) const;

  template <int dim2, class DoFHandlerType2, bool level_dof_access2>
  bool
  operator!=(
    const DoFAccessor<dim2, DoFHandlerType2, level_dof_access2> &) const;

  void set_dof_handler(DoFHandlerType<1, spacedim> *dh);

  void
  set_dof_index(
    const unsigned int            i,
    const types::global_dof_index index,
    const unsigned int fe_index = AccessorData::default_fe_index) const;

  void
  set_vertex_dof_index(
    const unsigned int            vertex,
    const unsigned int            i,
    const types::global_dof_index index,
    const unsigned int fe_index = AccessorData::default_fe_index) const;

  template <typename>
  friend class TriaRawIterator;


  template <int, int>
  friend class DoFHandler;
  template <int, int>
  friend class hp::DoFHandler;

  friend struct ::internal::DoFHandlerImplementation::Policy::
    Implementation;
  friend struct ::internal::DoFHandlerImplementation::Implementation;
  friend struct ::internal::hp::DoFHandlerImplementation::Implementation;
  friend struct ::internal::DoFCellAccessorImplementation::Implementation;
};



/* -------------------------------------------------------------------------- */


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

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

  DoFInvalidAccessor(const DoFInvalidAccessor<structdim, dim, spacedim> &);

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

  void
  set_dof_index(const unsigned int            i,
                const types::global_dof_index index,
                const unsigned int            fe_index =
                  DoFHandler<dim, spacedim>::default_fe_index) const;
};



/* -------------------------------------------------------------------------- */


template <typename DoFHandlerType, bool level_dof_access>
class DoFCellAccessor : public DoFAccessor<DoFHandlerType::dimension,
                                           DoFHandlerType,
                                           level_dof_access>
{
public:
  static const unsigned int dim = DoFHandlerType::dimension;

  static const unsigned int spacedim = DoFHandlerType::space_dimension;


  using AccessorData = DoFHandlerType;

  using BaseClass =
    DoFAccessor<DoFHandlerType::dimension, DoFHandlerType, level_dof_access>;

  using Container = DoFHandlerType;

  using face_iterator = TriaIterator<DoFAccessor<DoFHandlerType::dimension - 1,
                                                 DoFHandlerType,
                                                 level_dof_access>>;

  DoFCellAccessor(const Triangulation<DoFHandlerType::dimension,
                                      DoFHandlerType::space_dimension> *tria,
                  const int                                             level,
                  const int                                             index,
                  const AccessorData *local_data);

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

  template <int dim2, class DoFHandlerType2, bool level_dof_access2>
  explicit DoFCellAccessor(
    const DoFAccessor<dim2, DoFHandlerType2, level_dof_access2> &);

  DoFCellAccessor<DoFHandlerType, level_dof_access> &
  operator=(const DoFCellAccessor<DoFHandlerType, level_dof_access> &da) =
    delete;

  TriaIterator<DoFCellAccessor<DoFHandlerType, level_dof_access>>
  parent() const;

  TriaIterator<DoFCellAccessor<DoFHandlerType, level_dof_access>>
  neighbor(const unsigned int i) const;

  TriaIterator<DoFCellAccessor<DoFHandlerType, level_dof_access>>
  periodic_neighbor(const unsigned int i) const;

  TriaIterator<DoFCellAccessor<DoFHandlerType, level_dof_access>>
  neighbor_or_periodic_neighbor(const unsigned int i) const;

  TriaIterator<DoFCellAccessor<DoFHandlerType, level_dof_access>>
  child(const unsigned int i) const;

  face_iterator
  face(const unsigned int i) const;

  TriaIterator<DoFCellAccessor<DoFHandlerType, level_dof_access>>
  neighbor_child_on_subface(const unsigned int face_no,
                            const unsigned int subface_no) const;

  TriaIterator<DoFCellAccessor<DoFHandlerType, level_dof_access>>
  periodic_neighbor_child_on_subface(const unsigned int face_no,
                                     const unsigned int subface_no) const;

  template <class InputVector, typename number>
  void
  get_dof_values(const InputVector &values, Vector<number> &local_values) const;

  template <class InputVector, typename ForwardIterator>
  void
  get_dof_values(const InputVector &values,
                 ForwardIterator    local_values_begin,
                 ForwardIterator    local_values_end) const;

  template <class InputVector, typename ForwardIterator>
  void
  get_dof_values(
    const AffineConstraints<typename InputVector::value_type> &constraints,
    const InputVector &                                        values,
    ForwardIterator local_values_begin,
    ForwardIterator local_values_end) const;

  template <class OutputVector, typename number>
  void
  set_dof_values(const Vector<number> &local_values,
                 OutputVector &        values) const;

  template <class InputVector, typename number>
  void
  get_interpolated_dof_values(
    const InputVector &values,
    Vector<number> &   interpolated_values,
    const unsigned int fe_index = DoFHandlerType::default_fe_index) const;

  template <class OutputVector, typename number>
  void
  set_dof_values_by_interpolation(
    const Vector<number> &local_values,
    OutputVector &        values,
    const unsigned int    fe_index = DoFHandlerType::default_fe_index) const;

  template <typename number, typename OutputVector>
  void
  distribute_local_to_global(const Vector<number> &local_source,
                             OutputVector &        global_destination) const;

  template <typename ForwardIterator, typename OutputVector>
  void
  distribute_local_to_global(ForwardIterator local_source_begin,
                             ForwardIterator local_source_end,
                             OutputVector &  global_destination) const;

  template <typename ForwardIterator, typename OutputVector>
  void
  distribute_local_to_global(
    const AffineConstraints<typename OutputVector::value_type> &constraints,
    ForwardIterator local_source_begin,
    ForwardIterator local_source_end,
    OutputVector &  global_destination) const;

  template <typename number, typename OutputMatrix>
  void
  distribute_local_to_global(const FullMatrix<number> &local_source,
                             OutputMatrix &global_destination) const;

  template <typename number, typename OutputMatrix, typename OutputVector>
  void
  distribute_local_to_global(const FullMatrix<number> &local_matrix,
                             const Vector<number> &    local_vector,
                             OutputMatrix &            global_matrix,
                             OutputVector &            global_vector) const;

  void
  get_active_or_mg_dof_indices(
    std::vector<types::global_dof_index> &dof_indices) const;

  void
  get_dof_indices(std::vector<types::global_dof_index> &dof_indices) const;

  void
  get_mg_dof_indices(std::vector<types::global_dof_index> &dof_indices) const;

  const FiniteElement<DoFHandlerType::dimension,
                      DoFHandlerType::space_dimension> &
  get_fe() const;

  unsigned int
  active_fe_index() const;

  void
  set_active_fe_index(const unsigned int i) const;
  void
  set_dof_indices(const std::vector<types::global_dof_index> &dof_indices);

  void
  set_mg_dof_indices(const std::vector<types::global_dof_index> &dof_indices);

  void
  update_cell_dof_indices_cache() const;

  unsigned int
  future_fe_index() const;

  void
  set_future_fe_index(const unsigned int i) const;

  bool
  future_fe_index_set() const;

  void
  clear_future_fe_index() const;
private:
  template <int dim, int spacedim>
  friend class DoFHandler;
  friend struct ::internal::DoFCellAccessorImplementation::Implementation;
};


template <int sd, typename DoFHandlerType, bool level_dof_access>
inline bool
DoFAccessor<sd, DoFHandlerType, level_dof_access>::is_level_cell()
{
  return level_dof_access;
}



template <int structdim, int dim, int spacedim>
template <typename OtherAccessor>
DoFInvalidAccessor<structdim, dim, spacedim>::DoFInvalidAccessor(
  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."));
}



DEAL_II_NAMESPACE_CLOSE

// include more templates
#include "dof_accessor.templates.h"


#endif