mg_transfer.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2001 - 2018 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_mg_transfer_h
#define dealii_mg_transfer_h

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

#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/la_parallel_vector.h>

#include <deal.II/lac/vector_memory.h>

#include <deal.II/multigrid/mg_base.h>
#include <deal.II/multigrid/mg_constrained_dofs.h>
#include <deal.II/base/mg_level_object.h>

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

#include <memory>


DEAL_II_NAMESPACE_OPEN


namespace internal
{
  template <typename VectorType>
  struct MatrixSelector
  {
    typedef ::SparsityPattern Sparsity;
    typedef ::SparseMatrix<typename VectorType::value_type> Matrix;

    template <typename SparsityPatternType, typename DoFHandlerType>
    static void reinit(Matrix &matrix, Sparsity &sparsity, int level, const SparsityPatternType &sp, const DoFHandlerType &)
    {
      sparsity.copy_from (sp);
      (void)level;
      matrix.reinit (sparsity);
    }
  };

#ifdef DEAL_II_WITH_TRILINOS
  template <typename Number>
  struct MatrixSelector<LinearAlgebra::distributed::Vector<Number> >
  {
    typedef ::TrilinosWrappers::SparsityPattern Sparsity;
    typedef ::TrilinosWrappers::SparseMatrix Matrix;

    template <typename SparsityPatternType, typename DoFHandlerType>
    static void reinit(Matrix &matrix, Sparsity &, int level, const SparsityPatternType &sp, DoFHandlerType &dh)
    {
      const parallel::Triangulation<DoFHandlerType::dimension,DoFHandlerType::space_dimension> *dist_tria =
        dynamic_cast<const parallel::Triangulation<DoFHandlerType::dimension,DoFHandlerType::space_dimension>*>
        (&(dh.get_triangulation()));
      MPI_Comm communicator = dist_tria != nullptr ?
                              dist_tria->get_communicator() :
                              MPI_COMM_SELF;

      matrix.reinit(dh.locally_owned_mg_dofs(level+1),
                    dh.locally_owned_mg_dofs(level),
                    sp, communicator, true);
    }

  };

  template <>
  struct MatrixSelector<::TrilinosWrappers::MPI::Vector>
  {
    typedef ::TrilinosWrappers::SparsityPattern Sparsity;
    typedef ::TrilinosWrappers::SparseMatrix Matrix;

    template <typename SparsityPatternType, typename DoFHandlerType>
    static void reinit(Matrix &matrix, Sparsity &, int level, const SparsityPatternType &sp, DoFHandlerType &dh)
    {
      const parallel::Triangulation<DoFHandlerType::dimension,DoFHandlerType::space_dimension> *dist_tria =
        dynamic_cast<const parallel::Triangulation<DoFHandlerType::dimension,DoFHandlerType::space_dimension>*>
        (&(dh.get_triangulation()));
      MPI_Comm communicator = dist_tria != nullptr ?
                              dist_tria->get_communicator() :
                              MPI_COMM_SELF;
      matrix.reinit(dh.locally_owned_mg_dofs(level+1),
                    dh.locally_owned_mg_dofs(level),
                    sp, communicator, true);
    }

  };

#ifdef DEAL_II_WITH_MPI
  template <>
  struct MatrixSelector<::LinearAlgebra::EpetraWrappers::Vector>
  {
    typedef ::TrilinosWrappers::SparsityPattern Sparsity;
    typedef ::TrilinosWrappers::SparseMatrix Matrix;

    template <typename SparsityPatternType, typename DoFHandlerType>
    static void reinit(Matrix &matrix, Sparsity &, int level, const SparsityPatternType &sp, DoFHandlerType &dh)
    {
      const parallel::Triangulation<DoFHandlerType::dimension,DoFHandlerType::space_dimension> *dist_tria =
        dynamic_cast<const parallel::Triangulation<DoFHandlerType::dimension,DoFHandlerType::space_dimension>*>
        (&(dh.get_triangulation()));
      MPI_Comm communicator = dist_tria != nullptr ?
                              dist_tria->get_communicator() :
                              MPI_COMM_SELF;
      matrix.reinit(dh.locally_owned_mg_dofs(level+1),
                    dh.locally_owned_mg_dofs(level),
                    sp, communicator, true);
    }
  };
#endif

#else
  // ! DEAL_II_WITH_TRILINOS
  template <typename Number>
  struct MatrixSelector<LinearAlgebra::distributed::Vector<Number> >
  {
    typedef ::SparsityPattern Sparsity;
    typedef ::SparseMatrix<Number> Matrix;

    template <typename SparsityPatternType, typename DoFHandlerType>
    static void reinit(Matrix &, Sparsity &, int, const SparsityPatternType &, const DoFHandlerType &)
    {
      AssertThrow(false, ExcNotImplemented(
                    "ERROR: MGTransferPrebuilt with LinearAlgebra::distributed::Vector currently "
                    "needs deal.II to be configured with Trilinos."));
    }
  };

#endif
}

/*
 * MGTransferBase is defined in mg_base.h
 */




template <typename VectorType>
class MGLevelGlobalTransfer : public MGTransferBase<VectorType>
{
public:

  void clear ();

  template <int dim, class InVector, int spacedim>
  void
  copy_to_mg (const DoFHandler<dim,spacedim> &mg_dof,
              MGLevelObject<VectorType>      &dst,
              const InVector                 &src) const;

  template <int dim, class OutVector, int spacedim>
  void
  copy_from_mg (const DoFHandler<dim,spacedim>  &mg_dof,
                OutVector                       &dst,
                const MGLevelObject<VectorType> &src) const;

  template <int dim, class OutVector, int spacedim>
  void
  copy_from_mg_add (const DoFHandler<dim,spacedim>  &mg_dof,
                    OutVector                       &dst,
                    const MGLevelObject<VectorType> &src) const;

  void
  set_component_to_block_map (const std::vector<unsigned int> &map);

  std::size_t memory_consumption () const;

  void print_indices(std::ostream &os) const;

protected:

  template <int dim, int spacedim>
  void fill_and_communicate_copy_indices(const DoFHandler<dim,spacedim> &mg_dof);

  std::vector<types::global_dof_index> sizes;

  std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > >
  copy_indices;

  std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > >
  copy_indices_global_mine;

  std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > >
  copy_indices_level_mine;

  bool perform_plain_copy;

  std::vector<unsigned int> component_to_block_map;

  SmartPointer<const MGConstrainedDoFs, MGLevelGlobalTransfer<VectorType> > mg_constrained_dofs;
};



template <typename Number>
class MGLevelGlobalTransfer<LinearAlgebra::distributed::Vector<Number> > : public MGTransferBase<LinearAlgebra::distributed::Vector<Number> >
{
public:

  void clear ();

  template <int dim, typename Number2, int spacedim>
  void
  copy_to_mg (const DoFHandler<dim,spacedim>                        &mg_dof,
              MGLevelObject<LinearAlgebra::distributed::Vector<Number> > &dst,
              const LinearAlgebra::distributed::Vector<Number2>          &src) const;

  template <int dim, typename Number2, int spacedim>
  void
  copy_from_mg (const DoFHandler<dim,spacedim>                              &mg_dof,
                LinearAlgebra::distributed::Vector<Number2>                      &dst,
                const MGLevelObject<LinearAlgebra::distributed::Vector<Number> > &src) const;

  template <int dim, typename Number2, int spacedim>
  void
  copy_from_mg_add (const DoFHandler<dim,spacedim>                              &mg_dof,
                    LinearAlgebra::distributed::Vector<Number2>                      &dst,
                    const MGLevelObject<LinearAlgebra::distributed::Vector<Number> > &src) const;

  void
  set_component_to_block_map (const std::vector<unsigned int> &map);

  std::size_t memory_consumption () const;

  void print_indices(std::ostream &os) const;

protected:

  template <int dim, typename Number2, int spacedim>
  void
  copy_to_mg (const DoFHandler<dim,spacedim>                        &mg_dof,
              MGLevelObject<LinearAlgebra::distributed::Vector<Number> > &dst,
              const LinearAlgebra::distributed::Vector<Number2>          &src,
              const bool solution_transfer) const;

  template <int dim, int spacedim>
  void fill_and_communicate_copy_indices(const DoFHandler<dim,spacedim> &mg_dof);

  std::vector<types::global_dof_index> sizes;

  std::vector<std::vector<std::pair<unsigned int, unsigned int> > >
  copy_indices;


  std::vector<std::vector<std::pair<unsigned int, unsigned int> > >
  solution_copy_indices;

  std::vector<std::vector<std::pair<unsigned int, unsigned int> > >
  copy_indices_global_mine;

  std::vector<std::vector<std::pair<unsigned int, unsigned int> > >
  solution_copy_indices_global_mine;

  std::vector<std::vector<std::pair<unsigned int, unsigned int> > >
  copy_indices_level_mine;

  std::vector<std::vector<std::pair<unsigned int, unsigned int> > >
  solution_copy_indices_level_mine;

  bool perform_plain_copy;

  bool perform_renumbered_plain_copy;

  std::vector<unsigned int> component_to_block_map;

  SmartPointer<const MGConstrainedDoFs, MGLevelGlobalTransfer<LinearAlgebra::distributed::Vector<Number> > > mg_constrained_dofs;

  mutable LinearAlgebra::distributed::Vector<Number> ghosted_global_vector;

  mutable LinearAlgebra::distributed::Vector<Number> solution_ghosted_global_vector;

  mutable MGLevelObject<LinearAlgebra::distributed::Vector<Number> > ghosted_level_vector;

  mutable MGLevelObject<LinearAlgebra::distributed::Vector<Number> > solution_ghosted_level_vector;

};



template <typename VectorType>
class MGTransferPrebuilt : public MGLevelGlobalTransfer<VectorType>
{
public:
  MGTransferPrebuilt () = default;

  MGTransferPrebuilt (const MGConstrainedDoFs &mg_constrained_dofs);

  DEAL_II_DEPRECATED
  MGTransferPrebuilt (const ConstraintMatrix &constraints,
                      const MGConstrainedDoFs &mg_constrained_dofs);

  virtual ~MGTransferPrebuilt () = default;

  void initialize_constraints (const MGConstrainedDoFs &mg_constrained_dofs);

  DEAL_II_DEPRECATED
  void initialize_constraints (const ConstraintMatrix &constraints,
                               const MGConstrainedDoFs &mg_constrained_dofs);

  void clear ();

  template <int dim, int spacedim>
  void build_matrices (const DoFHandler<dim,spacedim> &mg_dof);

  virtual void prolongate (const unsigned int to_level,
                           VectorType         &dst,
                           const VectorType   &src) const;

  virtual void restrict_and_add (const unsigned int from_level,
                                 VectorType         &dst,
                                 const VectorType   &src) const;

  DeclException0(ExcNoProlongation);

  DeclException0(ExcMatricesNotBuilt);

  std::size_t memory_consumption () const;

  void print_matrices(std::ostream &os) const;

private:

  std::vector<std::shared_ptr<typename internal::MatrixSelector<VectorType>::Sparsity> > prolongation_sparsities;

  std::vector<std::shared_ptr<typename internal::MatrixSelector<VectorType>::Matrix> > prolongation_matrices;

  std::vector<std::vector<bool> > interface_dofs;
};


DEAL_II_NAMESPACE_CLOSE

#endif