mg_transfer_matrix_free.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 2017 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_matrix_free_h
#define dealii_mg_transfer_matrix_free_h

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

#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/la_parallel_block_vector.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/multigrid/mg_transfer.h>
#include <deal.II/multigrid/mg_transfer_internal.h>
#include <deal.II/base/vectorization.h>

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


DEAL_II_NAMESPACE_OPEN



template <int dim, typename Number>
class MGTransferMatrixFree : public MGLevelGlobalTransfer<LinearAlgebra::distributed::Vector<Number> >
{
public:
  MGTransferMatrixFree ();

  MGTransferMatrixFree (const MGConstrainedDoFs &mg_constrained_dofs);

  virtual ~MGTransferMatrixFree () = default;

  void initialize_constraints (const MGConstrainedDoFs &mg_constrained_dofs);

  void clear ();

  void build (const DoFHandler<dim,dim> &mg_dof);

  virtual void prolongate (const unsigned int                           to_level,
                           LinearAlgebra::distributed::Vector<Number>       &dst,
                           const LinearAlgebra::distributed::Vector<Number> &src) const;

  virtual void restrict_and_add (const unsigned int from_level,
                                 LinearAlgebra::distributed::Vector<Number>       &dst,
                                 const LinearAlgebra::distributed::Vector<Number> &src) const;

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

  DeclException0(ExcNoProlongation);

  std::size_t memory_consumption () const;

private:

  unsigned int fe_degree;

  bool element_is_continuous;

  unsigned int n_components;

  unsigned int n_child_cell_dofs;

  std::vector<std::vector<unsigned int> > level_dof_indices;

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

  std::vector<unsigned int> n_owned_level_cells;

  AlignedVector<VectorizedArray<Number> > prolongation_matrix_1d;

  mutable AlignedVector<VectorizedArray<Number> > evaluation_data;

  std::vector<AlignedVector<VectorizedArray<Number> > > weights_on_refined;

  std::vector<std::vector<std::vector<unsigned short> > > dirichlet_indices;

  template <int degree>
  void do_prolongate_add(const unsigned int                           to_level,
                         LinearAlgebra::distributed::Vector<Number>       &dst,
                         const LinearAlgebra::distributed::Vector<Number> &src) const;

  template <int degree>
  void do_restrict_add(const unsigned int                           from_level,
                       LinearAlgebra::distributed::Vector<Number>       &dst,
                       const LinearAlgebra::distributed::Vector<Number> &src) const;
};


template <int dim, typename Number>
class MGTransferBlockMatrixFree : public MGTransferBase<LinearAlgebra::distributed::BlockVector<Number>>
{
public:
  MGTransferBlockMatrixFree () = default;

  MGTransferBlockMatrixFree (const MGConstrainedDoFs &mg_constrained_dofs);

  MGTransferBlockMatrixFree (const std::vector<MGConstrainedDoFs> &mg_constrained_dofs);

  virtual ~MGTransferBlockMatrixFree () = default;

  void initialize_constraints (const MGConstrainedDoFs &mg_constrained_dofs);

  void initialize_constraints (const std::vector<MGConstrainedDoFs> &mg_constrained_dofs);

  void clear ();

  void build (const DoFHandler<dim,dim> &mg_dof);

  void build (const std::vector<const DoFHandler<dim,dim>*> &mg_dof);

  virtual void prolongate (const unsigned int                                    to_level,
                           LinearAlgebra::distributed::BlockVector<Number>       &dst,
                           const LinearAlgebra::distributed::BlockVector<Number> &src) const;

  virtual void restrict_and_add (const unsigned int from_level,
                                 LinearAlgebra::distributed::BlockVector<Number>       &dst,
                                 const LinearAlgebra::distributed::BlockVector<Number> &src) const;

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

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

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

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

  std::size_t memory_consumption () const;

  static const bool supports_dof_handler_vector = true;

private:

  std::vector<MGTransferMatrixFree<dim,Number>> matrix_free_transfer_vector;

  const bool same_for_all;
};


//------------------------ templated functions -------------------------
#ifndef DOXYGEN


template <int dim, typename Number>
template <typename Number2, int spacedim>
void
MGTransferMatrixFree<dim,Number>::
interpolate_to_mg(const DoFHandler<dim,spacedim> &dof_handler,
                  MGLevelObject<LinearAlgebra::distributed::Vector<Number>> &dst,
                  const LinearAlgebra::distributed::Vector<Number2> &src) const
{
  const unsigned int min_level = dst.min_level();
  const unsigned int max_level = dst.max_level();

  Assert (max_level == dof_handler.get_triangulation().n_global_levels()-1,
          ExcDimensionMismatch(max_level,dof_handler.get_triangulation().n_global_levels()-1));

  const parallel::Triangulation<dim,spacedim> *p_tria =
    (dynamic_cast<const parallel::Triangulation<dim,spacedim>*>
     (&dof_handler.get_triangulation()));
  MPI_Comm mpi_communicator = p_tria != nullptr ? p_tria->get_communicator() : MPI_COMM_SELF;

  // resize the dst vector if it's empty or has incorrect size
  MGLevelObject<IndexSet> relevant_dofs(min_level,max_level);
  for (unsigned int level = min_level; level <= max_level; ++level)
    {
      DoFTools::extract_locally_relevant_level_dofs(dof_handler, level,
                                                    relevant_dofs[level]);
      if (dst[level].size() != dof_handler.locally_owned_mg_dofs(level).size() ||
          dst[level].local_size() != dof_handler.locally_owned_mg_dofs(level).n_elements())
        dst[level].reinit(dof_handler.locally_owned_mg_dofs(level),
                          relevant_dofs[level],
                          mpi_communicator);
    }

  const FiniteElement<dim,spacedim> &fe = dof_handler.get_fe();

  // copy fine level vector to active cells in MG hierarchy
  this->copy_to_mg (dof_handler, dst, src, true);

  // FIXME: maybe need to store hanging nodes constraints per level?
  // MGConstrainedDoFs does NOT keep this info right now, only periodicity constraints...
  dst[max_level].update_ghost_values();
  // do the transfer from level to level-1:
  for (unsigned int level=max_level; level > min_level; --level)
    {
      // auxiliary vector which always has ghost elements
      LinearAlgebra::distributed::Vector<Number>
      ghosted_vector(dof_handler.locally_owned_mg_dofs(level),
                     relevant_dofs[level],
                     mpi_communicator);
      ghosted_vector = dst[level];
      ghosted_vector.update_ghost_values();

      std::vector<Number> dof_values_coarse(fe.dofs_per_cell);
      Vector<Number> dof_values_fine(fe.dofs_per_cell);
      Vector<Number> tmp(fe.dofs_per_cell);
      std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
      typename DoFHandler<dim>::cell_iterator cell=dof_handler.begin(level-1);
      typename DoFHandler<dim>::cell_iterator endc=dof_handler.end(level-1);
      for ( ; cell != endc; ++cell)
        if (cell->is_locally_owned_on_level())
          {
            // if we get to a cell without children (== active), we can
            // skip it as there values should be already set by the
            // equivalent of copy_to_mg()
            if (!cell->has_children())
              continue;

            std::fill(dof_values_coarse.begin(), dof_values_coarse.end(), 0.);
            for (unsigned int child=0; child<cell->n_children(); ++child)
              {
                cell->child(child)->get_mg_dof_indices(dof_indices);
                for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
                  dof_values_fine(i) = ghosted_vector(dof_indices[i]);
                fe.get_restriction_matrix(child, cell->refinement_case()).vmult (tmp, dof_values_fine);
                for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
                  if (fe.restriction_is_additive(i))
                    dof_values_coarse[i] += tmp[i];
                  else if (tmp(i) != 0.)
                    dof_values_coarse[i] = tmp[i];
              }
            cell->get_mg_dof_indices(dof_indices);
            for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
              dst[level-1](dof_indices[i]) = dof_values_coarse[i];
          }

      dst[level-1].compress(VectorOperation::insert);
      dst[level-1].update_ghost_values();
    }
}



template <int dim, typename Number>
template <typename Number2, int spacedim>
void
MGTransferBlockMatrixFree<dim,Number>::
copy_to_mg (const DoFHandler<dim,spacedim> &mg_dof,
            MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &dst,
            const LinearAlgebra::distributed::BlockVector<Number2>         &src) const
{
  AssertDimension(matrix_free_transfer_vector.size(),1);
  Assert(same_for_all,
         ExcMessage("This object was initialized with support for usage with one "
                    "DoFHandler for each block, but this method assumes that "
                    "the same DoFHandler is used for all the blocks!"));
  const std::vector<const DoFHandler<dim, spacedim>*> mg_dofs (src.n_blocks(),&mg_dof);

  copy_to_mg(mg_dofs, dst, src);
}

template <int dim, typename Number>
template <typename Number2, int spacedim>
void
MGTransferBlockMatrixFree<dim,Number>::
copy_to_mg (const std::vector<const DoFHandler<dim,spacedim>*> &mg_dof,
            MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &dst,
            const LinearAlgebra::distributed::BlockVector<Number2>         &src) const
{
  const unsigned int n_blocks  = src.n_blocks();
  AssertDimension(mg_dof.size(), n_blocks);

  if (n_blocks==0)
    return;

  const unsigned int min_level = dst.min_level();
  const unsigned int max_level = dst.max_level();

  // this function is normally called within the Multigrid class with
  // dst == defect level block vector. At first run this vector is not
  // initialized. Do this below:
  {
    const parallel::Triangulation<dim,spacedim> *tria =
      (dynamic_cast<const parallel::Triangulation<dim,spacedim>*>
       (&(mg_dof[0]->get_triangulation())));
    for (unsigned int i=1; i<n_blocks; ++i)
      AssertThrow((dynamic_cast<const parallel::Triangulation<dim,spacedim>*>
                   (&(mg_dof[0]->get_triangulation())) == tria),
                  ExcMessage("The DoFHandler use different Triangulations!"));

    for (unsigned int level = min_level; level <= max_level; ++level)
      {
        dst[level].reinit(n_blocks);
        bool collect_size = false;
        for (unsigned int b = 0; b < n_blocks; ++b)
          {
            LinearAlgebra::distributed::Vector<Number> &v = dst[level].block(b);
            if (v.size() != mg_dof[b]->locally_owned_mg_dofs(level).size() ||
                v.local_size() != mg_dof[b]->locally_owned_mg_dofs(level).n_elements())
              {
                v.reinit(mg_dof[b]->locally_owned_mg_dofs(level), tria != nullptr ?
                         tria->get_communicator() : MPI_COMM_SELF);
                collect_size = true;
              }
            else
              v = 0.;
          }
        if (collect_size)
          dst[level].collect_sizes ();
      }
  }

  // FIXME: this a quite ugly as we need a temporary object:
  MGLevelObject<LinearAlgebra::distributed::Vector<Number>> dst_non_block(min_level, max_level);

  for (unsigned int b = 0; b < n_blocks; ++b)
    {
      for (unsigned int l = min_level; l <= max_level; ++l)
        dst_non_block[l].reinit(dst[l].block(b));
      const unsigned int data_block = same_for_all ? 0 : b;
      matrix_free_transfer_vector[data_block].copy_to_mg(*mg_dof[b], dst_non_block, src.block(b));

      for (unsigned int l = min_level; l <= max_level; ++l)
        dst[l].block(b) = dst_non_block[l];
    }
}

template <int dim, typename Number>
template <typename Number2, int spacedim>
void
MGTransferBlockMatrixFree<dim,Number>::
copy_from_mg (const DoFHandler<dim,spacedim>                        &mg_dof,
              LinearAlgebra::distributed::BlockVector<Number2>      &dst,
              const MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &src) const
{
  AssertDimension(matrix_free_transfer_vector.size(),1);
  const std::vector<const DoFHandler<dim, spacedim>*> mg_dofs (dst.n_blocks(), &mg_dof);

  copy_from_mg(mg_dofs, dst, src);
}

template <int dim, typename Number>
template <typename Number2, int spacedim>
void
MGTransferBlockMatrixFree<dim,Number>::
copy_from_mg (const std::vector<const DoFHandler<dim,spacedim>*> &mg_dof,
              LinearAlgebra::distributed::BlockVector<Number2>         &dst,
              const MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &src) const
{
  const unsigned int n_blocks  = dst.n_blocks();
  AssertDimension(mg_dof.size(), n_blocks);

  if (n_blocks==0)
    return;

  const unsigned int min_level = src.min_level();
  const unsigned int max_level = src.max_level();

  for (unsigned int l = min_level; l <= max_level; ++l)
    AssertDimension(src[l].n_blocks(), dst.n_blocks());

  // FIXME: this a quite ugly as we need a temporary object:
  MGLevelObject<LinearAlgebra::distributed::Vector<Number>> src_non_block(min_level, max_level);

  for (unsigned int b = 0; b < n_blocks; ++b)
    {
      for (unsigned int l = min_level; l <= max_level; ++l)
        {
          src_non_block[l].reinit(src[l].block(b));
          src_non_block[l] = src[l].block(b);
        }
      const unsigned int data_block = same_for_all ? 0 : b;
      matrix_free_transfer_vector[data_block].copy_from_mg(*mg_dof[b], dst.block(b), src_non_block);
    }
}



#endif // DOXYGEN


DEAL_II_NAMESPACE_CLOSE

#endif