simple.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2010 - 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------


#ifndef dealii_mesh_worker_simple_h
#define dealii_mesh_worker_simple_h

#include <deal.II/algorithms/any_data.h>

#include <deal.II/base/mg_level_object.h>
#include <deal.II/base/smartpointer.h>

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

#include <deal.II/meshworker/dof_info.h>
#include <deal.II/meshworker/functional.h>

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

/*
 * The header containing the classes MeshWorker::Assember::MatrixSimple,
 * MeshWorker::Assember::MGMatrixSimple, MeshWorker::Assember::ResidualSimple,
 * and MeshWorker::Assember::SystemSimple.
 */

DEAL_II_NAMESPACE_OPEN

namespace MeshWorker
{
  namespace Assembler
  {
    template <typename VectorType>
    class ResidualSimple
    {
    public:
      void
      initialize(AnyData &results);

      void
      initialize(
        const AffineConstraints<typename VectorType::value_type> &constraints);

      void
      initialize_local_blocks(const BlockIndices &);

      template <class DOFINFO>
      void
      initialize_info(DOFINFO &info, bool face) const;

      template <class DOFINFO>
      void
      assemble(const DOFINFO &info);

      template <class DOFINFO>
      void
      assemble(const DOFINFO &info1, const DOFINFO &info2);

    protected:
      AnyData residuals;

      SmartPointer<const AffineConstraints<typename VectorType::value_type>,
                   ResidualSimple<VectorType>>
        constraints;
    };


    template <typename MatrixType>
    class MatrixSimple
    {
    public:
      MatrixSimple(double threshold = 1.e-12);

      void
      initialize(MatrixType &m);

      void
      initialize(std::vector<MatrixType> &m);

      void
      initialize(
        const AffineConstraints<typename MatrixType::value_type> &constraints);

      template <class DOFINFO>
      void
      initialize_info(DOFINFO &info, bool face) const;

      template <class DOFINFO>
      void
      assemble(const DOFINFO &info);

      template <class DOFINFO>
      void
      assemble(const DOFINFO &info1, const DOFINFO &info2);

    protected:
      std::vector<SmartPointer<MatrixType, MatrixSimple<MatrixType>>> matrix;

      const double threshold;

    private:
      void
      assemble(const FullMatrix<double> &                  M,
               const unsigned int                          index,
               const std::vector<types::global_dof_index> &i1,
               const std::vector<types::global_dof_index> &i2);

      SmartPointer<const AffineConstraints<typename MatrixType::value_type>,
                   MatrixSimple<MatrixType>>
        constraints;
    };


    template <typename MatrixType>
    class MGMatrixSimple
    {
    public:
      MGMatrixSimple(double threshold = 1.e-12);

      void
      initialize(MGLevelObject<MatrixType> &m);

      void
      initialize(const MGConstrainedDoFs &mg_constrained_dofs);

      void
      initialize_fluxes(MGLevelObject<MatrixType> &flux_up,
                        MGLevelObject<MatrixType> &flux_down);

      void
      initialize_interfaces(MGLevelObject<MatrixType> &interface_in,
                            MGLevelObject<MatrixType> &interface_out);
      template <class DOFINFO>
      void
      initialize_info(DOFINFO &info, bool face) const;

      template <class DOFINFO>
      void
      assemble(const DOFINFO &info);

      template <class DOFINFO>
      void
      assemble(const DOFINFO &info1, const DOFINFO &info2);

    private:
      void
      assemble(MatrixType &                                G,
               const FullMatrix<double> &                  M,
               const std::vector<types::global_dof_index> &i1,
               const std::vector<types::global_dof_index> &i2);

      void
      assemble(MatrixType &                                G,
               const FullMatrix<double> &                  M,
               const std::vector<types::global_dof_index> &i1,
               const std::vector<types::global_dof_index> &i2,
               const unsigned int                          level);

      void
      assemble_up(MatrixType &                                G,
                  const FullMatrix<double> &                  M,
                  const std::vector<types::global_dof_index> &i1,
                  const std::vector<types::global_dof_index> &i2,
                  const unsigned int level = numbers::invalid_unsigned_int);
      void
      assemble_down(MatrixType &                                G,
                    const FullMatrix<double> &                  M,
                    const std::vector<types::global_dof_index> &i1,
                    const std::vector<types::global_dof_index> &i2,
                    const unsigned int level = numbers::invalid_unsigned_int);

      void
      assemble_in(MatrixType &                                G,
                  const FullMatrix<double> &                  M,
                  const std::vector<types::global_dof_index> &i1,
                  const std::vector<types::global_dof_index> &i2,
                  const unsigned int level = numbers::invalid_unsigned_int);

      void
      assemble_out(MatrixType &                                G,
                   const FullMatrix<double> &                  M,
                   const std::vector<types::global_dof_index> &i1,
                   const std::vector<types::global_dof_index> &i2,
                   const unsigned int level = numbers::invalid_unsigned_int);

      SmartPointer<MGLevelObject<MatrixType>, MGMatrixSimple<MatrixType>>
        matrix;

      SmartPointer<MGLevelObject<MatrixType>, MGMatrixSimple<MatrixType>>
        flux_up;

      SmartPointer<MGLevelObject<MatrixType>, MGMatrixSimple<MatrixType>>
        flux_down;

      SmartPointer<MGLevelObject<MatrixType>, MGMatrixSimple<MatrixType>>
        interface_in;

      SmartPointer<MGLevelObject<MatrixType>, MGMatrixSimple<MatrixType>>
        interface_out;
      SmartPointer<const MGConstrainedDoFs, MGMatrixSimple<MatrixType>>
        mg_constrained_dofs;

      const double threshold;
    };


    template <typename MatrixType, typename VectorType>
    class SystemSimple : private MatrixSimple<MatrixType>,
                         private ResidualSimple<VectorType>
    {
    public:
      SystemSimple(double threshold = 1.e-12);

      void
      initialize(MatrixType &m, VectorType &rhs);

      void
      initialize(
        const AffineConstraints<typename VectorType::value_type> &constraints);

      template <class DOFINFO>
      void
      initialize_info(DOFINFO &info, bool face) const;

      template <class DOFINFO>
      void
      assemble(const DOFINFO &info);

      template <class DOFINFO>
      void
      assemble(const DOFINFO &info1, const DOFINFO &info2);

    private:
      void
      assemble(const FullMatrix<double> &                  M,
               const Vector<double> &                      vector,
               const unsigned int                          index,
               const std::vector<types::global_dof_index> &indices);

      void
      assemble(const FullMatrix<double> &                  M,
               const Vector<double> &                      vector,
               const unsigned int                          index,
               const std::vector<types::global_dof_index> &i1,
               const std::vector<types::global_dof_index> &i2);
    };


    //----------------------------------------------------------------------//

    template <typename VectorType>
    inline void
    ResidualSimple<VectorType>::initialize(AnyData &results)
    {
      residuals = results;
    }

    template <typename VectorType>
    inline void
    ResidualSimple<VectorType>::initialize(
      const AffineConstraints<typename VectorType::value_type> &c)
    {
      constraints = &c;
    }


    template <typename MatrixType>
    inline void
    ResidualSimple<MatrixType>::initialize_local_blocks(const BlockIndices &)
    {}


    template <typename VectorType>
    template <class DOFINFO>
    inline void
    ResidualSimple<VectorType>::initialize_info(DOFINFO &info, bool) const
    {
      info.initialize_vectors(residuals.size());
    }


    template <typename VectorType>
    template <class DOFINFO>
    inline void
    ResidualSimple<VectorType>::assemble(const DOFINFO &info)
    {
      for (unsigned int k = 0; k < residuals.size(); ++k)
        {
          VectorType *v = residuals.entry<VectorType *>(k);
          for (unsigned int i = 0; i != info.vector(k).n_blocks(); ++i)
            {
              const std::vector<types::global_dof_index> &ldi =
                info.vector(k).n_blocks() == 1 ? info.indices :
                                                 info.indices_by_block[i];

              if (constraints != nullptr)
                constraints->distribute_local_to_global(info.vector(k).block(i),
                                                        ldi,
                                                        *v);
              else
                v->add(ldi, info.vector(k).block(i));
            }
        }
    }

    template <typename VectorType>
    template <class DOFINFO>
    inline void
    ResidualSimple<VectorType>::assemble(const DOFINFO &info1,
                                         const DOFINFO &info2)
    {
      assemble(info1);
      assemble(info2);
    }


    //----------------------------------------------------------------------//

    template <typename MatrixType>
    inline MatrixSimple<MatrixType>::MatrixSimple(double threshold)
      : threshold(threshold)
    {}


    template <typename MatrixType>
    inline void
    MatrixSimple<MatrixType>::initialize(MatrixType &m)
    {
      matrix.resize(1);
      matrix[0] = &m;
    }


    template <typename MatrixType>
    inline void
    MatrixSimple<MatrixType>::initialize(std::vector<MatrixType> &m)
    {
      matrix.resize(m.size());
      for (unsigned int i = 0; i < m.size(); ++i)
        matrix[i] = &m[i];
    }


    template <typename MatrixType>
    inline void
    MatrixSimple<MatrixType>::initialize(
      const AffineConstraints<typename MatrixType::value_type> &c)
    {
      constraints = &c;
    }


    template <typename MatrixType>
    template <class DOFINFO>
    inline void
    MatrixSimple<MatrixType>::initialize_info(DOFINFO &info, bool face) const
    {
      Assert(matrix.size() != 0, ExcNotInitialized());

      const unsigned int n = info.indices_by_block.size();

      if (n == 0)
        info.initialize_matrices(matrix.size(), face);
      else
        {
          info.initialize_matrices(matrix.size() * n * n, face);
          unsigned int k = 0;
          for (unsigned int m = 0; m < matrix.size(); ++m)
            for (unsigned int i = 0; i < n; ++i)
              for (unsigned int j = 0; j < n; ++j, ++k)
                {
                  info.matrix(k, false).row    = i;
                  info.matrix(k, false).column = j;
                  if (face)
                    {
                      info.matrix(k, true).row    = i;
                      info.matrix(k, true).column = j;
                    }
                }
        }
    }



    template <typename MatrixType>
    inline void
    MatrixSimple<MatrixType>::assemble(
      const FullMatrix<double> &                  M,
      const unsigned int                          index,
      const std::vector<types::global_dof_index> &i1,
      const std::vector<types::global_dof_index> &i2)
    {
      AssertDimension(M.m(), i1.size());
      AssertDimension(M.n(), i2.size());

      if (constraints == nullptr)
        {
          for (unsigned int j = 0; j < i1.size(); ++j)
            for (unsigned int k = 0; k < i2.size(); ++k)
              if (std::fabs(M(j, k)) >= threshold)
                matrix[index]->add(i1[j], i2[k], M(j, k));
        }
      else
        constraints->distribute_local_to_global(M, i1, i2, *matrix[index]);
    }


    template <typename MatrixType>
    template <class DOFINFO>
    inline void
    MatrixSimple<MatrixType>::assemble(const DOFINFO &info)
    {
      Assert(!info.level_cell, ExcMessage("Cell may not access level dofs"));
      const unsigned int n = info.indices_by_block.size();

      if (n == 0)
        for (unsigned int m = 0; m < matrix.size(); ++m)
          assemble(info.matrix(m, false).matrix, m, info.indices, info.indices);
      else
        {
          for (unsigned int m = 0; m < matrix.size(); ++m)
            for (unsigned int k = 0; k < n * n; ++k)
              {
                assemble(
                  info.matrix(k + m * n * n, false).matrix,
                  m,
                  info.indices_by_block[info.matrix(k + m * n * n, false).row],
                  info.indices_by_block[info.matrix(k + m * n * n, false)
                                          .column]);
              }
        }
    }


    template <typename MatrixType>
    template <class DOFINFO>
    inline void
    MatrixSimple<MatrixType>::assemble(const DOFINFO &info1,
                                       const DOFINFO &info2)
    {
      Assert(!info1.level_cell, ExcMessage("Cell may not access level dofs"));
      Assert(!info2.level_cell, ExcMessage("Cell may not access level dofs"));
      AssertDimension(info1.indices_by_block.size(),
                      info2.indices_by_block.size());

      const unsigned int n = info1.indices_by_block.size();

      if (n == 0)
        {
          for (unsigned int m = 0; m < matrix.size(); ++m)
            {
              assemble(info1.matrix(m, false).matrix,
                       m,
                       info1.indices,
                       info1.indices);
              assemble(info1.matrix(m, true).matrix,
                       m,
                       info1.indices,
                       info2.indices);
              assemble(info2.matrix(m, false).matrix,
                       m,
                       info2.indices,
                       info2.indices);
              assemble(info2.matrix(m, true).matrix,
                       m,
                       info2.indices,
                       info1.indices);
            }
        }
      else
        {
          for (unsigned int m = 0; m < matrix.size(); ++m)
            for (unsigned int k = 0; k < n * n; ++k)
              {
                const unsigned int row = info1.matrix(k + m * n * n, false).row;
                const unsigned int column =
                  info1.matrix(k + m * n * n, false).column;

                assemble(info1.matrix(k + m * n * n, false).matrix,
                         m,
                         info1.indices_by_block[row],
                         info1.indices_by_block[column]);
                assemble(info1.matrix(k + m * n * n, true).matrix,
                         m,
                         info1.indices_by_block[row],
                         info2.indices_by_block[column]);
                assemble(info2.matrix(k + m * n * n, false).matrix,
                         m,
                         info2.indices_by_block[row],
                         info2.indices_by_block[column]);
                assemble(info2.matrix(k + m * n * n, true).matrix,
                         m,
                         info2.indices_by_block[row],
                         info1.indices_by_block[column]);
              }
        }
    }


    //----------------------------------------------------------------------//

    template <typename MatrixType>
    inline MGMatrixSimple<MatrixType>::MGMatrixSimple(double threshold)
      : threshold(threshold)
    {}


    template <typename MatrixType>
    inline void
    MGMatrixSimple<MatrixType>::initialize(MGLevelObject<MatrixType> &m)
    {
      matrix = &m;
    }

    template <typename MatrixType>
    inline void
    MGMatrixSimple<MatrixType>::initialize(const MGConstrainedDoFs &c)
    {
      mg_constrained_dofs = &c;
    }


    template <typename MatrixType>
    inline void
    MGMatrixSimple<MatrixType>::initialize_fluxes(
      MGLevelObject<MatrixType> &up,
      MGLevelObject<MatrixType> &down)
    {
      flux_up   = &up;
      flux_down = &down;
    }


    template <typename MatrixType>
    inline void
    MGMatrixSimple<MatrixType>::initialize_interfaces(
      MGLevelObject<MatrixType> &in,
      MGLevelObject<MatrixType> &out)
    {
      interface_in  = &in;
      interface_out = &out;
    }


    template <typename MatrixType>
    template <class DOFINFO>
    inline void
    MGMatrixSimple<MatrixType>::initialize_info(DOFINFO &info, bool face) const
    {
      const unsigned int n = info.indices_by_block.size();

      if (n == 0)
        info.initialize_matrices(1, face);
      else
        {
          info.initialize_matrices(n * n, face);
          unsigned int k = 0;
          for (unsigned int i = 0; i < n; ++i)
            for (unsigned int j = 0; j < n; ++j, ++k)
              {
                info.matrix(k, false).row    = i;
                info.matrix(k, false).column = j;
                if (face)
                  {
                    info.matrix(k, true).row    = i;
                    info.matrix(k, true).column = j;
                  }
              }
        }
    }


    template <typename MatrixType>
    inline void
    MGMatrixSimple<MatrixType>::assemble(
      MatrixType &                                G,
      const FullMatrix<double> &                  M,
      const std::vector<types::global_dof_index> &i1,
      const std::vector<types::global_dof_index> &i2)
    {
      AssertDimension(M.m(), i1.size());
      AssertDimension(M.n(), i2.size());
      Assert(mg_constrained_dofs == 0, ExcInternalError());
      // TODO: Possibly remove this function all together

      for (unsigned int j = 0; j < i1.size(); ++j)
        for (unsigned int k = 0; k < i2.size(); ++k)
          if (std::fabs(M(j, k)) >= threshold)
            G.add(i1[j], i2[k], M(j, k));
    }


    template <typename MatrixType>
    inline void
    MGMatrixSimple<MatrixType>::assemble(
      MatrixType &                                G,
      const FullMatrix<double> &                  M,
      const std::vector<types::global_dof_index> &i1,
      const std::vector<types::global_dof_index> &i2,
      const unsigned int                          level)
    {
      AssertDimension(M.m(), i1.size());
      AssertDimension(M.n(), i2.size());

      if (mg_constrained_dofs == nullptr)
        {
          for (unsigned int j = 0; j < i1.size(); ++j)
            for (unsigned int k = 0; k < i2.size(); ++k)
              if (std::fabs(M(j, k)) >= threshold)
                G.add(i1[j], i2[k], M(j, k));
        }
      else
        {
          for (unsigned int j = 0; j < i1.size(); ++j)
            for (unsigned int k = 0; k < i2.size(); ++k)
              {
                // Only enter the local values into the global matrix,
                //  if the value is larger than the threshold
                if (std::fabs(M(j, k)) < threshold)
                  continue;

                // Do not enter, if either the row or the column
                // corresponds to an index on the refinement edge. The
                // level problems are solved with homogeneous
                // Dirichlet boundary conditions, therefore we
                // eliminate these rows and columns. The corresponding
                // matrix entries are entered by assemble_in() and
                // assemble_out().
                if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) ||
                    mg_constrained_dofs->at_refinement_edge(level, i2[k]))
                  continue;

                // At the boundary, only enter the term on the
                // diagonal, but not the coupling terms
                if ((mg_constrained_dofs->is_boundary_index(level, i1[j]) ||
                     mg_constrained_dofs->is_boundary_index(level, i2[k])) &&
                    (i1[j] != i2[k]))
                  continue;

                G.add(i1[j], i2[k], M(j, k));
              }
        }
    }


    template <typename MatrixType>
    inline void
    MGMatrixSimple<MatrixType>::assemble_up(
      MatrixType &                                G,
      const FullMatrix<double> &                  M,
      const std::vector<types::global_dof_index> &i1,
      const std::vector<types::global_dof_index> &i2,
      const unsigned int                          level)
    {
      AssertDimension(M.n(), i1.size());
      AssertDimension(M.m(), i2.size());

      if (mg_constrained_dofs == nullptr)
        {
          for (unsigned int j = 0; j < i1.size(); ++j)
            for (unsigned int k = 0; k < i2.size(); ++k)
              if (std::fabs(M(k, j)) >= threshold)
                G.add(i1[j], i2[k], M(k, j));
        }
      else
        {
          for (unsigned int j = 0; j < i1.size(); ++j)
            for (unsigned int k = 0; k < i2.size(); ++k)
              if (std::fabs(M(k, j)) >= threshold)
                if (!mg_constrained_dofs->at_refinement_edge(level, i2[k]))
                  G.add(i1[j], i2[k], M(k, j));
        }
    }

    template <typename MatrixType>
    inline void
    MGMatrixSimple<MatrixType>::assemble_down(
      MatrixType &                                G,
      const FullMatrix<double> &                  M,
      const std::vector<types::global_dof_index> &i1,
      const std::vector<types::global_dof_index> &i2,
      const unsigned int                          level)
    {
      AssertDimension(M.m(), i1.size());
      AssertDimension(M.n(), i2.size());

      if (mg_constrained_dofs == nullptr)
        {
          for (unsigned int j = 0; j < i1.size(); ++j)
            for (unsigned int k = 0; k < i2.size(); ++k)
              if (std::fabs(M(j, k)) >= threshold)
                G.add(i1[j], i2[k], M(j, k));
        }
      else
        {
          for (unsigned int j = 0; j < i1.size(); ++j)
            for (unsigned int k = 0; k < i2.size(); ++k)
              if (std::fabs(M(j, k)) >= threshold)
                if (!mg_constrained_dofs->at_refinement_edge(level, i2[k]))
                  G.add(i1[j], i2[k], M(j, k));
        }
    }

    template <typename MatrixType>
    inline void
    MGMatrixSimple<MatrixType>::assemble_in(
      MatrixType &                                G,
      const FullMatrix<double> &                  M,
      const std::vector<types::global_dof_index> &i1,
      const std::vector<types::global_dof_index> &i2,
      const unsigned int                          level)
    {
      AssertDimension(M.m(), i1.size());
      AssertDimension(M.n(), i2.size());
      Assert(mg_constrained_dofs != nullptr, ExcInternalError());

      for (unsigned int j = 0; j < i1.size(); ++j)
        for (unsigned int k = 0; k < i2.size(); ++k)
          if (std::fabs(M(j, k)) >= threshold)
            // Enter values into matrix only if j corresponds to a
            // degree of freedom on the refinement edge, k does
            // not, and both are not on the boundary. This is part
            // the difference between the complete matrix with no
            // boundary condition at the refinement edge and and
            // the matrix assembled above by assemble().

            // Thus the logic is: enter the row if it is
            // constrained by hanging node constraints (actually,
            // the whole refinement edge), but not if it is
            // constrained by a boundary constraint.
            if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
                !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
              {
                if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
                     !mg_constrained_dofs->is_boundary_index(level, i2[k])) ||
                    (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
                     mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
                     i1[j] == i2[k]))
                  G.add(i1[j], i2[k], M(j, k));
              }
    }


    template <typename MatrixType>
    inline void
    MGMatrixSimple<MatrixType>::assemble_out(
      MatrixType &                                G,
      const FullMatrix<double> &                  M,
      const std::vector<types::global_dof_index> &i1,
      const std::vector<types::global_dof_index> &i2,
      const unsigned int                          level)
    {
      AssertDimension(M.n(), i1.size());
      AssertDimension(M.m(), i2.size());
      Assert(mg_constrained_dofs != nullptr, ExcInternalError());

      for (unsigned int j = 0; j < i1.size(); ++j)
        for (unsigned int k = 0; k < i2.size(); ++k)
          if (std::fabs(M(k, j)) >= threshold)
            if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
                !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
              {
                if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
                     !mg_constrained_dofs->is_boundary_index(level, i2[k])) ||
                    (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
                     mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
                     i1[j] == i2[k]))
                  G.add(i1[j], i2[k], M(k, j));
              }
    }


    template <typename MatrixType>
    template <class DOFINFO>
    inline void
    MGMatrixSimple<MatrixType>::assemble(const DOFINFO &info)
    {
      Assert(info.level_cell, ExcMessage("Cell must access level dofs"));
      const unsigned int level = info.cell->level();

      if (info.indices_by_block.size() == 0)
        {
          assemble((*matrix)[level],
                   info.matrix(0, false).matrix,
                   info.indices,
                   info.indices,
                   level);
          if (mg_constrained_dofs != nullptr)
            {
              assemble_in((*interface_in)[level],
                          info.matrix(0, false).matrix,
                          info.indices,
                          info.indices,
                          level);
              assemble_out((*interface_out)[level],
                           info.matrix(0, false).matrix,
                           info.indices,
                           info.indices,
                           level);
            }
        }
      else
        for (unsigned int k = 0; k < info.n_matrices(); ++k)
          {
            const unsigned int row    = info.matrix(k, false).row;
            const unsigned int column = info.matrix(k, false).column;

            assemble((*matrix)[level],
                     info.matrix(k, false).matrix,
                     info.indices_by_block[row],
                     info.indices_by_block[column],
                     level);

            if (mg_constrained_dofs != nullptr)
              {
                assemble_in((*interface_in)[level],
                            info.matrix(k, false).matrix,
                            info.indices_by_block[row],
                            info.indices_by_block[column],
                            level);
                assemble_out((*interface_out)[level],
                             info.matrix(k, false).matrix,
                             info.indices_by_block[column],
                             info.indices_by_block[row],
                             level);
              }
          }
    }


    template <typename MatrixType>
    template <class DOFINFO>
    inline void
    MGMatrixSimple<MatrixType>::assemble(const DOFINFO &info1,
                                         const DOFINFO &info2)
    {
      Assert(info1.level_cell, ExcMessage("Cell must access level dofs"));
      Assert(info2.level_cell, ExcMessage("Cell must access level dofs"));
      const unsigned int level1 = info1.cell->level();
      const unsigned int level2 = info2.cell->level();

      if (info1.indices_by_block.size() == 0)
        {
          if (level1 == level2)
            {
              assemble((*matrix)[level1],
                       info1.matrix(0, false).matrix,
                       info1.indices,
                       info1.indices,
                       level1);
              assemble((*matrix)[level1],
                       info1.matrix(0, true).matrix,
                       info1.indices,
                       info2.indices,
                       level1);
              assemble((*matrix)[level1],
                       info2.matrix(0, false).matrix,
                       info2.indices,
                       info2.indices,
                       level1);
              assemble((*matrix)[level1],
                       info2.matrix(0, true).matrix,
                       info2.indices,
                       info1.indices,
                       level1);
            }
          else
            {
              Assert(level1 > level2, ExcInternalError());
              // Do not add info2.M1,
              // which is done by
              // the coarser cell
              assemble((*matrix)[level1],
                       info1.matrix(0, false).matrix,
                       info1.indices,
                       info1.indices,
                       level1);
              if (level1 > 0)
                {
                  assemble_up((*flux_up)[level1],
                              info1.matrix(0, true).matrix,
                              info2.indices,
                              info1.indices,
                              level1);
                  assemble_down((*flux_down)[level1],
                                info2.matrix(0, true).matrix,
                                info2.indices,
                                info1.indices,
                                level1);
                }
            }
        }
      else
        for (unsigned int k = 0; k < info1.n_matrices(); ++k)
          {
            const unsigned int row    = info1.matrix(k, false).row;
            const unsigned int column = info1.matrix(k, false).column;

            if (level1 == level2)
              {
                assemble((*matrix)[level1],
                         info1.matrix(k, false).matrix,
                         info1.indices_by_block[row],
                         info1.indices_by_block[column],
                         level1);
                assemble((*matrix)[level1],
                         info1.matrix(k, true).matrix,
                         info1.indices_by_block[row],
                         info2.indices_by_block[column],
                         level1);
                assemble((*matrix)[level1],
                         info2.matrix(k, false).matrix,
                         info2.indices_by_block[row],
                         info2.indices_by_block[column],
                         level1);
                assemble((*matrix)[level1],
                         info2.matrix(k, true).matrix,
                         info2.indices_by_block[row],
                         info1.indices_by_block[column],
                         level1);
              }
            else
              {
                Assert(level1 > level2, ExcInternalError());
                // Do not add info2.M1,
                // which is done by
                // the coarser cell
                assemble((*matrix)[level1],
                         info1.matrix(k, false).matrix,
                         info1.indices_by_block[row],
                         info1.indices_by_block[column],
                         level1);
                if (level1 > 0)
                  {
                    assemble_up((*flux_up)[level1],
                                info1.matrix(k, true).matrix,
                                info2.indices_by_block[column],
                                info1.indices_by_block[row],
                                level1);
                    assemble_down((*flux_down)[level1],
                                  info2.matrix(k, true).matrix,
                                  info2.indices_by_block[row],
                                  info1.indices_by_block[column],
                                  level1);
                  }
              }
          }
    }

    //----------------------------------------------------------------------//

    template <typename MatrixType, typename VectorType>
    SystemSimple<MatrixType, VectorType>::SystemSimple(double t)
      : MatrixSimple<MatrixType>(t)
    {}


    template <typename MatrixType, typename VectorType>
    inline void
    SystemSimple<MatrixType, VectorType>::initialize(MatrixType &m,
                                                     VectorType &rhs)
    {
      AnyData     data;
      VectorType *p = &rhs;
      data.add(p, "right hand side");

      MatrixSimple<MatrixType>::initialize(m);
      ResidualSimple<VectorType>::initialize(data);
    }

    template <typename MatrixType, typename VectorType>
    inline void
    SystemSimple<MatrixType, VectorType>::initialize(
      const AffineConstraints<typename VectorType::value_type> &c)
    {
      ResidualSimple<VectorType>::initialize(c);
    }


    template <typename MatrixType, typename VectorType>
    template <class DOFINFO>
    inline void
    SystemSimple<MatrixType, VectorType>::initialize_info(DOFINFO &info,
                                                          bool     face) const
    {
      MatrixSimple<MatrixType>::initialize_info(info, face);
      ResidualSimple<VectorType>::initialize_info(info, face);
    }

    template <typename MatrixType, typename VectorType>
    inline void
    SystemSimple<MatrixType, VectorType>::assemble(
      const FullMatrix<double> &                  M,
      const Vector<double> &                      vector,
      const unsigned int                          index,
      const std::vector<types::global_dof_index> &indices)
    {
      AssertDimension(M.m(), indices.size());
      AssertDimension(M.n(), indices.size());

      AnyData     residuals = ResidualSimple<VectorType>::residuals;
      VectorType *v         = residuals.entry<VectorType *>(index);

      if (ResidualSimple<VectorType>::constraints == nullptr)
        {
          for (unsigned int i = 0; i < indices.size(); ++i)
            (*v)(indices[i]) += vector(i);

          for (unsigned int j = 0; j < indices.size(); ++j)
            for (unsigned int k = 0; k < indices.size(); ++k)
              if (std::fabs(M(j, k)) >= MatrixSimple<MatrixType>::threshold)
                MatrixSimple<MatrixType>::matrix[index]->add(indices[j],
                                                             indices[k],
                                                             M(j, k));
        }
      else
        {
          ResidualSimple<VectorType>::constraints->distribute_local_to_global(
            M,
            vector,
            indices,
            *MatrixSimple<MatrixType>::matrix[index],
            *v,
            true);
        }
    }

    template <typename MatrixType, typename VectorType>
    inline void
    SystemSimple<MatrixType, VectorType>::assemble(
      const FullMatrix<double> &                  M,
      const Vector<double> &                      vector,
      const unsigned int                          index,
      const std::vector<types::global_dof_index> &i1,
      const std::vector<types::global_dof_index> &i2)
    {
      AssertDimension(M.m(), i1.size());
      AssertDimension(M.n(), i2.size());

      AnyData     residuals = ResidualSimple<VectorType>::residuals;
      VectorType *v         = residuals.entry<VectorType *>(index);

      if (ResidualSimple<VectorType>::constraints == nullptr)
        {
          for (unsigned int j = 0; j < i1.size(); ++j)
            for (unsigned int k = 0; k < i2.size(); ++k)
              if (std::fabs(M(j, k)) >= MatrixSimple<MatrixType>::threshold)
                MatrixSimple<MatrixType>::matrix[index]->add(i1[j],
                                                             i2[k],
                                                             M(j, k));
        }
      else
        {
          ResidualSimple<VectorType>::constraints->distribute_local_to_global(
            vector, i1, i2, *v, M, false);
          ResidualSimple<VectorType>::constraints->distribute_local_to_global(
            M, i1, i2, *MatrixSimple<MatrixType>::matrix[index]);
        }
    }


    template <typename MatrixType, typename VectorType>
    template <class DOFINFO>
    inline void
    SystemSimple<MatrixType, VectorType>::assemble(const DOFINFO &info)
    {
      AssertDimension(MatrixSimple<MatrixType>::matrix.size(),
                      ResidualSimple<VectorType>::residuals.size());
      Assert(!info.level_cell, ExcMessage("Cell may not access level dofs"));
      const unsigned int n = info.indices_by_block.size();

      if (n == 0)
        {
          for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
               ++m)
            assemble(info.matrix(m, false).matrix,
                     info.vector(m).block(0),
                     m,
                     info.indices);
        }
      else
        {
          for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
               ++m)
            for (unsigned int k = 0; k < n * n; ++k)
              {
                const unsigned int row = info.matrix(k + m * n * n, false).row;
                const unsigned int column =
                  info.matrix(k + m * n * n, false).column;

                if (row == column)
                  assemble(info.matrix(k + m * n * n, false).matrix,
                           info.vector(m).block(row),
                           m,
                           info.indices_by_block[row]);
                else
                  assemble(info.matrix(k + m * n * n, false).matrix,
                           info.vector(m).block(row),
                           m,
                           info.indices_by_block[row],
                           info.indices_by_block[column]);
              }
        }
    }


    template <typename MatrixType, typename VectorType>
    template <class DOFINFO>
    inline void
    SystemSimple<MatrixType, VectorType>::assemble(const DOFINFO &info1,
                                                   const DOFINFO &info2)
    {
      Assert(!info1.level_cell, ExcMessage("Cell may not access level dofs"));
      Assert(!info2.level_cell, ExcMessage("Cell may not access level dofs"));
      AssertDimension(info1.indices_by_block.size(),
                      info2.indices_by_block.size());

      const unsigned int n = info1.indices_by_block.size();

      if (n == 0)
        {
          for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
               ++m)
            {
              assemble(info1.matrix(m, false).matrix,
                       info1.vector(m).block(0),
                       m,
                       info1.indices);
              assemble(info1.matrix(m, true).matrix,
                       info1.vector(m).block(0),
                       m,
                       info1.indices,
                       info2.indices);
              assemble(info2.matrix(m, false).matrix,
                       info2.vector(m).block(0),
                       m,
                       info2.indices);
              assemble(info2.matrix(m, true).matrix,
                       info2.vector(m).block(0),
                       m,
                       info2.indices,
                       info1.indices);
            }
        }
      else
        {
          for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
               ++m)
            for (unsigned int k = 0; k < n * n; ++k)
              {
                const unsigned int row = info1.matrix(k + m * n * n, false).row;
                const unsigned int column =
                  info1.matrix(k + m * n * n, false).column;

                if (row == column)
                  {
                    assemble(info1.matrix(k + m * n * n, false).matrix,
                             info1.vector(m).block(row),
                             m,
                             info1.indices_by_block[row]);
                    assemble(info2.matrix(k + m * n * n, false).matrix,
                             info2.vector(m).block(row),
                             m,
                             info2.indices_by_block[row]);
                  }
                else
                  {
                    assemble(info1.matrix(k + m * n * n, false).matrix,
                             info1.vector(m).block(row),
                             m,
                             info1.indices_by_block[row],
                             info1.indices_by_block[column]);
                    assemble(info2.matrix(k + m * n * n, false).matrix,
                             info2.vector(m).block(row),
                             m,
                             info2.indices_by_block[row],
                             info2.indices_by_block[column]);
                  }
                assemble(info1.matrix(k + m * n * n, true).matrix,
                         info1.vector(m).block(row),
                         m,
                         info1.indices_by_block[row],
                         info2.indices_by_block[column]);
                assemble(info2.matrix(k + m * n * n, true).matrix,
                         info2.vector(m).block(row),
                         m,
                         info2.indices_by_block[row],
                         info1.indices_by_block[column]);
              }
        }
    }
  } // namespace Assembler
} // namespace MeshWorker

DEAL_II_NAMESPACE_CLOSE

#endif