assembler.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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------


#ifndef dealii_mesh_worker_assembler_h
#define dealii_mesh_worker_assembler_h

#include <deal.II/base/smartpointer.h>
#include <deal.II/base/mg_level_object.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/meshworker/dof_info.h>
#include <deal.II/meshworker/functional.h>
#include <deal.II/meshworker/simple.h>
#include <deal.II/multigrid/mg_constrained_dofs.h>
#include <deal.II/algorithms/any_data.h>


DEAL_II_NAMESPACE_OPEN

namespace MeshWorker
{
  namespace Assembler
  {
    template <typename VectorType>
    class ResidualLocalBlocksToGlobalBlocks
    {
    public:
      void initialize(const BlockInfo *block_info,
                      AnyData &residuals);
      void initialize(const ConstraintMatrix &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(VectorType                                 &global,
                    const BlockVector<double>                  &local,
                    const std::vector<types::global_dof_index> &dof);

      AnyData residuals;

      SmartPointer<const BlockInfo,
                   ResidualLocalBlocksToGlobalBlocks<VectorType> > block_info;
      SmartPointer<const ConstraintMatrix,
                   ResidualLocalBlocksToGlobalBlocks<VectorType> > constraints;
    };


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

      void initialize(const BlockInfo *block_info,
                      MatrixBlockVector<MatrixType> &matrices);

      void initialize(const ConstraintMatrix &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(MatrixBlock<MatrixType>                    &global,
                    const FullMatrix<number>                   &local,
                    const unsigned int                          block_row,
                    const unsigned int                          block_col,
                    const std::vector<types::global_dof_index> &dof1,
                    const std::vector<types::global_dof_index> &dof2);

      SmartPointer<MatrixBlockVector<MatrixType>,
                   MatrixLocalBlocksToGlobalBlocks<MatrixType, number> > matrices;

      SmartPointer<const BlockInfo,
                   MatrixLocalBlocksToGlobalBlocks<MatrixType, number> > block_info;
      SmartPointer<const ConstraintMatrix,
                   MatrixLocalBlocksToGlobalBlocks<MatrixType,number> > constraints;

      const double threshold;

    };

    template <typename MatrixType, typename number = double>
    class MGMatrixLocalBlocksToGlobalBlocks
    {
    public:
      typedef MGMatrixBlockVector<MatrixType> MatrixPtrVector;
      typedef SmartPointer<MatrixPtrVector, MGMatrixLocalBlocksToGlobalBlocks<MatrixType,number> >
      MatrixPtrVectorPtr;

      MGMatrixLocalBlocksToGlobalBlocks(double threshold = 1.e-12);

      void initialize(const BlockInfo *block_info,
                      MatrixPtrVector &matrices);

      void initialize(const MGConstrainedDoFs &mg_constrained_dofs);

      void initialize_edge_flux(MatrixPtrVector &up, MatrixPtrVector &down);

      void initialize_interfaces (MatrixPtrVector &interface_in, MatrixPtrVector &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                                 &global,
                    const FullMatrix<number>                   &local,
                    const unsigned int                          block_row,
                    const unsigned int                          block_col,
                    const std::vector<types::global_dof_index> &dof1,
                    const std::vector<types::global_dof_index> &dof2,
                    const unsigned int                          level1,
                    const unsigned int                          level2,
                    bool                                        transpose = false);

      void assemble_fluxes(MatrixType                                 &global,
                           const FullMatrix<number>                   &local,
                           const unsigned int                          block_row,
                           const unsigned int                          block_col,
                           const std::vector<types::global_dof_index> &dof1,
                           const std::vector<types::global_dof_index> &dof2,
                           const unsigned int                          level1,
                           const unsigned int                          level2);

      void assemble_up(MatrixType                                 &global,
                       const FullMatrix<number>                   &local,
                       const unsigned int                          block_row,
                       const unsigned int                          block_col,
                       const std::vector<types::global_dof_index> &dof1,
                       const std::vector<types::global_dof_index> &dof2,
                       const unsigned int                          level1,
                       const unsigned int                          level2);

      void assemble_down(MatrixType                                 &global,
                         const FullMatrix<number>                   &local,
                         const unsigned int                          block_row,
                         const unsigned int                          block_col,
                         const std::vector<types::global_dof_index> &dof1,
                         const std::vector<types::global_dof_index> &dof2,
                         const unsigned int                          level1,
                         const unsigned int                          level2);

      void assemble_in(MatrixType                                 &global,
                       const FullMatrix<number>                   &local,
                       const unsigned int                          block_row,
                       const unsigned int                          block_col,
                       const std::vector<types::global_dof_index> &dof1,
                       const std::vector<types::global_dof_index> &dof2,
                       const unsigned int                          level1,
                       const unsigned int                          level2);

      void assemble_out(MatrixType                                 &global,
                        const FullMatrix<number>                   &local,
                        const unsigned int                          block_row,
                        const unsigned int                          block_col,
                        const std::vector<types::global_dof_index> &dof1,
                        const std::vector<types::global_dof_index> &dof2,
                        const unsigned int                          level1,
                        const unsigned int                          level2);

      MatrixPtrVectorPtr matrices;

      MatrixPtrVectorPtr flux_down;

      MatrixPtrVectorPtr flux_up;

      MatrixPtrVectorPtr interface_out;

      MatrixPtrVectorPtr interface_in;

      SmartPointer<const BlockInfo, MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number> > block_info;

      SmartPointer<const MGConstrainedDoFs,MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number> > mg_constrained_dofs;


      const double threshold;

    };

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

    template <typename VectorType>
    inline void
    ResidualLocalBlocksToGlobalBlocks<VectorType>::initialize(const BlockInfo *b,
                                                              AnyData         &m)
    {
      block_info = b;
      residuals = m;
    }

    template <typename VectorType>
    inline void
    ResidualLocalBlocksToGlobalBlocks<VectorType>::initialize
    (const ConstraintMatrix &c)
    {
      constraints = &c;
    }


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

    template <typename VectorType>
    inline void
    ResidualLocalBlocksToGlobalBlocks<VectorType>::assemble
    (VectorType                                 &global,
     const BlockVector<double>                  &local,
     const std::vector<types::global_dof_index> &dof)
    {
      if (constraints == 0)
        {
          for (unsigned int b=0; b<local.n_blocks(); ++b)
            for (unsigned int j=0; j<local.block(b).size(); ++j)
              {
                // The coordinates of
                // the current entry in
                // DoFHandler
                // numbering, which
                // differs from the
                // block-wise local
                // numbering we use in
                // our local vectors
                const unsigned int jcell = this->block_info->local().local_to_global(b, j);
                global(dof[jcell]) += local.block(b)(j);
              }
        }
      else
        constraints->distribute_local_to_global(local, dof, global);
    }


    template <typename VectorType>
    template <class DOFINFO>
    inline void
    ResidualLocalBlocksToGlobalBlocks<VectorType>::assemble
    (const DOFINFO &info)
    {
      for (unsigned int i=0; i<residuals.size(); ++i)
        assemble(*(residuals.entry<VectorType>(i)), info.vector(i), info.indices);
    }


    template <typename VectorType>
    template <class DOFINFO>
    inline void
    ResidualLocalBlocksToGlobalBlocks<VectorType>::assemble
    (const DOFINFO &info1,
     const DOFINFO &info2)
    {
      for (unsigned int i=0; i<residuals.size(); ++i)
        {
          assemble(*(residuals.entry<VectorType>(i)), info1.vector(i), info1.indices);
          assemble(*(residuals.entry<VectorType>(i)), info2.vector(i), info2.indices);
        }
    }


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

    template <typename MatrixType, typename number>
    inline
    MatrixLocalBlocksToGlobalBlocks<MatrixType, number>::MatrixLocalBlocksToGlobalBlocks
    (double threshold)
      :
      threshold(threshold)
    {}


    template <typename MatrixType, typename number>
    inline void
    MatrixLocalBlocksToGlobalBlocks<MatrixType, number>::initialize
    (const BlockInfo               *b,
     MatrixBlockVector<MatrixType> &m)
    {
      block_info = b;
      matrices = &m;
    }



    template <typename MatrixType, typename number>
    inline void
    MatrixLocalBlocksToGlobalBlocks<MatrixType, number>::initialize
    (const ConstraintMatrix &c)
    {
      constraints = &c;
    }



    template <typename MatrixType,typename number>
    template <class DOFINFO>
    inline void
    MatrixLocalBlocksToGlobalBlocks<MatrixType, number>::initialize_info
    (DOFINFO &info,
     bool face) const
    {
      info.initialize_matrices(*matrices, face);
    }



    template <typename MatrixType, typename number>
    inline void
    MatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble
    (MatrixBlock<MatrixType>                    &global,
     const FullMatrix<number>                   &local,
     const unsigned int                          block_row,
     const unsigned int                          block_col,
     const std::vector<types::global_dof_index> &dof1,
     const std::vector<types::global_dof_index> &dof2)
    {
      if (constraints == nullptr)
        {
          for (unsigned int j=0; j<local.n_rows(); ++j)
            for (unsigned int k=0; k<local.n_cols(); ++k)
              if (std::fabs(local(j,k)) >= threshold)
                {
                  // The coordinates of
                  // the current entry in
                  // DoFHandler
                  // numbering, which
                  // differs from the
                  // block-wise local
                  // numbering we use in
                  // our local matrices
                  const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
                  const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);

                  global.add(dof1[jcell], dof2[kcell], local(j,k));
                }
        }
      else
        {
          const BlockIndices &bi = this->block_info->local();
          std::vector<types::global_dof_index> sliced_row_indices (bi.block_size(block_row));
          for (unsigned int i=0; i<sliced_row_indices.size(); ++i)
            sliced_row_indices[i] = dof1[bi.block_start(block_row)+i];

          std::vector<types::global_dof_index> sliced_col_indices (bi.block_size(block_col));
          for (unsigned int i=0; i<sliced_col_indices.size(); ++i)
            sliced_col_indices[i] = dof2[bi.block_start(block_col)+i];

          constraints->distribute_local_to_global(local,
                                                  sliced_row_indices, sliced_col_indices, global);
        }
    }


    template <typename MatrixType, typename number>
    template <class DOFINFO>
    inline void
    MatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble (const DOFINFO &info)
    {
      for (unsigned int i=0; i<matrices->size(); ++i)
        {
          // Row and column index of
          // the block we are dealing with
          const types::global_dof_index row = matrices->block(i).row;
          const types::global_dof_index col = matrices->block(i).column;

          assemble(matrices->block(i), info.matrix(i,false).matrix, row, col, info.indices, info.indices);
        }
    }


    template <typename MatrixType, typename number>
    template <class DOFINFO>
    inline void
    MatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble (const DOFINFO &info1,
        const DOFINFO &info2)
    {
      for (unsigned int i=0; i<matrices->size(); ++i)
        {
          // Row and column index of
          // the block we are dealing with
          const types::global_dof_index row = matrices->block(i).row;
          const types::global_dof_index col = matrices->block(i).column;

          assemble(matrices->block(i), info1.matrix(i,false).matrix, row, col, info1.indices, info1.indices);
          assemble(matrices->block(i), info1.matrix(i,true).matrix, row, col, info1.indices, info2.indices);
          assemble(matrices->block(i), info2.matrix(i,false).matrix, row, col, info2.indices, info2.indices);
          assemble(matrices->block(i), info2.matrix(i,true).matrix, row, col, info2.indices, info1.indices);
        }
    }


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

    template <typename MatrixType, typename number>
    inline
    MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::MGMatrixLocalBlocksToGlobalBlocks
    (double threshold)
      :
      threshold(threshold)
    {}


    template <typename MatrixType, typename number>
    inline void
    MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::initialize
    (const BlockInfo *b,
     MatrixPtrVector &m)
    {
      block_info = b;
      AssertDimension(block_info->local().size(), block_info->global().size());
      matrices = &m;
    }


    template <typename MatrixType, typename number>
    inline void
    MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::initialize
    (const MGConstrainedDoFs &mg_c)
    {
      mg_constrained_dofs = &mg_c;
    }


    template <typename MatrixType,typename number>
    template <class DOFINFO>
    inline void
    MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::initialize_info
    (DOFINFO &info,
     bool face) const
    {
      info.initialize_matrices(*matrices, face);
    }



    template <typename MatrixType, typename number>
    inline void
    MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::initialize_edge_flux
    (MatrixPtrVector &up,
     MatrixPtrVector &down)
    {
      flux_up = up;
      flux_down = down;
    }


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


    template <typename MatrixType, typename number>
    inline void
    MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble
    (MatrixType                                 &global,
     const FullMatrix<number>                   &local,
     const unsigned int                          block_row,
     const unsigned int                          block_col,
     const std::vector<types::global_dof_index> &dof1,
     const std::vector<types::global_dof_index> &dof2,
     const unsigned int                          level1,
     const unsigned int                          level2,
     bool                                        transpose)
    {
      for (unsigned int j=0; j<local.n_rows(); ++j)
        for (unsigned int k=0; k<local.n_cols(); ++k)
          if (std::fabs(local(j,k)) >= threshold)
            {
              // The coordinates of
              // the current entry in
              // DoFHandler
              // numbering, which
              // differs from the
              // block-wise local
              // numbering we use in
              // our local matrices
              const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
              const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);

              // The global dof
              // indices to assemble
              // in. Since we may
              // have face matrices
              // coupling two
              // different cells, we
              // provide two sets of
              // dof indices.
              const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
              const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;

              if (mg_constrained_dofs == 0)
                {
                  if (transpose)
                    global.add(kglobal, jglobal, local(j,k));
                  else
                    global.add(jglobal, kglobal, local(j,k));
                }
              else
                {
                  if (!mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
                      !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
                    {
                      if (mg_constrained_dofs->set_boundary_values())
                        {
                          if ((!mg_constrained_dofs->is_boundary_index(level1, jglobal) &&
                               !mg_constrained_dofs->is_boundary_index(level2, kglobal))
                              ||
                              (mg_constrained_dofs->is_boundary_index(level1, jglobal) &&
                               mg_constrained_dofs->is_boundary_index(level2, kglobal) &&
                               jglobal == kglobal))
                            {
                              if (transpose)
                                global.add(kglobal, jglobal, local(j,k));
                              else
                                global.add(jglobal, kglobal, local(j,k));
                            }
                        }
                      else
                        {
                          if (transpose)
                            global.add(kglobal, jglobal, local(j,k));
                          else
                            global.add(jglobal, kglobal, local(j,k));
                        }
                    }
                }
            }
    }


    template <typename MatrixType, typename number>
    inline void
    MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble_fluxes
    (MatrixType                                 &global,
     const FullMatrix<number>                   &local,
     const unsigned int                          block_row,
     const unsigned int                          block_col,
     const std::vector<types::global_dof_index> &dof1,
     const std::vector<types::global_dof_index> &dof2,
     const unsigned int                          level1,
     const unsigned int                          level2)
    {
      for (unsigned int j=0; j<local.n_rows(); ++j)
        for (unsigned int k=0; k<local.n_cols(); ++k)
          if (std::fabs(local(j,k)) >= threshold)
            {
              // The coordinates of
              // the current entry in
              // DoFHandler
              // numbering, which
              // differs from the
              // block-wise local
              // numbering we use in
              // our local matrices
              const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
              const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);

              // The global dof
              // indices to assemble
              // in. Since we may
              // have face matrices
              // coupling two
              // different cells, we
              // provide two sets of
              // dof indices.
              const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
              const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;

              if (mg_constrained_dofs == 0)
                global.add(jglobal, kglobal, local(j,k));
              else
                {
                  if (!mg_constrained_dofs->non_refinement_edge_index(level1, jglobal) &&
                      !mg_constrained_dofs->non_refinement_edge_index(level2, kglobal))
                    {
                      if (!mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
                          !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
                        global.add(jglobal, kglobal, local(j,k));
                    }
                }
            }
    }

    template <typename MatrixType, typename number>
    inline void
    MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble_up
    (MatrixType                                 &global,
     const FullMatrix<number>                   &local,
     const unsigned int                          block_row,
     const unsigned int                          block_col,
     const std::vector<types::global_dof_index> &dof1,
     const std::vector<types::global_dof_index> &dof2,
     const unsigned int                          level1,
     const unsigned int                          level2)
    {
      for (unsigned int j=0; j<local.n_rows(); ++j)
        for (unsigned int k=0; k<local.n_cols(); ++k)
          if (std::fabs(local(j,k)) >= threshold)
            {
              // The coordinates of
              // the current entry in
              // DoFHandler
              // numbering, which
              // differs from the
              // block-wise local
              // numbering we use in
              // our local matrices
              const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
              const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);

              // The global dof
              // indices to assemble
              // in. Since we may
              // have face matrices
              // coupling two
              // different cells, we
              // provide two sets of
              // dof indices.
              const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
              const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;

              if (mg_constrained_dofs == 0)
                global.add(jglobal, kglobal, local(j,k));
              else
                {
                  if (!mg_constrained_dofs->non_refinement_edge_index(level1, jglobal) &&
                      !mg_constrained_dofs->non_refinement_edge_index(level2, kglobal))
                    {
                      if (!mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
                          !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
                        global.add(jglobal, kglobal, local(j,k));
                    }
                }
            }
    }

    template <typename MatrixType, typename number>
    inline void
    MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble_down
    (MatrixType                                 &global,
     const FullMatrix<number>                   &local,
     const unsigned int                          block_row,
     const unsigned int                          block_col,
     const std::vector<types::global_dof_index> &dof1,
     const std::vector<types::global_dof_index> &dof2,
     const unsigned int                          level1,
     const unsigned int                          level2)
    {
      for (unsigned int j=0; j<local.n_rows(); ++j)
        for (unsigned int k=0; k<local.n_cols(); ++k)
          if (std::fabs(local(k,j)) >= threshold)
            {
              // The coordinates of
              // the current entry in
              // DoFHandler
              // numbering, which
              // differs from the
              // block-wise local
              // numbering we use in
              // our local matrices
              const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
              const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);

              // The global dof
              // indices to assemble
              // in. Since we may
              // have face matrices
              // coupling two
              // different cells, we
              // provide two sets of
              // dof indices.
              const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
              const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;

              if (mg_constrained_dofs == 0)
                global.add(jglobal, kglobal, local(k,j));
              else
                {
                  if (!mg_constrained_dofs->non_refinement_edge_index(level1, jglobal) &&
                      !mg_constrained_dofs->non_refinement_edge_index(level2, kglobal))
                    {
                      if (!mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
                          !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
                        global.add(jglobal, kglobal, local(k,j));
                    }
                }
            }
    }

    template <typename MatrixType, typename number>
    inline void
    MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble_in
    (MatrixType                                 &global,
     const FullMatrix<number>                   &local,
     const unsigned int                          block_row,
     const unsigned int                          block_col,
     const std::vector<types::global_dof_index> &dof1,
     const std::vector<types::global_dof_index> &dof2,
     const unsigned int                          level1,
     const unsigned int                          level2)
    {
//      AssertDimension(local.n(), dof1.size());
//      AssertDimension(local.m(), dof2.size());

      for (unsigned int j=0; j<local.n_rows(); ++j)
        for (unsigned int k=0; k<local.n_cols(); ++k)
          if (std::fabs(local(j,k)) >= threshold)
            {
              // The coordinates of
              // the current entry in
              // DoFHandler
              // numbering, which
              // differs from the
              // block-wise local
              // numbering we use in
              // our local matrices
              const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
              const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);

              // The global dof
              // indices to assemble
              // in. Since we may
              // have face matrices
              // coupling two
              // different cells, we
              // provide two sets of
              // dof indices.
              const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
              const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;

              if (mg_constrained_dofs == 0)
                global.add(jglobal, kglobal, local(j,k));
              else
                {
                  if (mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
                      !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
                    {
                      if (mg_constrained_dofs->set_boundary_values())
                        {
                          if ((!mg_constrained_dofs->is_boundary_index(level1, jglobal) &&
                               !mg_constrained_dofs->is_boundary_index(level2, kglobal))
                              ||
                              (mg_constrained_dofs->is_boundary_index(level1, jglobal) &&
                               mg_constrained_dofs->is_boundary_index(level2, kglobal) &&
                               jglobal == kglobal))
                            global.add(jglobal, kglobal, local(j,k));
                        }
                      else
                        global.add(jglobal, kglobal, local(j,k));
                    }
                }
            }
    }

    template <typename MatrixType, typename number>
    inline void
    MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble_out
    (MatrixType                                 &global,
     const FullMatrix<number>                   &local,
     const unsigned int                          block_row,
     const unsigned int                          block_col,
     const std::vector<types::global_dof_index> &dof1,
     const std::vector<types::global_dof_index> &dof2,
     const unsigned int                          level1,
     const unsigned int                          level2)
    {
//      AssertDimension(local.n(), dof1.size());
//      AssertDimension(local.m(), dof2.size());

      for (unsigned int j=0; j<local.n_rows(); ++j)
        for (unsigned int k=0; k<local.n_cols(); ++k)
          if (std::fabs(local(k,j)) >= threshold)
            {
              // The coordinates of
              // the current entry in
              // DoFHandler
              // numbering, which
              // differs from the
              // block-wise local
              // numbering we use in
              // our local matrices
              const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
              const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);

              // The global dof
              // indices to assemble
              // in. Since we may
              // have face matrices
              // coupling two
              // different cells, we
              // provide two sets of
              // dof indices.
              const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
              const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;

              if (mg_constrained_dofs == 0)
                global.add(jglobal, kglobal, local(k,j));
              else
                {
                  if (mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
                      !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
                    {
                      if (mg_constrained_dofs->set_boundary_values())
                        {
                          if ((!mg_constrained_dofs->is_boundary_index(level1, jglobal) &&
                               !mg_constrained_dofs->is_boundary_index(level2, kglobal))
                              ||
                              (mg_constrained_dofs->is_boundary_index(level1, jglobal) &&
                               mg_constrained_dofs->is_boundary_index(level2, kglobal) &&
                               jglobal == kglobal))
                            global.add(jglobal, kglobal, local(k,j));
                        }
                      else
                        global.add(jglobal, kglobal, local(k,j));
                    }
                }
            }
    }


    template <typename MatrixType, typename number>
    template <class DOFINFO>
    inline void
    MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble(const DOFINFO &info)
    {
      const unsigned int level = info.cell->level();

      for (unsigned int i=0; i<matrices->size(); ++i)
        {
          // Row and column index of
          // the block we are dealing with
          const unsigned int row = matrices->block(i)[level].row;
          const unsigned int col = matrices->block(i)[level].column;

          assemble(matrices->block(i)[level].matrix, info.matrix(i,false).matrix, row, col,
                   info.indices, info.indices, level, level);
          if (mg_constrained_dofs != 0)
            {
              if (interface_in != 0)
                assemble_in(interface_in->block(i)[level], info.matrix(i,false).matrix, row, col,
                            info.indices, info.indices, level, level);
              if (interface_out != 0)
                assemble_in(interface_out->block(i)[level], info.matrix(i,false).matrix, row, col,
                            info.indices, info.indices, level, level);

              assemble_in(matrices->block_in(i)[level], info.matrix(i,false).matrix, row, col,
                          info.indices, info.indices, level, level);
              assemble_out(matrices->block_out(i)[level], info.matrix(i,false).matrix, row, col,
                           info.indices, info.indices, level, level);
            }
        }
    }


    template <typename MatrixType, typename number>
    template <class DOFINFO>
    inline void
    MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble
    (const DOFINFO &info1,
     const DOFINFO &info2)
    {
      const unsigned int level1 = info1.cell->level();
      const unsigned int level2 = info2.cell->level();

      for (unsigned int i=0; i<matrices->size(); ++i)
        {
          MGLevelObject<MatrixBlock<MatrixType> > &o = matrices->block(i);

          // Row and column index of
          // the block we are dealing with
          const unsigned int row = o[level1].row;
          const unsigned int col = o[level1].column;

          if (level1 == level2)
            {
              if (mg_constrained_dofs == 0)
                {
                  assemble(o[level1].matrix, info1.matrix(i,false).matrix, row, col,
                           info1.indices, info1.indices, level1, level1);
                  assemble(o[level1].matrix, info1.matrix(i,true).matrix, row, col,
                           info1.indices, info2.indices, level1, level2);
                  assemble(o[level1].matrix, info2.matrix(i,false).matrix, row, col,
                           info2.indices, info2.indices, level2, level2);
                  assemble(o[level1].matrix, info2.matrix(i,true).matrix, row, col,
                           info2.indices, info1.indices, level2, level1);
                }
              else
                {
                  assemble_fluxes(o[level1].matrix, info1.matrix(i,false).matrix, row, col,
                                  info1.indices, info1.indices, level1, level1);
                  assemble_fluxes(o[level1].matrix, info1.matrix(i,true).matrix, row, col,
                                  info1.indices, info2.indices, level1, level2);
                  assemble_fluxes(o[level1].matrix, info2.matrix(i,false).matrix, row, col,
                                  info2.indices, info2.indices, level2, level2);
                  assemble_fluxes(o[level1].matrix, info2.matrix(i,true).matrix, row, col,
                                  info2.indices, info1.indices, level2, level1);
                }
            }
          else
            {
              Assert(level1 > level2, ExcNotImplemented());
              if (flux_up->size() != 0)
                {
                  // Do not add M22,
                  // which is done by
                  // the coarser cell
                  assemble_fluxes(o[level1].matrix, info1.matrix(i,false).matrix, row, col,
                                  info1.indices, info1.indices, level1, level1);
                  assemble_up(flux_up->block(i)[level1].matrix, info1.matrix(i,true).matrix, row, col,
                              info1.indices, info2.indices, level1, level2);
                  assemble_down(flux_down->block(i)[level1].matrix, info2.matrix(i,true).matrix, row, col,
                                info2.indices, info1.indices, level2, level1);
                }
            }
        }
    }
  }
}

DEAL_II_NAMESPACE_CLOSE

#endif