operators.h

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


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

#include <deal.II/lac/diagonal_matrix.h>
#include <deal.II/lac/la_parallel_vector.h>

#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>

#include <deal.II/multigrid/mg_constrained_dofs.h>


DEAL_II_NAMESPACE_OPEN


namespace MatrixFreeOperators
{
  namespace BlockHelper
  {
    // 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 &)
    {}
  } // namespace BlockHelper

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

    using size_type = typename VectorType::size_type;

    Base();

    virtual ~Base() override = 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:
    using value_type = typename OperatorType::value_type;

    using size_type = typename OperatorType::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
    transform_from_q_points_to_basis(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:
    using value_type = typename Base<dim, VectorType>::value_type;

    using size_type = typename Base<dim, VectorType>::size_type;

    MassOperator();

    virtual void
    compute_diagonal() override;

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

    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:
    using value_type = typename Base<dim, VectorType>::value_type;

    using size_type = typename Base<dim, VectorType>::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 VectorType &,
      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_component_on_cell =
      Utilities::pow(fe_degree + 1, dim);
    Assert(inverse_jxw.size() > 0 &&
             inverse_jxw.size() % dofs_per_component_on_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)
    for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
      inverse_jxw[q] = 1. / fe_eval.JxW(q);
    // copy values to rest of vector
    for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
      for (unsigned int i = 0; i < dofs_per_component_on_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_component =
      Utilities::pow(fe_degree + 1, dim);
    Assert(inverse_coefficients.size() > 0 &&
             inverse_coefficients.size() % dofs_per_component == 0,
           ExcMessage(
             "Expected diagonal to be a multiple of scalar dof per cells"));
    if (inverse_coefficients.size() != dofs_per_component)
      AssertDimension(n_actual_components * dofs_per_component,
                      inverse_coefficients.size());

    Assert(dim >= 1 || 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_component ? dofs_per_component : 0;
    const VectorizedArray<Number> *inv_coefficient =
      inverse_coefficients.data();
    VectorizedArray<Number> temp_data_field[dofs_per_component];
    for (unsigned int d = 0; d < n_actual_components; ++d)
      {
        const VectorizedArray<Number> *in  = in_array + d * dofs_per_component;
        VectorizedArray<Number> *      out = out_array + d * dofs_per_component;
        // 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);
        if (dim > 1)
          {
            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_component; ++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_component; ++q)
                out[q] *= inv_coefficient[q];

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

        inv_coefficient += shift_coefficient;
      }
  }



  template <int dim, int fe_degree, int n_components, typename Number>
  inline void
  CellwiseInverseMassMatrix<dim, fe_degree, n_components, Number>::
    transform_from_q_points_to_basis(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);
    internal::EvaluatorTensorProduct<internal::evaluate_evenodd,
                                     dim,
                                     fe_degree + 1,
                                     fe_degree + 1,
                                     VectorizedArray<Number>>
      evaluator(AlignedVector<VectorizedArray<Number>>(),
                AlignedVector<VectorizedArray<Number>>(),
                inverse_shape);

    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;

        if (dim == 3)
          {
            evaluator.template hessians<2, true, false>(in, out);
            evaluator.template hessians<1, true, false>(out, out);
            evaluator.template hessians<0, true, false>(out, out);
          }
        if (dim == 2)
          {
            evaluator.template hessians<1, true, false>(in, out);
            evaluator.template hessians<0, true, false>(out, out);
          }
        if (dim == 1)
          evaluator.template hessians<0, true, false>(in, out);
      }
  }



  //----------------- 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(BlockHelper::n_blocks(vec), selected_rows.size());
    for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
      {
        const unsigned int index = selected_rows[i];
        if (!BlockHelper::subblock(vec, index)
               .partitioners_are_compatible(
                 *data->get_dof_info(index).vector_partitioner))
          data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);

        Assert(BlockHelper::subblock(vec, index)
                 .partitioners_are_globally_compatible(
                   *data->get_dof_info(index).vector_partitioner),
               ExcInternalError());
      }
    BlockHelper::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(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 (const auto interface_index : interface_indices)
          if (locally_owned.is_element(interface_index))
            edge_constrained_indices[j].push_back(
              locally_owned.index_within_set(interface_index));
        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 < BlockHelper::n_blocks(dst); ++j)
      {
        const std::vector<unsigned int> &constrained_dofs =
          data->get_constrained_dofs(selected_rows[j]);
        for (const auto constrained_dof : constrained_dofs)
          BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
        for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
          BlockHelper::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
  {
    using Number = typename Base<dim, VectorType>::value_type;
    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
  {
    using Number = typename Base<dim, VectorType>::value_type;
    for (unsigned int i = 0; i < BlockHelper::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 (BlockHelper::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(
          BlockHelper::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(
          BlockHelper::subblock(src, i));
        BlockHelper::subblock(const_cast<VectorType &>(src), i)
          .reinit(data->get_dof_info(mf_component).vector_partitioner);
        BlockHelper::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
  {
    using Number = typename Base<dim, VectorType>::value_type;
    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 < BlockHelper::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>(
              BlockHelper::subblock(src, j).local_element(
                edge_constrained_indices[j][i]),
              BlockHelper::subblock(dst, j).local_element(
                edge_constrained_indices[j][i]));
            BlockHelper::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(BlockHelper::n_blocks(dst), BlockHelper::n_blocks(src));
    AssertDimension(BlockHelper::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 < BlockHelper::n_blocks(dst); ++j)
      {
        const std::vector<unsigned int> &constrained_dofs =
          data->get_constrained_dofs(selected_rows[j]);
        for (const auto constrained_dof : constrained_dofs)
          BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
            BlockHelper::subblock(src, j).local_element(constrained_dof);
      }

    // reset edge constrained values, multiply by unit matrix and add into
    // destination
    for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
      {
        for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
          {
            BlockHelper::subblock(const_cast<VectorType &>(src), j)
              .local_element(edge_constrained_indices[j][i]) =
              edge_constrained_values[j][i].first;
            BlockHelper::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
  {
    using Number = typename Base<dim, VectorType>::value_type;
    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 < BlockHelper::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>(
            BlockHelper::subblock(src, j).local_element(
              edge_constrained_indices[j][i]),
            BlockHelper::subblock(dst, j).local_element(
              edge_constrained_indices[j][i]));
          BlockHelper::subblock(const_cast<VectorType &>(src), j)
            .local_element(edge_constrained_indices[j][i]) = 0.;
        }

    apply_add(dst, src);

    for (unsigned int j = 0; j < BlockHelper::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)
              BlockHelper::subblock(dst, j).local_element(c) = 0.;
            ++c;

            // reset the src values
            BlockHelper::subblock(const_cast<VectorType &>(src), j)
              .local_element(edge_constrained_indices[j][i]) =
              edge_constrained_values[j][i].first;
          }
        for (; c < BlockHelper::subblock(dst, j).local_size(); ++c)
          BlockHelper::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
  {
    using Number = typename Base<dim, VectorType>::value_type;
    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 < BlockHelper::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)
              BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
            ++c;
          }
        for (; c < BlockHelper::subblock(src_cpy, j).local_size(); ++c)
          BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
      }

    apply_add(dst, src_cpy);

    for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
      for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
        BlockHelper::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
  {
    using Number = typename Base<dim, VectorType>::value_type;
    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()
  {
    using Number = typename Base<dim, VectorType>::value_type;
    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) =
        Number(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
  {
    using Number = typename Base<dim, VectorType>::value_type;
    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()
  {
    using Number = typename Base<dim, VectorType>::value_type;
    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);

    this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
                          this,
                          diagonal_vector,
                          /*unused*/ diagonal_vector);
    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 Implementation
  {
    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;
    }
  } // namespace Implementation



  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(Implementation::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
  {
    using Number = typename Base<dim, VectorType>::value_type;
    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 VectorType &,
      const std::pair<unsigned int, unsigned int> &cell_range) const
  {
    using Number = typename Base<dim, VectorType>::value_type;
    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