operators.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 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_matrix_free_operators_h
#define dealii_matrix_free_operators_h


#include <deal.II/base/exceptions.h>
#include <deal.II/base/vectorization.h>
#include <deal.II/base/subscriptor.h>

#include <deal.II/lac/diagonal_matrix.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/multigrid/mg_constrained_dofs.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/fe_evaluation.h>


DEAL_II_NAMESPACE_OPEN


namespace MatrixFreeOperators
{

  namespace
  {
    // workaroud for unifying non-block vector and block vector implementations
    // a non-block vector has one block and the only subblock is the vector itself
    template <typename VectorType>
    typename std::enable_if<IsBlockVector<VectorType>::value,
             unsigned int>::type
             n_blocks(const VectorType &vector)
    {
      return vector.n_blocks();
    }

    template <typename VectorType>
    typename std::enable_if<!IsBlockVector<VectorType>::value,
             unsigned int>::type
             n_blocks(const VectorType &)
    {
      return 1;
    }

    template <typename VectorType>
    typename std::enable_if<IsBlockVector<VectorType>::value,
             typename VectorType::BlockType &>::type
             subblock(VectorType &vector, unsigned int block_no)
    {
      return vector.block(block_no);
    }

    template <typename VectorType>
    typename std::enable_if<IsBlockVector<VectorType>::value,
             const typename VectorType::BlockType &>::type
             subblock(const VectorType &vector, unsigned int block_no)
    {
      AssertIndexRange (block_no, vector.n_blocks());
      return vector.block(block_no);
    }

    template <typename VectorType>
    typename std::enable_if<!IsBlockVector<VectorType>::value,
             VectorType &>::type
             subblock(VectorType &vector, unsigned int)
    {
      return vector;
    }

    template <typename VectorType>
    typename std::enable_if<!IsBlockVector<VectorType>::value,
             const VectorType &>::type
             subblock(const VectorType &vector, unsigned int)
    {
      return vector;
    }

    template <typename VectorType>
    typename std::enable_if<IsBlockVector<VectorType>::value,
             void>::type
             collect_sizes(VectorType &vector)
    {
      vector.collect_sizes();
    }

    template <typename VectorType>
    typename std::enable_if<!IsBlockVector<VectorType>::value,
             void>::type
             collect_sizes(const VectorType &)
    {}
  }

  template <int dim, typename VectorType = LinearAlgebra::distributed::Vector<double> >
  class Base : public Subscriptor
  {
  public:
    typedef typename VectorType::value_type value_type;

    typedef typename VectorType::size_type size_type;

    Base ();

    virtual ~Base() = default;

    virtual void clear();

    void initialize (std::shared_ptr<const MatrixFree<dim,value_type> > data,
                     const std::vector<unsigned int> &selected_row_blocks = std::vector<unsigned int>(),
                     const std::vector<unsigned int> &selected_column_blocks = std::vector<unsigned int>());

    void initialize (std::shared_ptr<const MatrixFree<dim,value_type> > data,
                     const MGConstrainedDoFs &mg_constrained_dofs,
                     const unsigned int level,
                     const std::vector<unsigned int> &selected_row_blocks = std::vector<unsigned int>());

    void
    initialize(std::shared_ptr<const MatrixFree<dim, value_type> > data_,
               const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
               const unsigned int level,
               const std::vector<unsigned int> &selected_row_blocks = std::vector<unsigned int>());

    size_type m () const;

    size_type n () const;

    void vmult_interface_down(VectorType &dst,
                              const VectorType &src) const;

    void vmult_interface_up(VectorType &dst,
                            const VectorType &src) const;

    void vmult (VectorType &dst,
                const VectorType &src) const;

    void Tvmult (VectorType &dst,
                 const VectorType &src) const;

    void vmult_add (VectorType &dst,
                    const VectorType &src) const;

    void Tvmult_add (VectorType &dst,
                     const VectorType &src) const;

    value_type el (const unsigned int row,
                   const unsigned int col) const;

    virtual std::size_t memory_consumption () const;

    void initialize_dof_vector (VectorType &vec) const;

    virtual void compute_diagonal () = 0;

    std::shared_ptr<const MatrixFree<dim,value_type> >
    get_matrix_free () const;

    const std::shared_ptr<DiagonalMatrix<VectorType> > &
    get_matrix_diagonal_inverse() const;

    const std::shared_ptr<DiagonalMatrix<VectorType> > &
    get_matrix_diagonal() const;


    void precondition_Jacobi(VectorType &dst,
                             const VectorType &src,
                             const value_type omega) const;

  protected:

    void preprocess_constraints(VectorType &dst, const VectorType &src) const;

    void postprocess_constraints(VectorType &dst, const VectorType &src) const;

    void set_constrained_entries_to_one (VectorType &dst) const;

    virtual void apply_add(VectorType &dst,
                           const VectorType &src) const = 0;

    virtual void Tapply_add(VectorType &dst,
                            const VectorType &src) const;

    std::shared_ptr<const MatrixFree<dim,value_type> > data;

    std::shared_ptr<DiagonalMatrix<VectorType > > diagonal_entries;

    std::shared_ptr<DiagonalMatrix<VectorType > > inverse_diagonal_entries;

    std::vector<unsigned int> selected_rows;

    std::vector<unsigned int> selected_columns;

  private:

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

    mutable std::vector<std::vector<std::pair<value_type,value_type> > >
    edge_constrained_values;

    bool have_interface_matrices;

    void mult_add (VectorType &dst,
                   const VectorType &src,
                   const bool transpose) const;

    void adjust_ghost_range_if_necessary(const VectorType &vec,
                                         const bool        is_row) const;
  };



  template <typename OperatorType>
  class MGInterfaceOperator : public Subscriptor
  {
  public:
    typedef typename OperatorType::value_type value_type;

    typedef typename OperatorType::size_type size_type;

    MGInterfaceOperator();

    void clear();

    void initialize (const OperatorType &operator_in);

    template <typename VectorType>
    void vmult (VectorType &dst,
                const VectorType &src) const;

    template <typename VectorType>
    void Tvmult (VectorType &dst,
                 const VectorType &src) const;

    template <typename VectorType>
    void initialize_dof_vector (VectorType &vec) const;


  private:
    SmartPointer<const OperatorType> mf_base_operator;
  };



  template <int dim, int fe_degree, int n_components = 1, typename Number = double>
  class CellwiseInverseMassMatrix
  {
  public:
    CellwiseInverseMassMatrix (const FEEvaluationBase<dim,n_components,Number> &fe_eval);

    void apply(const AlignedVector<VectorizedArray<Number> > &inverse_coefficient,
               const unsigned int             n_actual_components,
               const VectorizedArray<Number> *in_array,
               VectorizedArray<Number>       *out_array) const;

    void fill_inverse_JxW_values(AlignedVector<VectorizedArray<Number> > &inverse_jxw) const;

  private:
    const FEEvaluationBase<dim,n_components,Number> &fe_eval;

    AlignedVector<VectorizedArray<Number> > inverse_shape;
  };



  template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1, int n_components = 1, typename VectorType = LinearAlgebra::distributed::Vector<double> >
  class MassOperator : public Base<dim, VectorType>
  {
  public:
    typedef typename Base<dim,VectorType>::value_type value_type;

    typedef typename Base<dim,VectorType>::size_type size_type;

    MassOperator ();

    virtual void compute_diagonal ();

  private:
    virtual void apply_add (VectorType       &dst,
                            const VectorType &src) const;

    void local_apply_cell (const MatrixFree<dim,value_type>            &data,
                           VectorType                                  &dst,
                           const VectorType                            &src,
                           const std::pair<unsigned int,unsigned int>  &cell_range) const;
  };



  template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1, int n_components = 1, typename VectorType = LinearAlgebra::distributed::Vector<double> >
  class LaplaceOperator : public Base<dim, VectorType>
  {
  public:
    typedef typename Base<dim,VectorType>::value_type value_type;

    typedef typename Base<dim,VectorType>::size_type size_type;

    LaplaceOperator ();

    virtual void compute_diagonal ();

    void set_coefficient(const std::shared_ptr<Table<2, VectorizedArray<value_type> > > &scalar_coefficient );

    virtual void clear();

    std::shared_ptr< Table<2, VectorizedArray<value_type> > > get_coefficient();

  private:
    virtual void apply_add (VectorType       &dst,
                            const VectorType &src) const;

    void local_apply_cell (const MatrixFree<dim,value_type>            &data,
                           VectorType                                  &dst,
                           const VectorType                            &src,
                           const std::pair<unsigned int,unsigned int>  &cell_range) const;

    void local_diagonal_cell (const MatrixFree<dim,value_type>            &data,
                              VectorType                                  &dst,
                              const unsigned int &,
                              const std::pair<unsigned int,unsigned int>  &cell_range) const;

    void do_operation_on_cell(FEEvaluation<dim,fe_degree,n_q_points_1d,n_components,value_type> &phi,
                              const unsigned int cell) const;

    std::shared_ptr< Table<2, VectorizedArray<value_type> > > scalar_coefficient;
  };



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

  template <int dim, int fe_degree, int n_components, typename Number>
  inline
  CellwiseInverseMassMatrix<dim,fe_degree,n_components,Number>
  ::CellwiseInverseMassMatrix (const FEEvaluationBase<dim,n_components,Number> &fe_eval)
    :
    fe_eval (fe_eval)
  {
    FullMatrix<double> shapes_1d(fe_degree+1, fe_degree+1);
    for (unsigned int i=0, c=0; i<shapes_1d.m(); ++i)
      for (unsigned int j=0; j<shapes_1d.n(); ++j, ++c)
        shapes_1d(i,j) = fe_eval.get_shape_info().shape_values[c][0];
    shapes_1d.gauss_jordan();
    const unsigned int stride = (fe_degree+2)/2;
    inverse_shape.resize(stride*(fe_degree+1));
    for (unsigned int i=0; i<stride; ++i)
      for (unsigned int q=0; q<(fe_degree+2)/2; ++q)
        {
          inverse_shape[i*stride+q] =
            0.5 * (shapes_1d(i,q) + shapes_1d(i,fe_degree-q));
          inverse_shape[(fe_degree-i)*stride+q] =
            0.5 * (shapes_1d(i,q) - shapes_1d(i,fe_degree-q));
        }
    if (fe_degree % 2 == 0)
      for (unsigned int q=0; q<(fe_degree+2)/2; ++q)
        inverse_shape[fe_degree/2*stride+q] = shapes_1d(fe_degree/2,q);
  }



  template <int dim, int fe_degree, int n_components, typename Number>
  inline
  void
  CellwiseInverseMassMatrix<dim,fe_degree,n_components,Number>
  ::fill_inverse_JxW_values(AlignedVector<VectorizedArray<Number> > &inverse_jxw) const
  {
    constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim);
    Assert(inverse_jxw.size() > 0 &&
           inverse_jxw.size() % dofs_per_cell == 0,
           ExcMessage("Expected diagonal to be a multiple of scalar dof per cells"));

    // temporarily reduce size of inverse_jxw to dofs_per_cell to get JxW values
    // from fe_eval (will not reallocate any memory)
    const unsigned int previous_size = inverse_jxw.size();
    inverse_jxw.resize(dofs_per_cell);
    fe_eval.fill_JxW_values(inverse_jxw);

    // invert
    inverse_jxw.resize_fast(previous_size);
    for (unsigned int q=0; q<dofs_per_cell; ++q)
      inverse_jxw[q] = 1./inverse_jxw[q];
    // copy values to rest of vector
    for (unsigned int q=dofs_per_cell; q<previous_size; )
      for (unsigned int i=0; i<dofs_per_cell; ++i, ++q)
        inverse_jxw[q] = inverse_jxw[i];
  }



  template <int dim, int fe_degree, int n_components, typename Number>
  inline
  void
  CellwiseInverseMassMatrix<dim,fe_degree,n_components,Number>
  ::apply(const AlignedVector<VectorizedArray<Number> > &inverse_coefficients,
          const unsigned int             n_actual_components,
          const VectorizedArray<Number> *in_array,
          VectorizedArray<Number>       *out_array) const
  {
    constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim);
    Assert(inverse_coefficients.size() > 0 &&
           inverse_coefficients.size() % dofs_per_cell == 0,
           ExcMessage("Expected diagonal to be a multiple of scalar dof per cells"));
    if (inverse_coefficients.size() != dofs_per_cell)
      AssertDimension(n_actual_components * dofs_per_cell, inverse_coefficients.size());

    Assert(dim == 2 || dim == 3, ExcNotImplemented());

    internal::EvaluatorTensorProduct<internal::evaluate_evenodd,dim,fe_degree+1,
             fe_degree+1, VectorizedArray<Number> >
             evaluator(inverse_shape, inverse_shape, inverse_shape);

    const unsigned int shift_coefficient =
      inverse_coefficients.size() > dofs_per_cell ? dofs_per_cell : 0;
    const VectorizedArray<Number> *inv_coefficient = &inverse_coefficients[0];
    VectorizedArray<Number> temp_data_field[dofs_per_cell];
    for (unsigned int d=0; d<n_actual_components; ++d)
      {
        const VectorizedArray<Number> *in = in_array+d*dofs_per_cell;
        VectorizedArray<Number> *out = out_array+d*dofs_per_cell;
        // Need to select 'apply' method with hessian slot because values
        // assume symmetries that do not exist in the inverse shapes
        evaluator.template hessians<0,false,false> (in, temp_data_field);
        evaluator.template hessians<1,false,false> (temp_data_field, out);

        if (dim == 3)
          {
            evaluator.template hessians<2,false,false> (out, temp_data_field);
            for (unsigned int q=0; q<dofs_per_cell; ++q)
              temp_data_field[q] *= inv_coefficient[q];
            evaluator.template hessians<2,true,false> (temp_data_field, out);
          }
        else if (dim == 2)
          for (unsigned int q=0; q<dofs_per_cell; ++q)
            out[q] *= inv_coefficient[q];

        evaluator.template hessians<1,true,false>(out, temp_data_field);
        evaluator.template hessians<0,true,false>(temp_data_field, out);

        inv_coefficient += shift_coefficient;
      }
  }

  //----------------- Base operator -----------------------------
  template <int dim, typename VectorType>
  Base<dim,VectorType>::Base ()
    :
    Subscriptor(),
    have_interface_matrices(false)
  {}



  template <int dim, typename VectorType>
  typename Base<dim,VectorType>::size_type
  Base<dim,VectorType>::m () const
  {
    Assert(data.get() != nullptr,
           ExcNotInitialized());
    typename Base<dim, VectorType>::size_type total_size = 0;
    for (unsigned int i=0; i<selected_rows.size(); ++i)
      total_size += data->get_vector_partitioner(selected_rows[i])->size();
    return total_size;
  }



  template <int dim, typename VectorType>
  typename Base<dim,VectorType>::size_type
  Base<dim,VectorType>::n () const
  {
    Assert(data.get() != nullptr,
           ExcNotInitialized());
    typename Base<dim, VectorType>::size_type total_size = 0;
    for (unsigned int i=0; i<selected_columns.size(); ++i)
      total_size += data->get_vector_partitioner(selected_columns[i])->size();
    return total_size;
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::clear ()
  {
    data.reset();
    inverse_diagonal_entries.reset();
  }



  template <int dim, typename VectorType>
  typename Base<dim,VectorType>::value_type
  Base<dim,VectorType>::el (const unsigned int row,
                            const unsigned int col) const
  {
    (void) col;
    Assert (row == col, ExcNotImplemented());
    Assert (inverse_diagonal_entries.get() != nullptr &&
            inverse_diagonal_entries->m() > 0, ExcNotInitialized());
    return 1.0/(*inverse_diagonal_entries)(row,row);
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::initialize_dof_vector (VectorType &vec) const
  {
    Assert(data.get() != nullptr,
           ExcNotInitialized());
    AssertDimension(n_blocks(vec), selected_rows.size());
    for (unsigned int i = 0; i < n_blocks(vec); ++i)
      {
        const unsigned int index = selected_rows[i];
        if (!subblock(vec,index).partitioners_are_compatible
            (*data->get_dof_info(index).vector_partitioner))
          data->initialize_dof_vector(subblock(vec,index),index);

        Assert(subblock(vec,index).partitioners_are_globally_compatible
               (*data->get_dof_info(index).vector_partitioner),
               ExcInternalError());
      }
    collect_sizes(vec);
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::
  initialize (std::shared_ptr<const MatrixFree<dim,typename Base<dim,VectorType>::value_type> > data_,
              const std::vector<unsigned int> &given_row_selection,
              const std::vector<unsigned int> &given_column_selection)
  {
    data = data_;

    selected_rows.clear();
    selected_columns.clear();
    if (given_row_selection.empty())
      for (unsigned int i=0; i<data_->n_components(); ++i)
        selected_rows.push_back(i);
    else
      {
        for (unsigned int i=0; i<given_row_selection.size(); ++i)
          {
            AssertIndexRange(given_row_selection[i], data_->n_components());
            for (unsigned int j=0; j<given_row_selection.size(); ++j)
              if (j!=i)
                Assert(given_row_selection[j] != given_row_selection[i],
                       ExcMessage("Given row indices must be unique"));

            selected_rows.push_back(given_row_selection[i]);
          }
      }
    if (given_column_selection.size() == 0)
      selected_columns = selected_rows;
    else
      {
        for (unsigned int i=0; i<given_column_selection.size(); ++i)
          {
            AssertIndexRange(given_column_selection[i], data_->n_components());
            for (unsigned int j=0; j<given_column_selection.size(); ++j)
              if (j!=i)
                Assert(given_column_selection[j] != given_column_selection[i],
                       ExcMessage("Given column indices must be unique"));

            selected_columns.push_back(given_column_selection[i]);
          }
      }

    edge_constrained_indices.clear();
    edge_constrained_indices.resize(selected_rows.size());
    edge_constrained_values.clear();
    edge_constrained_values.resize(selected_rows.size());
    have_interface_matrices = false;
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::
  initialize (std::shared_ptr<const MatrixFree<dim,typename Base<dim,VectorType>::value_type> > data_,
              const MGConstrainedDoFs         &mg_constrained_dofs,
              const unsigned int               level,
              const std::vector<unsigned int> &given_row_selection)
  {
    std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(1, mg_constrained_dofs);
    initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
  }



  template <int dim, typename VectorType>
  void
  Base<dim, VectorType>::
  initialize (std::shared_ptr<const MatrixFree<dim, value_type> > data_,
              const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
              const unsigned int                    level,
              const std::vector<unsigned int>      &given_row_selection)
  {
    AssertThrow (level != numbers::invalid_unsigned_int,
                 ExcMessage("level is not set"));

    selected_rows.clear();
    selected_columns.clear();
    if (given_row_selection.empty())
      for (unsigned int i=0; i<data_->n_components(); ++i)
        selected_rows.push_back(i);
    else
      {
        for (unsigned int i=0; i<given_row_selection.size(); ++i)
          {
            AssertIndexRange(given_row_selection[i], data_->n_components());
            for (unsigned int j=0; j<given_row_selection.size(); ++j)
              if (j!=i)
                Assert(given_row_selection[j] != given_row_selection[i],
                       ExcMessage("Given row indices must be unique"));

            selected_rows.push_back(given_row_selection[i]);
          }
      }
    selected_columns = selected_rows;

    AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
    edge_constrained_indices.clear();
    edge_constrained_indices.resize(selected_rows.size());
    edge_constrained_values.clear();
    edge_constrained_values.resize(selected_rows.size());

    data = data_;

    for (unsigned int j = 0; j < selected_rows.size(); ++j)
      {
        if (data_->n_macro_cells() > 0)
          {
            AssertDimension(static_cast<int>(level),
                            data_->get_cell_iterator(0,0,j)->level());
          }

        // setup edge_constrained indices
        std::vector<types::global_dof_index> interface_indices;
        mg_constrained_dofs[j].get_refinement_edge_indices(level)
        .fill_index_vector(interface_indices);
        edge_constrained_indices[j].clear();
        edge_constrained_indices[j].reserve(interface_indices.size());
        edge_constrained_values[j].resize(interface_indices.size());
        const IndexSet &locally_owned = data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
        for (unsigned int i=0; i<interface_indices.size(); ++i)
          if (locally_owned.is_element(interface_indices[i]))
            edge_constrained_indices[j].push_back
            (locally_owned.index_within_set(interface_indices[i]));
        have_interface_matrices |=
          Utilities::MPI::max(static_cast<unsigned int>(edge_constrained_indices[j].size()),
                              data->get_vector_partitioner()->get_mpi_communicator()) > 0;
      }
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::set_constrained_entries_to_one (VectorType &dst) const
  {
    for (unsigned int j=0; j<n_blocks(dst); ++j)
      {
        const std::vector<unsigned int> &
        constrained_dofs = data->get_constrained_dofs(selected_rows[j]);
        for (unsigned int i=0; i<constrained_dofs.size(); ++i)
          subblock(dst,j).local_element(constrained_dofs[i]) = 1.;
        for (unsigned int i=0; i<edge_constrained_indices[j].size(); ++i)
          subblock(dst,j).local_element(edge_constrained_indices[j][i]) = 1.;
      }
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::vmult (VectorType       &dst,
                               const VectorType &src) const
  {
    typedef typename Base<dim,VectorType>::value_type Number;
    dst = Number(0.);
    vmult_add (dst, src);
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::vmult_add (VectorType &dst,
                                   const VectorType &src) const
  {
    mult_add (dst, src, false);
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::Tvmult_add (VectorType &dst,
                                    const VectorType &src) const
  {
    mult_add (dst, src, true);
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::adjust_ghost_range_if_necessary(const VectorType &src,
                                                        const bool is_row) const
  {
    typedef typename Base<dim,VectorType>::value_type Number;
    for (unsigned int i=0; i<n_blocks(src); ++i)
      {
        const unsigned int mf_component = is_row ? selected_rows[i] : selected_columns[i];
        // If both vectors use the same partitioner -> done
        if (subblock(src,i).get_partitioner().get() ==
            data->get_dof_info(mf_component).vector_partitioner.get())
          continue;

        // If not, assert that the local ranges are the same and reset to the
        // current partitioner
        Assert(subblock(src,i).get_partitioner()->local_size() ==
               data->get_dof_info(mf_component).vector_partitioner->local_size(),
               ExcMessage("The vector passed to the vmult() function does not have "
                          "the correct size for compatibility with MatrixFree."));

        // copy the vector content to a temporary vector so that it does not get
        // lost
        LinearAlgebra::distributed::Vector<Number> copy_vec(subblock(src,i));
        subblock(const_cast<VectorType &>(src),i).
        reinit(data->get_dof_info(mf_component).vector_partitioner);
        subblock(const_cast<VectorType &>(src),i).copy_locally_owned_data_from(copy_vec);
      }
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::preprocess_constraints (VectorType &dst,
                                                const VectorType &src) const
  {
    typedef typename Base<dim,VectorType>::value_type Number;
    adjust_ghost_range_if_necessary(src, false);
    adjust_ghost_range_if_necessary(dst, true);

    // set zero Dirichlet values on the input vector (and remember the src and
    // dst values because we need to reset them at the end)
    for (unsigned int j=0; j<n_blocks(dst); ++j)
      {
        for (unsigned int i=0; i<edge_constrained_indices[j].size(); ++i)
          {
            edge_constrained_values[j][i] =
              std::pair<Number,Number>
              (subblock(src,j).local_element(edge_constrained_indices[j][i]),
               subblock(dst,j).local_element(edge_constrained_indices[j][i]));
            subblock(const_cast<VectorType &>(src),j).local_element(edge_constrained_indices[j][i]) = 0.;
          }
      }
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::mult_add (VectorType &dst,
                                  const VectorType &src,
                                  const bool transpose) const
  {
    AssertDimension(dst.size(), src.size());
    AssertDimension(n_blocks(dst), n_blocks(src));
    AssertDimension(n_blocks(dst), selected_rows.size());
    preprocess_constraints (dst,src);
    if (transpose)
      Tapply_add(dst,src);
    else
      apply_add(dst,src);
    postprocess_constraints(dst,src);
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::postprocess_constraints (VectorType &dst,
                                                 const VectorType &src) const
  {
    for (unsigned int j=0; j<n_blocks(dst); ++j)
      {
        const std::vector<unsigned int> &
        constrained_dofs = data->get_constrained_dofs(selected_rows[j]);
        for (unsigned int i=0; i<constrained_dofs.size(); ++i)
          subblock(dst,j).local_element(constrained_dofs[i]) +=
            subblock(src,j).local_element(constrained_dofs[i]);
      }

    // reset edge constrained values, multiply by unit matrix and add into
    // destination
    for (unsigned int j=0; j<n_blocks(dst); ++j)
      {
        for (unsigned int i=0; i<edge_constrained_indices[j].size(); ++i)
          {
            subblock(const_cast<VectorType &>(src),j).local_element
            (edge_constrained_indices[j][i])
              = edge_constrained_values[j][i].first;
            subblock(dst,j).local_element(edge_constrained_indices[j][i])
              = edge_constrained_values[j][i].second +
                edge_constrained_values[j][i].first;
          }
      }
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::
  vmult_interface_down(VectorType &dst,
                       const VectorType &src) const
  {
    typedef typename Base<dim,VectorType>::value_type Number;
    AssertDimension(dst.size(), src.size());
    adjust_ghost_range_if_necessary(src, false);
    adjust_ghost_range_if_necessary(dst, true);

    dst = Number(0.);

    if (!have_interface_matrices)
      return;

    // set zero Dirichlet values on the input vector (and remember the src and
    // dst values because we need to reset them at the end)
    for (unsigned int j=0; j<n_blocks(dst); ++j)
      for (unsigned int i=0; i<edge_constrained_indices[j].size(); ++i)
        {
          edge_constrained_values[j][i] =
            std::pair<Number,Number>
            (subblock(src,j).local_element(edge_constrained_indices[j][i]),
             subblock(dst,j).local_element(edge_constrained_indices[j][i]));
          subblock(const_cast<VectorType &>(src),j)
          .local_element(edge_constrained_indices[j][i]) = 0.;
        }

    apply_add(dst,src);

    for (unsigned int j=0; j<n_blocks(dst); ++j)
      {
        unsigned int c=0;
        for (unsigned int i=0; i<edge_constrained_indices[j].size(); ++i)
          {
            for ( ; c<edge_constrained_indices[j][i]; ++c)
              subblock(dst,j).local_element(c) = 0.;
            ++c;

            // reset the src values
            subblock(const_cast<VectorType &>(src),j)
            .local_element(edge_constrained_indices[j][i])
              = edge_constrained_values[j][i].first;
          }
        for ( ; c<subblock(dst,j).local_size(); ++c)
          subblock(dst,j).local_element(c) = 0.;
      }
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::
  vmult_interface_up(VectorType &dst,
                     const VectorType &src) const
  {
    typedef typename Base<dim,VectorType>::value_type Number;
    AssertDimension(dst.size(), src.size());
    adjust_ghost_range_if_necessary(src, false);
    adjust_ghost_range_if_necessary(dst, true);

    dst = Number(0.);

    if (!have_interface_matrices)
      return;

    VectorType src_cpy (src);
    for (unsigned int j=0; j<n_blocks(dst); ++j)
      {
        unsigned int c=0;
        for (unsigned int i=0; i<edge_constrained_indices[j].size(); ++i)
          {
            for ( ; c<edge_constrained_indices[j][i]; ++c)
              subblock(src_cpy,j).local_element(c) = 0.;
            ++c;
          }
        for ( ; c<subblock(src_cpy,j).local_size(); ++c)
          subblock(src_cpy,j).local_element(c) = 0.;
      }

    apply_add(dst,src_cpy);

    for (unsigned int j=0; j<n_blocks(dst); ++j)
      for (unsigned int i=0; i<edge_constrained_indices[j].size(); ++i)
        subblock(dst,j).local_element(edge_constrained_indices[j][i]) = 0.;
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::Tvmult (VectorType       &dst,
                                const VectorType &src) const
  {
    typedef typename Base<dim,VectorType>::value_type Number;
    dst = Number(0.);
    Tvmult_add (dst,src);
  }



  template <int dim, typename VectorType>
  std::size_t
  Base<dim,VectorType>::memory_consumption () const
  {
    return inverse_diagonal_entries.get() != nullptr ?
           inverse_diagonal_entries->memory_consumption() : sizeof(*this);
  }



  template <int dim, typename VectorType>
  std::shared_ptr<const MatrixFree<dim,typename Base<dim,VectorType>::value_type> >
  Base<dim,VectorType>::get_matrix_free() const
  {
    return data;
  }



  template <int dim, typename VectorType>
  const std::shared_ptr<DiagonalMatrix<VectorType> > &
  Base<dim,VectorType>::get_matrix_diagonal_inverse() const
  {
    Assert(inverse_diagonal_entries.get() != nullptr &&
           inverse_diagonal_entries->m() > 0, ExcNotInitialized());
    return inverse_diagonal_entries;
  }



  template <int dim, typename VectorType>
  const std::shared_ptr<DiagonalMatrix<VectorType> > &
  Base<dim,VectorType>::get_matrix_diagonal() const
  {
    Assert(diagonal_entries.get() != nullptr &&
           diagonal_entries->m() > 0, ExcNotInitialized());
    return diagonal_entries;
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::Tapply_add(VectorType &dst,
                                   const VectorType &src) const
  {
    apply_add(dst,src);
  }



  template <int dim, typename VectorType>
  void
  Base<dim,VectorType>::precondition_Jacobi(VectorType &dst,
                                            const VectorType &src,
                                            const typename Base<dim,VectorType>::value_type omega) const
  {
    Assert(inverse_diagonal_entries.get() &&
           inverse_diagonal_entries->m() > 0, ExcNotInitialized());
    inverse_diagonal_entries->vmult(dst,src);
    dst *= omega;
  }



  //------------------------- MGInterfaceOperator ------------------------------

  template <typename OperatorType>
  MGInterfaceOperator<OperatorType>::MGInterfaceOperator ()
    :
    Subscriptor(),
    mf_base_operator(nullptr)
  {
  }



  template <typename OperatorType>
  void
  MGInterfaceOperator<OperatorType>::clear ()
  {
    mf_base_operator = nullptr;
  }



  template <typename OperatorType>
  void
  MGInterfaceOperator<OperatorType>::initialize (const OperatorType &operator_in)
  {
    mf_base_operator = &operator_in;
  }



  template <typename OperatorType>
  template <typename VectorType>
  void
  MGInterfaceOperator<OperatorType>::vmult (VectorType &dst,
                                            const VectorType &src) const
  {
#ifndef DEAL_II_MSVC
    static_assert (std::is_same<typename VectorType::value_type,value_type>::value,
                   "The vector type must be based on the same value type as this"
                   "operator");
#endif

    Assert(mf_base_operator != nullptr,
           ExcNotInitialized());

    mf_base_operator->vmult_interface_down(dst, src);
  }



  template <typename OperatorType>
  template <typename VectorType>
  void
  MGInterfaceOperator<OperatorType>::Tvmult (VectorType &dst,
                                             const VectorType &src) const
  {
#ifndef DEAL_II_MSVC
    static_assert (std::is_same<typename VectorType::value_type,value_type>::value,
                   "The vector type must be based on the same value type as this"
                   "operator");
#endif

    Assert(mf_base_operator != nullptr,
           ExcNotInitialized());

    mf_base_operator->vmult_interface_up(dst, src);
  }



  template <typename OperatorType>
  template <typename VectorType>
  void
  MGInterfaceOperator<OperatorType>::initialize_dof_vector (VectorType &vec) const
  {
    Assert(mf_base_operator != nullptr,
           ExcNotInitialized());

    mf_base_operator->initialize_dof_vector(vec);
  }



  //-----------------------------MassOperator----------------------------------

  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename VectorType>
  MassOperator<dim, fe_degree, n_q_points_1d, n_components, VectorType>::
  MassOperator ()
    :
    Base<dim, VectorType>()
  {}



  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename VectorType>
  void
  MassOperator<dim, fe_degree, n_q_points_1d, n_components, VectorType>::
  compute_diagonal()
  {
    typedef typename Base<dim,VectorType>::value_type Number;
    Assert((Base<dim, VectorType>::data.get() != nullptr), ExcNotInitialized());

    this->inverse_diagonal_entries.
    reset(new DiagonalMatrix<VectorType>());
    this->diagonal_entries.
    reset(new DiagonalMatrix<VectorType>());
    VectorType &inverse_diagonal_vector = this->inverse_diagonal_entries->get_vector();
    VectorType &diagonal_vector         = this->diagonal_entries->get_vector();
    this->initialize_dof_vector(inverse_diagonal_vector);
    this->initialize_dof_vector(diagonal_vector);
    inverse_diagonal_vector = Number(1.);
    apply_add(diagonal_vector, inverse_diagonal_vector);

    this->set_constrained_entries_to_one(diagonal_vector);
    inverse_diagonal_vector = diagonal_vector;

    const unsigned int local_size = inverse_diagonal_vector.local_size();
    for (unsigned int i=0; i<local_size; ++i)
      inverse_diagonal_vector.local_element(i)
        =1./inverse_diagonal_vector.local_element(i);

    inverse_diagonal_vector.update_ghost_values();
    diagonal_vector.update_ghost_values();
  }



  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename VectorType>
  void
  MassOperator<dim, fe_degree, n_q_points_1d, n_components, VectorType>::
  apply_add (VectorType       &dst,
             const VectorType &src) const
  {
    Base<dim, VectorType>::data->cell_loop (&MassOperator::local_apply_cell,
                                            this, dst, src);
  }



  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename VectorType>
  void
  MassOperator<dim, fe_degree, n_q_points_1d, n_components, VectorType>::
  local_apply_cell (const MatrixFree<dim,typename Base<dim,VectorType>::value_type> &data,
                    VectorType                                  &dst,
                    const VectorType                            &src,
                    const std::pair<unsigned int,unsigned int>  &cell_range) const
  {
    typedef typename Base<dim,VectorType>::value_type Number;
    FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number> phi(data, this->selected_rows[0]);
    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
      {
        phi.reinit (cell);
        phi.read_dof_values(src);
        phi.evaluate (true,false,false);
        for (unsigned int q=0; q<phi.n_q_points; ++q)
          phi.submit_value (phi.get_value(q), q);
        phi.integrate (true,false);
        phi.distribute_local_to_global (dst);
      }
  }


  //-----------------------------LaplaceOperator----------------------------------

  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename VectorType>
  LaplaceOperator<dim, fe_degree, n_q_points_1d, n_components, VectorType>::
  LaplaceOperator ()
    :
    Base<dim, VectorType>()
  {
  }



  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename VectorType>
  void
  LaplaceOperator<dim, fe_degree, n_q_points_1d, n_components, VectorType>::
  clear ()
  {
    Base<dim, VectorType>::clear();
    scalar_coefficient.reset();
  }



  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename VectorType>
  void
  LaplaceOperator<dim, fe_degree, n_q_points_1d, n_components, VectorType>::
  set_coefficient(const std::shared_ptr<Table<2, VectorizedArray<typename Base<dim,VectorType>::value_type> > > &scalar_coefficient_ )
  {
    scalar_coefficient = scalar_coefficient_;
  }



  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename VectorType>
  std::shared_ptr< Table<2, VectorizedArray< typename LaplaceOperator<dim, fe_degree, n_q_points_1d, n_components, VectorType>::value_type> > >
  LaplaceOperator<dim, fe_degree, n_q_points_1d, n_components, VectorType>::
  get_coefficient()
  {
    Assert (scalar_coefficient.get(),
            ExcNotInitialized());
    return scalar_coefficient;
  }



  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename VectorType>
  void
  LaplaceOperator<dim, fe_degree, n_q_points_1d, n_components, VectorType>::
  compute_diagonal()
  {
    typedef typename Base<dim,VectorType>::value_type Number;
    Assert((Base<dim, VectorType>::data.get() != nullptr), ExcNotInitialized());

    unsigned int dummy = 0;
    this->inverse_diagonal_entries.
    reset(new DiagonalMatrix<VectorType>());
    this->diagonal_entries.
    reset(new DiagonalMatrix<VectorType>());
    VectorType &inverse_diagonal_vector = this->inverse_diagonal_entries->get_vector();
    VectorType &diagonal_vector         = this->diagonal_entries->get_vector();
    this->initialize_dof_vector(inverse_diagonal_vector);
    this->initialize_dof_vector(diagonal_vector);

    this->data->cell_loop (&LaplaceOperator::local_diagonal_cell,
                           this, diagonal_vector, dummy);
    this->set_constrained_entries_to_one(diagonal_vector);

    inverse_diagonal_vector = diagonal_vector;

    for (unsigned int i=0; i<inverse_diagonal_vector.local_size(); ++i)
      if (std::abs(inverse_diagonal_vector.local_element(i)) > std::sqrt(std::numeric_limits<Number>::epsilon()))
        inverse_diagonal_vector.local_element(i) = 1./inverse_diagonal_vector.local_element(i);
      else
        inverse_diagonal_vector.local_element(i) = 1.;

    inverse_diagonal_vector.update_ghost_values();
    diagonal_vector.update_ghost_values();
  }



  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename VectorType>
  void
  LaplaceOperator<dim, fe_degree, n_q_points_1d, n_components, VectorType>::
  apply_add (VectorType       &dst,
             const VectorType &src) const
  {
    Base<dim, VectorType>::data->cell_loop (&LaplaceOperator::local_apply_cell,
                                            this, dst, src);
  }

  namespace
  {
    template <typename Number>
    bool
    non_negative(const VectorizedArray<Number> &n)
    {
      for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
        if (n[v] < 0.)
          return false;

      return true;
    }
  }



  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename VectorType>
  void
  LaplaceOperator<dim, fe_degree, n_q_points_1d, n_components, VectorType>::
  do_operation_on_cell(FEEvaluation<dim,fe_degree,n_q_points_1d,n_components,typename Base<dim,VectorType>::value_type> &phi,
                       const unsigned int cell) const
  {
    phi.evaluate (false,true,false);
    if (scalar_coefficient.get())
      {
        for (unsigned int q=0; q<phi.n_q_points; ++q)
          {
            Assert (non_negative((*scalar_coefficient)(cell,q)),
                    ExcMessage("Coefficient must be non-negative"));
            phi.submit_gradient ((*scalar_coefficient)(cell,q)*phi.get_gradient(q), q);
          }
      }
    else
      {
        for (unsigned int q=0; q<phi.n_q_points; ++q)
          {
            phi.submit_gradient (phi.get_gradient(q), q);
          }
      }
    phi.integrate (false,true);
  }




  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename VectorType>
  void
  LaplaceOperator<dim, fe_degree, n_q_points_1d, n_components, VectorType>::
  local_apply_cell (const MatrixFree<dim,typename Base<dim,VectorType>::value_type> &data,
                    VectorType       &dst,
                    const VectorType &src,
                    const std::pair<unsigned int,unsigned int>  &cell_range) const
  {
    typedef typename Base<dim,VectorType>::value_type Number;
    FEEvaluation<dim,fe_degree,n_q_points_1d,n_components,Number> phi (data, this->selected_rows[0]);
    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
      {
        phi.reinit (cell);
        phi.read_dof_values(src);
        do_operation_on_cell(phi,cell);
        phi.distribute_local_to_global (dst);
      }
  }


  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename VectorType>
  void
  LaplaceOperator<dim, fe_degree, n_q_points_1d, n_components, VectorType>::
  local_diagonal_cell (const MatrixFree<dim,typename Base<dim,VectorType>::value_type> &data,
                       VectorType                                       &dst,
                       const unsigned int &,
                       const std::pair<unsigned int,unsigned int>       &cell_range) const
  {
    typedef typename Base<dim,VectorType>::value_type Number;
    FEEvaluation<dim,fe_degree,n_q_points_1d,n_components,Number> phi (data, this->selected_rows[0]);
    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
      {
        phi.reinit (cell);
        VectorizedArray<Number> local_diagonal_vector[phi.static_dofs_per_cell];
        for (unsigned int i=0; i<phi.dofs_per_component; ++i)
          {
            for (unsigned int j=0; j<phi.dofs_per_component; ++j)
              phi.begin_dof_values()[j] = VectorizedArray<Number>();
            phi.begin_dof_values()[i] = 1.;
            do_operation_on_cell(phi,cell);
            local_diagonal_vector[i] = phi.begin_dof_values()[i];
          }
        for (unsigned int i=0; i<phi.dofs_per_component; ++i)
          for (unsigned int c=0; c<phi.n_components; ++c)
            phi.begin_dof_values()[i+c*phi.dofs_per_component] = local_diagonal_vector[i];
        phi.distribute_local_to_global (dst);
      }
  }


} // end of namespace MatrixFreeOperators


DEAL_II_NAMESPACE_CLOSE

#endif