]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add MGTransferBlockMatrixFree 4676/head
authorDenis Davydov <davydden@gmail.com>
Fri, 4 Aug 2017 08:27:39 +0000 (10:27 +0200)
committerDenis Davydov <davydden@gmail.com>
Fri, 18 Aug 2017 12:56:25 +0000 (14:56 +0200)
12 files changed:
doc/news/changes/minor/20170818DenisDavydov [new file with mode: 0644]
include/deal.II/multigrid/mg_transfer_matrix_free.h
source/multigrid/mg_transfer_matrix_free.cc
source/multigrid/mg_transfer_matrix_free.inst.in
tests/matrix_free/parallel_multigrid_adaptive_06.cc [new file with mode: 0644]
tests/matrix_free/parallel_multigrid_adaptive_06.with_mpi=true.with_p4est=true.mpirun=3.output [new file with mode: 0644]
tests/matrix_free/parallel_multigrid_adaptive_06ref.cc [new file with mode: 0644]
tests/matrix_free/parallel_multigrid_adaptive_06ref.with_mpi=true.with_p4est=true.mpirun=3.output [new file with mode: 0644]
tests/multigrid/transfer_05.cc [new file with mode: 0644]
tests/multigrid/transfer_05.with_mpi=true.with_p4est=true.mpirun=3.output [new file with mode: 0644]
tests/multigrid/transfer_matrix_free_02_block.cc [new file with mode: 0644]
tests/multigrid/transfer_matrix_free_02_block.with_mpi=true.with_p4est=true.with_trilinos=true.mpirun=5.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20170818DenisDavydov b/doc/news/changes/minor/20170818DenisDavydov
new file mode 100644 (file)
index 0000000..f79ad3c
--- /dev/null
@@ -0,0 +1,5 @@
+New: add MGTransferBlockMatrixFree which allows using geometric multi-grids with
+LinearAlgebra::distributed::BlockVector<Number> where each block uses the same
+index space of degrees of freedom.
+<br>
+(Denis Davydov, 2017/08/18)
index 43e3cccc58f94c013ef2de51e8f8eb1919ac7537..75c11578fd6613f4812773c0cdc3d0cccb25c570 100644 (file)
@@ -19,6 +19,7 @@
 #include <deal.II/base/config.h>
 
 #include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/la_parallel_block_vector.h>
 #include <deal.II/multigrid/mg_base.h>
 #include <deal.II/multigrid/mg_constrained_dofs.h>
 #include <deal.II/base/mg_level_object.h>
@@ -234,9 +235,229 @@ private:
 };
 
 
+/**
+ * Implementation of the MGTransferBase interface for which the transfer
+ * operations is implemented in a matrix-free way based on the interpolation
+ * matrices of the underlying finite element. This requires considerably less
+ * memory than MGTransferPrebuilt and can also be considerably faster than
+ * that variant.
+ *
+ * This class works with LinearAlgebra::distributed::BlockVector and
+ * performs exactly the same transfer operations for each block as
+ * MGTransferMatrixFree. This implies that each block should cover the
+ * same index space of DoFs and have the same partitioning.
+ *
+ * @author Denis Davydov
+ * @date 2017
+ */
+template <int dim, typename Number>
+class MGTransferBlockMatrixFree : public MGTransferBase<LinearAlgebra::distributed::BlockVector<Number>>
+{
+public:
+  /**
+   * Constructor without constraint matrices. Use this constructor only with
+   * discontinuous finite elements or with no local refinement.
+   */
+  MGTransferBlockMatrixFree ();
+
+  /**
+   * Constructor with constraints. Equivalent to the default constructor
+   * followed by initialize_constraints().
+   */
+  MGTransferBlockMatrixFree (const MGConstrainedDoFs &mg_constrained_dofs);
+
+  /**
+   * Destructor.
+   */
+  virtual ~MGTransferBlockMatrixFree ();
+
+  /**
+   * Initialize the constraints to be used in build().
+   */
+  void initialize_constraints (const MGConstrainedDoFs &mg_constrained_dofs);
+
+  /**
+   * Reset the object to the state it had right after the default constructor.
+   */
+  void clear ();
+
+  /**
+   * Actually build the information for the prolongation for each level.
+   */
+  void build (const DoFHandler<dim,dim> &mg_dof);
+
+  /**
+   * Prolongate a vector from level <tt>to_level-1</tt> to level
+   * <tt>to_level</tt> using the embedding matrices of the underlying finite
+   * element. The previous content of <tt>dst</tt> is overwritten.
+   *
+   * @param to_level The index of the level to prolongate to, which is the
+   * level of @p dst.
+   *
+   * @param src is a vector with as many elements as there are degrees of
+   * freedom on the coarser level involved.
+   *
+   * @param dst has as many elements as there are degrees of freedom on the
+   * finer level.
+   */
+  virtual void prolongate (const unsigned int                                    to_level,
+                           LinearAlgebra::distributed::BlockVector<Number>       &dst,
+                           const LinearAlgebra::distributed::BlockVector<Number> &src) const;
+
+  /**
+   * Restrict a vector from level <tt>from_level</tt> to level
+   * <tt>from_level-1</tt> using the transpose operation of the prolongate()
+   * method. If the region covered by cells on level <tt>from_level</tt> is
+   * smaller than that of level <tt>from_level-1</tt> (local refinement), then
+   * some degrees of freedom in <tt>dst</tt> are active and will not be
+   * altered. For the other degrees of freedom, the result of the restriction
+   * is added.
+   *
+   * @param from_level The index of the level to restrict from, which is the
+   * level of @p src.
+   *
+   * @param src is a vector with as many elements as there are degrees of
+   * freedom on the finer level involved.
+   *
+   * @param dst has as many elements as there are degrees of freedom on the
+   * coarser level.
+   */
+  virtual void restrict_and_add (const unsigned int from_level,
+                                 LinearAlgebra::distributed::BlockVector<Number>       &dst,
+                                 const LinearAlgebra::distributed::BlockVector<Number> &src) const;
+
+  /**
+   * Transfer from a block-vector on the global grid to block-vectors defined on each of the levels separately.
+   *
+   * This function will initialize @dst accordingly if needed as required by the Multigrid class.
+   */
+  template <typename Number2, int spacedim>
+  void
+  copy_to_mg (const DoFHandler<dim,spacedim> &mg_dof,
+              MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &dst,
+              const LinearAlgebra::distributed::BlockVector<Number2>         &src) const;
+
+  /**
+   * Transfer from multi-level block-vector to normal vector.
+   */
+  template <typename Number2, int spacedim>
+  void
+  copy_from_mg (const DoFHandler<dim,spacedim>                                       &mg_dof,
+                LinearAlgebra::distributed::BlockVector<Number2>                     &dst,
+                const MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &src) const;
+
+  /**
+   * Memory used by this object.
+   */
+  std::size_t memory_consumption () const;
+
+private:
+
+  /**
+   * A non-block matrix-free version of transfer operation.
+   */
+  MGTransferMatrixFree<dim,Number> matrix_free_transfer;
+};
+
+
 /*@}*/
 
 
+//------------------------ templated functions -------------------------
+#ifndef DOXYGEN
+
+template <int dim, typename Number>
+template <typename Number2, int spacedim>
+void
+MGTransferBlockMatrixFree<dim,Number>::
+copy_to_mg (const DoFHandler<dim,spacedim> &mg_dof,
+            MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &dst,
+            const LinearAlgebra::distributed::BlockVector<Number2>         &src) const
+{
+  const unsigned int n_blocks  = src.n_blocks();
+  const unsigned int min_level = dst.min_level();
+  const unsigned int max_level = dst.max_level();
+
+  // this function is normally called within the Multigrid class with
+  // dst == defect level block vector. At first run this vector is not
+  // initialized. Do this below:
+  {
+    const parallel::Triangulation<dim,spacedim> *tria =
+      (dynamic_cast<const parallel::Triangulation<dim,spacedim>*>
+       (&mg_dof.get_triangulation()));
+
+    for (unsigned int level = min_level; level <= max_level; ++level)
+      {
+        dst[level].reinit(n_blocks);
+        bool collect_size = false;
+        for (unsigned int b = 0; b < n_blocks; ++b)
+          {
+            LinearAlgebra::distributed::Vector<Number> &v = dst[level].block(b);
+            if (v.size() != mg_dof.locally_owned_mg_dofs(level).size() ||
+                v.local_size() != mg_dof.locally_owned_mg_dofs(level).n_elements())
+              {
+                v.reinit(mg_dof.locally_owned_mg_dofs(level), tria != nullptr ?
+                         tria->get_communicator() : MPI_COMM_SELF);
+                collect_size = true;
+              }
+            else
+              v = 0.;
+          }
+        if (collect_size)
+          dst[level].collect_sizes ();
+      }
+  }
+
+  // FIXME: this a quite ugly as we need a temporary object:
+  MGLevelObject<LinearAlgebra::distributed::Vector<Number>> dst_non_block(min_level, max_level);
+
+  for (unsigned int b = 0; b < n_blocks; ++b)
+    {
+      for (unsigned int l = min_level; l <= max_level; ++l)
+        dst_non_block[l].reinit(dst[l].block(b));
+
+      matrix_free_transfer.copy_to_mg(mg_dof, dst_non_block, src.block(b));
+
+      for (unsigned int l = min_level; l <= max_level; ++l)
+        dst[l].block(b) = dst_non_block[l];
+    }
+}
+
+template <int dim, typename Number>
+template <typename Number2, int spacedim>
+void
+MGTransferBlockMatrixFree<dim,Number>::
+copy_from_mg (const DoFHandler<dim,spacedim>                        &mg_dof,
+              LinearAlgebra::distributed::BlockVector<Number2>      &dst,
+              const MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &src) const
+{
+  const unsigned int n_blocks  = dst.n_blocks();
+  const unsigned int min_level = src.min_level();
+  const unsigned int max_level = src.max_level();
+
+  for (unsigned int l = min_level; l <= max_level; ++l)
+    AssertDimension(src[l].n_blocks(), dst.n_blocks());
+
+  // FIXME: this a quite ugly as we need a temporary object:
+  MGLevelObject<LinearAlgebra::distributed::Vector<Number>> src_non_block(min_level, max_level);
+
+  for (unsigned int b = 0; b < n_blocks; ++b)
+    {
+      for (unsigned int l = min_level; l <= max_level; ++l)
+        {
+          src_non_block[l].reinit(src[l].block(b));
+          src_non_block[l] = src[l].block(b);
+        }
+
+      matrix_free_transfer.copy_from_mg(mg_dof, dst.block(b), src_non_block);
+    }
+}
+
+
+
+#endif // DOXYGEN
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
index 877c41e258763ed8b43328b1d7177fa0f46f90c6..4ebfdc44921e032437b5f99e2948c681098afa0d 100644 (file)
@@ -553,6 +553,90 @@ MGTransferMatrixFree<dim,Number>::memory_consumption() const
 }
 
 
+template <int dim, typename Number>
+MGTransferBlockMatrixFree<dim,Number>::MGTransferBlockMatrixFree ()
+{}
+
+
+
+template <int dim, typename Number>
+MGTransferBlockMatrixFree<dim,Number>::MGTransferBlockMatrixFree (const MGConstrainedDoFs &mg_c)
+  :
+  matrix_free_transfer(mg_c)
+{
+}
+
+
+
+template <int dim, typename Number>
+MGTransferBlockMatrixFree<dim,Number>::~MGTransferBlockMatrixFree ()
+{}
+
+
+
+template <int dim, typename Number>
+void MGTransferBlockMatrixFree<dim,Number>::initialize_constraints
+(const MGConstrainedDoFs &mg_c)
+{
+  matrix_free_transfer.initialize_constraints(mg_c);
+}
+
+
+
+template <int dim, typename Number>
+void MGTransferBlockMatrixFree<dim,Number>::clear ()
+{
+  matrix_free_transfer.clear();
+}
+
+
+
+template <int dim, typename Number>
+void MGTransferBlockMatrixFree<dim,Number>::build
+(const DoFHandler<dim,dim>  &mg_dof)
+{
+  matrix_free_transfer.build(mg_dof);
+}
+
+
+
+template <int dim, typename Number>
+void MGTransferBlockMatrixFree<dim,Number>
+::prolongate (const unsigned int                           to_level,
+              LinearAlgebra::distributed::BlockVector<Number>       &dst,
+              const LinearAlgebra::distributed::BlockVector<Number> &src) const
+{
+  AssertDimension(dst.n_blocks(), src.n_blocks());
+  const unsigned int n_blocks = dst.n_blocks();
+  for (unsigned int b = 0; b < n_blocks; ++b)
+    matrix_free_transfer.prolongate(to_level, dst.block(b), src.block(b));
+}
+
+
+
+template <int dim, typename Number>
+void MGTransferBlockMatrixFree<dim,Number>
+::restrict_and_add (const unsigned int                           from_level,
+                    LinearAlgebra::distributed::BlockVector<Number>       &dst,
+                    const LinearAlgebra::distributed::BlockVector<Number> &src) const
+{
+  AssertDimension(dst.n_blocks(), src.n_blocks());
+  const unsigned int n_blocks = dst.n_blocks();
+  for (unsigned int b = 0; b < n_blocks; ++b)
+    matrix_free_transfer.restrict_and_add(from_level, dst.block(b), src.block(b));
+}
+
+
+
+template <int dim, typename Number>
+std::size_t
+MGTransferBlockMatrixFree<dim,Number>::memory_consumption() const
+{
+  return matrix_free_transfer.memory_consumption();
+}
+
+
+
 // explicit instantiation
 #include "mg_transfer_matrix_free.inst"
 
index 719e19d8e5c315e6060d6e26ab4f5e55e53cd93c..14d4cecb37d905774aea63add54957dbf5996147 100644 (file)
@@ -18,4 +18,5 @@
 for (deal_II_dimension : DIMENSIONS; S1 : REAL_SCALARS)
 {
     template class MGTransferMatrixFree< deal_II_dimension, S1 >;
+    template class MGTransferBlockMatrixFree< deal_II_dimension, S1 >;
 }
diff --git a/tests/matrix_free/parallel_multigrid_adaptive_06.cc b/tests/matrix_free/parallel_multigrid_adaptive_06.cc
new file mode 100644 (file)
index 0000000..30ad01e
--- /dev/null
@@ -0,0 +1,399 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 - 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// similar to parallel_multigrid_adaptive_06ref but using MGTransferBlockMatrixFree
+// to solve block-diagonal matrix with Laplace operator on diagonals.
+// As expected, when we use a block vector with a single block, we get
+// the same results as the reference, non-block solution, i.e.
+//
+// DEAL:2d:cg::Starting value 21.93
+// DEAL:2d:cg::Convergence step 7 value 1.961e-07
+//
+//
+// What this test is really for is block operations. As expected we see
+// exactly the same number of iterations and the ratio of \sqrt(2) in values, i.e.
+//
+// DEAL:2d:cg::Starting value 31.02
+// DEAL:2d:cg::Convergence step 7 value 2.774e-07
+
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/multigrid/multigrid.h>
+#include <deal.II/multigrid/mg_transfer.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+#include <deal.II/multigrid/mg_tools.h>
+#include <deal.II/multigrid/mg_coarse.h>
+#include <deal.II/multigrid/mg_smoother.h>
+#include <deal.II/multigrid/mg_matrix.h>
+
+#include <deal.II/matrix_free/operators.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+std::ofstream logfile("output");
+
+using namespace dealii::MatrixFreeOperators;
+
+/**
+ * Block-Laplace operator with Block vector.
+ */
+template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1, typename BlockVectorType = LinearAlgebra::distributed::BlockVector<double> >
+class BlockLaplace : public Subscriptor
+{
+public:
+  typedef typename BlockVectorType::value_type value_type;
+  typedef typename BlockVectorType::size_type size_type;
+
+  BlockLaplace ()
+    : Subscriptor()
+  {
+  }
+
+  void  initialize (std::shared_ptr<const MatrixFree<dim,value_type>> data)
+  {
+    laplace.initialize(data);
+  }
+
+  void  initialize (std::shared_ptr<const MatrixFree<dim,value_type>> data, const MGConstrainedDoFs &mg_constrained_dofs,
+                    const unsigned int level)
+  {
+    laplace.initialize(data, mg_constrained_dofs, level);
+  }
+
+  void  vmult_interface_down (BlockVectorType &dst, const BlockVectorType &src) const
+  {
+    for (unsigned int b = 0; b < src.n_blocks(); ++b)
+      laplace.vmult_interface_down(dst.block(b), src.block(b));
+  }
+
+  void  vmult_interface_up (BlockVectorType &dst, const BlockVectorType &src) const
+  {
+    for (unsigned int b = 0; b < src.n_blocks(); ++b)
+      laplace.vmult_interface_up(dst.block(b), src.block(b));
+  }
+
+  void  vmult (BlockVectorType &dst, const BlockVectorType &src) const
+  {
+    for (unsigned int b = 0; b < src.n_blocks(); ++b)
+      laplace.vmult(dst.block(b), src.block(b));
+  }
+
+  void  Tvmult (BlockVectorType &dst, const BlockVectorType &src) const
+  {
+    for (unsigned int b = 0; b < src.n_blocks(); ++b)
+      laplace.Tvmult(dst.block(b), src.block(b));
+  }
+
+  void  vmult_add (BlockVectorType &dst, const BlockVectorType &src) const
+  {
+    for (unsigned int b = 0; b < src.n_blocks(); ++b)
+      laplace.vmult_add(dst.block(b), src.block(b));
+  }
+
+  void  Tvmult_add (BlockVectorType &dst, const BlockVectorType &src) const
+  {
+    for (unsigned int b = 0; b < src.n_blocks(); ++b)
+      laplace.Tvmult_add(dst.block(b), src.block(b));
+  }
+
+  void precondition_Jacobi (BlockVectorType &dst, const BlockVectorType &src, const value_type omega) const
+  {
+    for (unsigned int b = 0; b < src.n_blocks(); ++b)
+      laplace.precondition_Jacobi(dst.block(b), src.block(b), omega);
+  }
+
+  void compute_diagonal ()
+  {
+    laplace.compute_diagonal ();
+  }
+
+
+  virtual void clear()
+  {
+    laplace.clear();
+  }
+
+private:
+
+  MatrixFreeOperators::LaplaceOperator<dim, fe_degree, n_q_points_1d, 1, typename BlockVectorType::BlockType> laplace;
+
+};
+
+
+
+template <typename MatrixType, typename Number>
+class MGCoarseIterative : public MGCoarseGridBase<LinearAlgebra::distributed::BlockVector<Number> >
+{
+public:
+  MGCoarseIterative() {}
+
+  void initialize(const MatrixType &matrix)
+  {
+    coarse_matrix = &matrix;
+  }
+
+  virtual void operator() (const unsigned int   level,
+                           LinearAlgebra::distributed::BlockVector<Number> &dst,
+                           const LinearAlgebra::distributed::BlockVector<Number> &src) const
+  {
+    ReductionControl solver_control (1e4, 1e-50, 1e-10);
+    SolverCG<LinearAlgebra::distributed::BlockVector<Number> > solver_coarse (solver_control);
+    solver_coarse.solve (*coarse_matrix, dst, src, PreconditionIdentity());
+  }
+
+  const MatrixType *coarse_matrix;
+};
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d, typename number>
+void do_test (const DoFHandler<dim>  &dof, const unsigned int nb)
+{
+  if (types_are_equal<number,float>::value == true)
+    {
+      deallog.push("float");
+    }
+  else
+    {
+    }
+
+  deallog << "Testing " << dof.get_fe().get_name();
+  deallog << std::endl;
+  deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
+
+  IndexSet locally_relevant_dofs;
+  DoFTools::extract_locally_relevant_dofs(dof, locally_relevant_dofs);
+
+  // Dirichlet BC
+  ZeroFunction<dim> zero_function;
+  typename FunctionMap<dim>::type dirichlet_boundary;
+  dirichlet_boundary[0] = &zero_function;
+
+  // fine-level constraints
+  ConstraintMatrix constraints;
+  constraints.reinit(locally_relevant_dofs);
+  DoFTools::make_hanging_node_constraints(dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, dirichlet_boundary,
+                                           constraints);
+  constraints.close();
+
+  // level constraints:
+  MGConstrainedDoFs mg_constrained_dofs;
+  mg_constrained_dofs.initialize(dof, dirichlet_boundary);
+
+  MappingQ<dim> mapping(fe_degree+1);
+
+  BlockLaplace<dim,fe_degree,n_q_points_1d,LinearAlgebra::distributed::BlockVector<number>> fine_matrix;
+  std::shared_ptr<MatrixFree<dim,number> > fine_level_data(new MatrixFree<dim,number> ());
+
+  typename MatrixFree<dim,number>::AdditionalData fine_level_additional_data;
+  fine_level_additional_data.tasks_parallel_scheme = MatrixFree<dim,number>::AdditionalData::none;
+  fine_level_additional_data.tasks_block_size = 3;
+  fine_level_data->reinit (mapping, dof, constraints, QGauss<1>(n_q_points_1d),
+                           fine_level_additional_data);
+
+  fine_matrix.initialize(fine_level_data);
+  fine_matrix.compute_diagonal();
+
+  LinearAlgebra::distributed::BlockVector<number> in(nb), sol(nb);
+  for (unsigned int b=0; b < nb; ++b)
+    {
+      fine_level_data->initialize_dof_vector(in.block(b));
+      fine_level_data->initialize_dof_vector(sol.block(b));
+    }
+
+  in.collect_sizes();
+  sol.collect_sizes();
+
+  // set constant rhs vector
+  {
+    // this is to make it consistent with parallel_multigrid_adaptive.cc
+    ConstraintMatrix hanging_node_constraints;
+    hanging_node_constraints.reinit(locally_relevant_dofs);
+    DoFTools::make_hanging_node_constraints(dof, hanging_node_constraints);
+    hanging_node_constraints.close();
+
+    for (unsigned int i=0; i<in.block(0).local_size(); ++i)
+      if (!hanging_node_constraints.is_constrained(in.block(0).get_partitioner()->local_to_global(i)))
+        in.block(0).local_element(i) = 1.;
+
+    for (unsigned int b = 1; b < nb; ++b)
+      in.block(b) = in.block(0);
+  }
+
+  // set up multigrid in analogy to step-37
+  typedef BlockLaplace<dim,fe_degree,n_q_points_1d,LinearAlgebra::distributed::BlockVector<number>> LevelMatrixType;
+
+  MGLevelObject<LevelMatrixType> mg_matrices;
+  MGLevelObject<MatrixFree<dim,number> > mg_level_data;
+  mg_matrices.resize(0, dof.get_triangulation().n_global_levels()-1);
+  mg_level_data.resize(0, dof.get_triangulation().n_global_levels()-1);
+  for (unsigned int level = 0; level<dof.get_triangulation().n_global_levels(); ++level)
+    {
+      typename MatrixFree<dim,number>::AdditionalData mg_additional_data;
+      mg_additional_data.tasks_parallel_scheme = MatrixFree<dim,number>::AdditionalData::none;
+      mg_additional_data.tasks_block_size = 3;
+      mg_additional_data.level_mg_handler = level;
+
+      ConstraintMatrix level_constraints;
+      IndexSet relevant_dofs;
+      DoFTools::extract_locally_relevant_level_dofs(dof, level,
+                                                    relevant_dofs);
+      level_constraints.reinit(relevant_dofs);
+      level_constraints.add_lines(mg_constrained_dofs.get_boundary_indices(level));
+      level_constraints.close();
+
+      mg_level_data[level].reinit (mapping, dof, level_constraints, QGauss<1>(n_q_points_1d),
+                                   mg_additional_data);
+      mg_matrices[level].initialize(std::make_shared<MatrixFree<dim,number> >(
+                                      mg_level_data[level]),
+                                    mg_constrained_dofs,
+                                    level);
+      mg_matrices[level].compute_diagonal();
+    }
+
+  MGLevelObject<MGInterfaceOperator<LevelMatrixType> > mg_interface_matrices;
+  mg_interface_matrices.resize(0, dof.get_triangulation().n_global_levels()-1);
+  for (unsigned int level=0; level<dof.get_triangulation().n_global_levels(); ++level)
+    mg_interface_matrices[level].initialize(mg_matrices[level]);
+
+  MGTransferBlockMatrixFree<dim,number> mg_transfer(mg_constrained_dofs);
+  mg_transfer.build(dof);
+
+  MGCoarseIterative<LevelMatrixType,number> mg_coarse;
+  mg_coarse.initialize(mg_matrices[0]);
+
+  typedef PreconditionJacobi<LevelMatrixType> SMOOTHER;
+  MGSmootherPrecondition<LevelMatrixType, SMOOTHER, LinearAlgebra::distributed::BlockVector<number> >
+  mg_smoother;
+
+  mg_smoother.initialize(mg_matrices, typename SMOOTHER::AdditionalData(0.8));
+
+  mg::Matrix<LinearAlgebra::distributed::BlockVector<number> >
+  mg_matrix(mg_matrices);
+  mg::Matrix<LinearAlgebra::distributed::BlockVector<number> >
+  mg_interface(mg_interface_matrices);
+
+  Multigrid<LinearAlgebra::distributed::BlockVector<number> > mg(
+    mg_matrix,
+    mg_coarse,
+    mg_transfer,
+    mg_smoother,
+    mg_smoother);
+  mg.set_edge_matrices(mg_interface, mg_interface);
+  PreconditionMG<dim, LinearAlgebra::distributed::BlockVector<number>,
+                 MGTransferBlockMatrixFree<dim,number>>
+                 preconditioner(dof, mg, mg_transfer);
+
+  {
+    // avoid output from inner (coarse-level) solver
+    deallog.depth_file(3);
+
+    ReductionControl control(30, 1e-20, 1e-7);
+    SolverCG<LinearAlgebra::distributed::BlockVector<number> > solver(control);
+    solver.solve(fine_matrix, sol, in, preconditioner);
+  }
+
+  if (types_are_equal<number,float>::value == true)
+    deallog.pop();
+
+  fine_matrix.clear();
+  for (unsigned int level = 0; level<dof.get_triangulation().n_global_levels(); ++level)
+    mg_matrices[level].clear();
+}
+
+
+
+template <int dim, int fe_degree>
+void test (const unsigned int nbands = 1)
+{
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD,
+                                                 Triangulation<dim>::limit_level_difference_at_vertices,
+                                                 parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+  GridGenerator::hyper_cube (tria);
+  tria.refine_global(6-dim);
+  const unsigned int n_runs = fe_degree == 1 ? 6-dim : 5-dim;
+  for (unsigned int i=0; i<n_runs; ++i)
+    {
+      for (typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active();
+           cell != tria.end(); ++cell)
+        if (cell->is_locally_owned() &&
+            (cell->center().norm() < 0.5 && (cell->level() < 5 ||
+                                             cell->center().norm() > 0.45)
+             ||
+             (dim == 2 && cell->center().norm() > 1.2)))
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+      FE_Q<dim> fe (fe_degree);
+      DoFHandler<dim> dof (tria);
+      dof.distribute_dofs(fe);
+      dof.distribute_mg_dofs(fe);
+
+      do_test<dim, fe_degree, fe_degree+1, double> (dof, nbands);
+    }
+}
+
+
+
+int main (int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    {
+      deallog.attach(logfile);
+      deallog << std::setprecision (4);
+    }
+
+  {
+    deallog.push("2d");
+    test<2,1>();
+    test<2,3>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3,1>();
+    test<3,2>();
+    deallog.pop();
+  }
+
+  // 2 blocks
+  {
+    deallog.push("2d");
+    test<2,1>(2);
+    test<2,3>(2);
+    deallog.pop();
+    deallog.push("3d");
+    test<3,1>(2);
+    test<3,2>(2);
+    deallog.pop();
+  }
+
+}
diff --git a/tests/matrix_free/parallel_multigrid_adaptive_06.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/matrix_free/parallel_multigrid_adaptive_06.with_mpi=true.with_p4est=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..322c18c
--- /dev/null
@@ -0,0 +1,97 @@
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 507
+DEAL:2d:cg::Starting value 21.93
+DEAL:2d:cg::Convergence step 7 value 1.961e-07
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 881
+DEAL:2d:cg::Starting value 27.77
+DEAL:2d:cg::Convergence step 8 value 2.956e-07
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 2147
+DEAL:2d:cg::Starting value 43.26
+DEAL:2d:cg::Convergence step 7 value 2.063e-06
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 6517
+DEAL:2d:cg::Starting value 76.92
+DEAL:2d:cg::Convergence step 7 value 7.253e-06
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 4259
+DEAL:2d:cg::Starting value 64.26
+DEAL:2d:cg::Convergence step 9 value 1.536e-06
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 7457
+DEAL:2d:cg::Starting value 83.11
+DEAL:2d:cg::Convergence step 9 value 8.264e-06
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 18535
+DEAL:2d:cg::Starting value 131.0
+DEAL:2d:cg::Convergence step 9 value 1.071e-05
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 1103
+DEAL:3d:cg::Starting value 31.21
+DEAL:3d:cg::Convergence step 8 value 3.711e-07
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 3694
+DEAL:3d:cg::Starting value 55.00
+DEAL:3d:cg::Convergence step 8 value 5.151e-06
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 9566
+DEAL:3d:cg::Starting value 76.18
+DEAL:3d:cg::Convergence step 9 value 6.124e-06
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Number of degrees of freedom: 7494
+DEAL:3d:cg::Starting value 82.90
+DEAL:3d:cg::Convergence step 9 value 3.485e-06
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Number of degrees of freedom: 26548
+DEAL:3d:cg::Starting value 152.6
+DEAL:3d:cg::Convergence step 10 value 2.494e-06
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 507
+DEAL:2d:cg::Starting value 31.02
+DEAL:2d:cg::Convergence step 7 value 2.774e-07
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 881
+DEAL:2d:cg::Starting value 39.27
+DEAL:2d:cg::Convergence step 8 value 4.180e-07
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 2147
+DEAL:2d:cg::Starting value 61.17
+DEAL:2d:cg::Convergence step 7 value 2.918e-06
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 6517
+DEAL:2d:cg::Starting value 108.8
+DEAL:2d:cg::Convergence step 7 value 1.026e-05
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 4259
+DEAL:2d:cg::Starting value 90.87
+DEAL:2d:cg::Convergence step 9 value 2.172e-06
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 7457
+DEAL:2d:cg::Starting value 117.5
+DEAL:2d:cg::Convergence step 9 value 1.169e-05
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 18535
+DEAL:2d:cg::Starting value 185.2
+DEAL:2d:cg::Convergence step 9 value 1.515e-05
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 1103
+DEAL:3d:cg::Starting value 44.14
+DEAL:3d:cg::Convergence step 8 value 5.248e-07
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 3694
+DEAL:3d:cg::Starting value 77.78
+DEAL:3d:cg::Convergence step 8 value 7.285e-06
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 9566
+DEAL:3d:cg::Starting value 107.7
+DEAL:3d:cg::Convergence step 9 value 8.661e-06
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Number of degrees of freedom: 7494
+DEAL:3d:cg::Starting value 117.2
+DEAL:3d:cg::Convergence step 9 value 4.929e-06
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Number of degrees of freedom: 26548
+DEAL:3d:cg::Starting value 215.8
+DEAL:3d:cg::Convergence step 10 value 3.527e-06
diff --git a/tests/matrix_free/parallel_multigrid_adaptive_06ref.cc b/tests/matrix_free/parallel_multigrid_adaptive_06ref.cc
new file mode 100644 (file)
index 0000000..2823d82
--- /dev/null
@@ -0,0 +1,315 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// same as parallel_multigrid_adaptive_02, but with Jacobi smoother.
+// The only reason for this test is to provide a reference for
+// parallel_multigrid_adaptive_06
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/multigrid/multigrid.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+#include <deal.II/multigrid/mg_tools.h>
+#include <deal.II/multigrid/mg_coarse.h>
+#include <deal.II/multigrid/mg_smoother.h>
+#include <deal.II/multigrid/mg_matrix.h>
+
+#include <deal.II/matrix_free/operators.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+std::ofstream logfile("output");
+
+using namespace dealii::MatrixFreeOperators;
+
+template <int dim, typename LAPLACEOPERATOR>
+class MGTransferMF : public MGTransferMatrixFree<dim, typename LAPLACEOPERATOR::value_type>
+{
+public:
+  MGTransferMF(const MGLevelObject<LAPLACEOPERATOR> &laplace,
+               const MGConstrainedDoFs &mg_constrained_dofs)
+    :
+    MGTransferMatrixFree<dim, typename LAPLACEOPERATOR::value_type>(mg_constrained_dofs),
+    laplace_operator (laplace)
+  {
+  }
+
+  /**
+   * Overload copy_to_mg from MGTransferPrebuilt to get the vectors compatible
+   * with MatrixFree and bypass the crude vector initialization in
+   * MGTransferPrebuilt
+   */
+  template <class InVector, int spacedim>
+  void
+  copy_to_mg (const DoFHandler<dim,spacedim> &mg_dof_handler,
+              MGLevelObject<LinearAlgebra::distributed::Vector<typename LAPLACEOPERATOR::value_type> > &dst,
+              const InVector &src) const
+  {
+    for (unsigned int level=dst.min_level();
+         level<=dst.max_level(); ++level)
+      laplace_operator[level].initialize_dof_vector(dst[level]);
+    MGTransferMatrixFree<dim, typename LAPLACEOPERATOR::value_type>::
+    copy_to_mg(mg_dof_handler, dst, src);
+  }
+
+private:
+  const MGLevelObject<LAPLACEOPERATOR> &laplace_operator;
+};
+
+
+
+template <typename MatrixType, typename Number>
+class MGCoarseIterative : public MGCoarseGridBase<LinearAlgebra::distributed::Vector<Number> >
+{
+public:
+  MGCoarseIterative() {}
+
+  void initialize(const MatrixType &matrix)
+  {
+    coarse_matrix = &matrix;
+  }
+
+  virtual void operator() (const unsigned int   level,
+                           LinearAlgebra::distributed::Vector<double> &dst,
+                           const LinearAlgebra::distributed::Vector<double> &src) const
+  {
+    ReductionControl solver_control (1e4, 1e-50, 1e-10);
+    SolverCG<LinearAlgebra::distributed::Vector<double> > solver_coarse (solver_control);
+    solver_coarse.solve (*coarse_matrix, dst, src, PreconditionIdentity());
+  }
+
+  const MatrixType *coarse_matrix;
+};
+
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d, typename number>
+void do_test (const DoFHandler<dim>  &dof)
+{
+  if (types_are_equal<number,float>::value == true)
+    {
+      deallog.push("float");
+    }
+  else
+    {
+    }
+
+  deallog << "Testing " << dof.get_fe().get_name();
+  deallog << std::endl;
+  deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
+
+  IndexSet locally_relevant_dofs;
+  DoFTools::extract_locally_relevant_dofs(dof, locally_relevant_dofs);
+
+  // Dirichlet BC
+  ZeroFunction<dim> zero_function;
+  typename FunctionMap<dim>::type dirichlet_boundary;
+  dirichlet_boundary[0] = &zero_function;
+
+  // fine-level constraints
+  ConstraintMatrix constraints;
+  constraints.reinit(locally_relevant_dofs);
+  DoFTools::make_hanging_node_constraints(dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, dirichlet_boundary,
+                                           constraints);
+  constraints.close();
+
+  // level constraints:
+  MGConstrainedDoFs mg_constrained_dofs;
+  mg_constrained_dofs.initialize(dof, dirichlet_boundary);
+
+  MappingQ<dim> mapping(fe_degree+1);
+
+  LaplaceOperator<dim,fe_degree,n_q_points_1d,1,LinearAlgebra::distributed::Vector<number> > fine_matrix;
+  std::shared_ptr<MatrixFree<dim,number> > fine_level_data(new MatrixFree<dim,number> ());
+
+  typename MatrixFree<dim,number>::AdditionalData fine_level_additional_data;
+  fine_level_additional_data.tasks_parallel_scheme = MatrixFree<dim,number>::AdditionalData::none;
+  fine_level_additional_data.tasks_block_size = 3;
+  fine_level_data->reinit (mapping, dof, constraints, QGauss<1>(n_q_points_1d),
+                           fine_level_additional_data);
+
+  fine_matrix.initialize(fine_level_data);
+  fine_matrix.compute_diagonal();
+
+
+  LinearAlgebra::distributed::Vector<number> in, sol;
+  fine_matrix.initialize_dof_vector(in);
+  fine_matrix.initialize_dof_vector(sol);
+
+  // set constant rhs vector
+  {
+    // this is to make it consistent with parallel_multigrid_adaptive.cc
+    ConstraintMatrix hanging_node_constraints;
+    hanging_node_constraints.reinit(locally_relevant_dofs);
+    DoFTools::make_hanging_node_constraints(dof, hanging_node_constraints);
+    hanging_node_constraints.close();
+
+    for (unsigned int i=0; i<in.local_size(); ++i)
+      if (!hanging_node_constraints.is_constrained(in.get_partitioner()->local_to_global(i)))
+        in.local_element(i) = 1.;
+  }
+
+  // set up multigrid in analogy to step-37
+  typedef LaplaceOperator<dim,fe_degree,n_q_points_1d,1,LinearAlgebra::distributed::Vector<number> > LevelMatrixType;
+
+  MGLevelObject<LevelMatrixType> mg_matrices;
+  MGLevelObject<MatrixFree<dim,number> > mg_level_data;
+  mg_matrices.resize(0, dof.get_triangulation().n_global_levels()-1);
+  mg_level_data.resize(0, dof.get_triangulation().n_global_levels()-1);
+  for (unsigned int level = 0; level<dof.get_triangulation().n_global_levels(); ++level)
+    {
+      typename MatrixFree<dim,number>::AdditionalData mg_additional_data;
+      mg_additional_data.tasks_parallel_scheme = MatrixFree<dim,number>::AdditionalData::none;
+      mg_additional_data.tasks_block_size = 3;
+      mg_additional_data.level_mg_handler = level;
+
+      ConstraintMatrix level_constraints;
+      IndexSet relevant_dofs;
+      DoFTools::extract_locally_relevant_level_dofs(dof, level,
+                                                    relevant_dofs);
+      level_constraints.reinit(relevant_dofs);
+      level_constraints.add_lines(mg_constrained_dofs.get_boundary_indices(level));
+      level_constraints.close();
+
+      mg_level_data[level].reinit (mapping, dof, level_constraints, QGauss<1>(n_q_points_1d),
+                                   mg_additional_data);
+      mg_matrices[level].initialize(std::make_shared<MatrixFree<dim,number> >(
+                                      mg_level_data[level]),
+                                    mg_constrained_dofs,
+                                    level);
+      mg_matrices[level].compute_diagonal();
+    }
+  MGLevelObject<MGInterfaceOperator<LevelMatrixType> > mg_interface_matrices;
+  mg_interface_matrices.resize(0, dof.get_triangulation().n_global_levels()-1);
+  for (unsigned int level=0; level<dof.get_triangulation().n_global_levels(); ++level)
+    mg_interface_matrices[level].initialize(mg_matrices[level]);
+
+  MGTransferMF<dim,LevelMatrixType> mg_transfer(mg_matrices,
+                                                mg_constrained_dofs);
+  mg_transfer.build(dof);
+
+  MGCoarseIterative<LevelMatrixType,number> mg_coarse;
+  mg_coarse.initialize(mg_matrices[0]);
+
+  typedef PreconditionJacobi<LevelMatrixType> SMOOTHER;
+  MGSmootherPrecondition<LevelMatrixType, SMOOTHER, LinearAlgebra::distributed::Vector<double> >
+  mg_smoother;
+
+  mg_smoother.initialize(mg_matrices, typename SMOOTHER::AdditionalData(0.8));
+
+  mg::Matrix<LinearAlgebra::distributed::Vector<double> >
+  mg_matrix(mg_matrices);
+  mg::Matrix<LinearAlgebra::distributed::Vector<double> >
+  mg_interface(mg_interface_matrices);
+
+  Multigrid<LinearAlgebra::distributed::Vector<double> > mg(dof,
+                                                            mg_matrix,
+                                                            mg_coarse,
+                                                            mg_transfer,
+                                                            mg_smoother,
+                                                            mg_smoother);
+  mg.set_edge_matrices(mg_interface, mg_interface);
+  PreconditionMG<dim, LinearAlgebra::distributed::Vector<double>,
+                 MGTransferMF<dim,LevelMatrixType> >
+                 preconditioner(dof, mg, mg_transfer);
+
+  {
+    // avoid output from inner (coarse-level) solver
+    deallog.depth_file(3);
+
+    ReductionControl control(30, 1e-20, 1e-7);
+    SolverCG<LinearAlgebra::distributed::Vector<double> > solver(control);
+    solver.solve(fine_matrix, sol, in, preconditioner);
+  }
+
+  if (types_are_equal<number,float>::value == true)
+    deallog.pop();
+
+  fine_matrix.clear();
+  for (unsigned int level = 0; level<dof.get_triangulation().n_global_levels(); ++level)
+    mg_matrices[level].clear();
+}
+
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD,
+                                                 Triangulation<dim>::limit_level_difference_at_vertices,
+                                                 parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+  GridGenerator::hyper_cube (tria);
+  tria.refine_global(6-dim);
+  const unsigned int n_runs = fe_degree == 1 ? 6-dim : 5-dim;
+  for (unsigned int i=0; i<n_runs; ++i)
+    {
+      for (typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active();
+           cell != tria.end(); ++cell)
+        if (cell->is_locally_owned() &&
+            (cell->center().norm() < 0.5 && (cell->level() < 5 ||
+                                             cell->center().norm() > 0.45)
+             ||
+             (dim == 2 && cell->center().norm() > 1.2)))
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+      FE_Q<dim> fe (fe_degree);
+      DoFHandler<dim> dof (tria);
+      dof.distribute_dofs(fe);
+      dof.distribute_mg_dofs(fe);
+
+      do_test<dim, fe_degree, fe_degree+1, double> (dof);
+    }
+}
+
+
+
+int main (int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    {
+      deallog.attach(logfile);
+      deallog << std::setprecision (4);
+    }
+
+  {
+    deallog.push("2d");
+    test<2,1>();
+    test<2,3>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3,1>();
+    test<3,2>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/parallel_multigrid_adaptive_06ref.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/matrix_free/parallel_multigrid_adaptive_06ref.with_mpi=true.with_p4est=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..8ca8213
--- /dev/null
@@ -0,0 +1,49 @@
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 507
+DEAL:2d:cg::Starting value 21.93
+DEAL:2d:cg::Convergence step 7 value 1.961e-07
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 881
+DEAL:2d:cg::Starting value 27.77
+DEAL:2d:cg::Convergence step 8 value 2.956e-07
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 2147
+DEAL:2d:cg::Starting value 43.26
+DEAL:2d:cg::Convergence step 7 value 2.063e-06
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 6517
+DEAL:2d:cg::Starting value 76.92
+DEAL:2d:cg::Convergence step 7 value 7.253e-06
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 4259
+DEAL:2d:cg::Starting value 64.26
+DEAL:2d:cg::Convergence step 9 value 1.536e-06
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 7457
+DEAL:2d:cg::Starting value 83.11
+DEAL:2d:cg::Convergence step 9 value 8.264e-06
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 18535
+DEAL:2d:cg::Starting value 131.0
+DEAL:2d:cg::Convergence step 9 value 1.071e-05
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 1103
+DEAL:3d:cg::Starting value 31.21
+DEAL:3d:cg::Convergence step 8 value 3.711e-07
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 3694
+DEAL:3d:cg::Starting value 55.00
+DEAL:3d:cg::Convergence step 8 value 5.151e-06
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 9566
+DEAL:3d:cg::Starting value 76.18
+DEAL:3d:cg::Convergence step 9 value 6.124e-06
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Number of degrees of freedom: 7494
+DEAL:3d:cg::Starting value 82.90
+DEAL:3d:cg::Convergence step 9 value 3.485e-06
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Number of degrees of freedom: 26548
+DEAL:3d:cg::Starting value 152.6
+DEAL:3d:cg::Convergence step 10 value 2.494e-06
diff --git a/tests/multigrid/transfer_05.cc b/tests/multigrid/transfer_05.cc
new file mode 100644 (file)
index 0000000..8a927f7
--- /dev/null
@@ -0,0 +1,198 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2006 - 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+// check MGTransferBlockMatrixFree::copy_[from|to]_mg by comparision to non-block
+// MGTransferMatrixFree
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/function_map.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/multigrid/mg_tools.h>
+#include <deal.II/multigrid/mg_transfer.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
+#include <algorithm>
+
+using namespace std;
+
+
+template <int dim, typename Number>
+void check(const unsigned int fe_degree)
+{
+  FE_Q<dim> fe(fe_degree);
+  deallog << "FE: " << fe.get_name() << std::endl;
+
+  // run a few different sizes...
+  unsigned int sizes [] = {1, 2, 3};
+  for (unsigned int cycle=0; cycle<sizeof(sizes)/sizeof(unsigned int); ++cycle)
+    {
+      unsigned int n_refinements = 0;
+      unsigned int n_subdiv = sizes[cycle];
+      if (n_subdiv > 1)
+        while (n_subdiv%2 == 0)
+          {
+            n_refinements += 1;
+            n_subdiv /= 2;
+          }
+      n_refinements += 3-dim;
+      if (fe_degree < 3)
+        n_refinements += 1;
+
+      parallel::distributed::Triangulation<dim>
+      tr(MPI_COMM_WORLD,
+         Triangulation<dim>::limit_level_difference_at_vertices,
+         parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+      GridGenerator::subdivided_hyper_cube(tr, n_subdiv);
+      tr.refine_global(n_refinements);
+
+      // adaptive refinement into a circle
+      for (typename Triangulation<dim>::active_cell_iterator cell=tr.begin_active(); cell != tr.end(); ++cell)
+        if (cell->is_locally_owned() &&
+            cell->center().norm() < 0.5)
+          cell->set_refine_flag();
+      tr.execute_coarsening_and_refinement();
+      for (typename Triangulation<dim>::active_cell_iterator cell=tr.begin_active(); cell != tr.end(); ++cell)
+        if (cell->is_locally_owned() &&
+            cell->center().norm() > 0.3 && cell->center().norm() < 0.4)
+          cell->set_refine_flag();
+      tr.execute_coarsening_and_refinement();
+      for (typename Triangulation<dim>::active_cell_iterator cell=tr.begin_active(); cell != tr.end(); ++cell)
+        if (cell->is_locally_owned() &&
+            cell->center().norm() > 0.33 && cell->center().norm() < 0.37)
+          cell->set_refine_flag();
+      tr.execute_coarsening_and_refinement();
+
+      deallog << "no. cells: " << tr.n_global_active_cells() << std::endl;
+
+      DoFHandler<dim> mgdof(tr);
+      mgdof.distribute_dofs(fe);
+      mgdof.distribute_mg_dofs(fe);
+
+      MGConstrainedDoFs mg_constrained_dofs;
+      ZeroFunction<dim> zero_function;
+      typename FunctionMap<dim>::type dirichlet_boundary;
+      dirichlet_boundary[0] = &zero_function;
+      mg_constrained_dofs.initialize(mgdof, dirichlet_boundary);
+
+      // build reference non-block
+      MGTransferMatrixFree<dim, Number> transfer_ref(mg_constrained_dofs);
+      transfer_ref.build(mgdof);
+
+      // build matrix-free block transfer
+      MGTransferBlockMatrixFree<dim, Number> transfer(mg_constrained_dofs);
+      transfer.build(mgdof);
+
+      const unsigned int nb = 3;
+
+      MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> lbv(0, tr.n_global_levels()-1);
+      LinearAlgebra::distributed::BlockVector<Number> bv(nb);
+
+      MGLevelObject<LinearAlgebra::distributed::Vector<Number>> lv(0, tr.n_global_levels()-1);
+      LinearAlgebra::distributed::Vector<Number> v(nb);
+
+      // initialize
+      v.reinit(mgdof.locally_owned_dofs(), MPI_COMM_WORLD);
+      for (unsigned int b = 0; b < nb; ++b)
+        bv.block(b).reinit(mgdof.locally_owned_dofs(), MPI_COMM_WORLD);
+
+      for (unsigned int l=lbv.min_level(); l<=lbv.max_level(); ++l)
+        {
+          lv[l].reinit(mgdof.locally_owned_mg_dofs(l), MPI_COMM_WORLD);
+
+          lbv[l].reinit(nb);
+          for (unsigned int b = 0; b < nb; ++b)
+            lbv[l].block(b).reinit(mgdof.locally_owned_mg_dofs(l), MPI_COMM_WORLD);
+
+          lbv[l].collect_sizes();
+
+          // set values:
+          for (unsigned int b = 0; b < nb; ++b)
+            for (unsigned int i=0; i<lbv[l].block(b).local_size(); ++i)
+              lbv[l].block(b).local_element(i) = (double)Testing::rand()/RAND_MAX;
+
+          lbv[l].compress(VectorOperation::insert);
+        }
+
+      // check copy_from_mg
+      transfer.copy_from_mg(mgdof, bv, lbv);
+      for (unsigned int b=0; b < nb; ++b)
+        {
+          for (unsigned int l=lv.min_level(); l<=lv.max_level(); ++l)
+            lv[l] = lbv[l].block(b);
+
+          transfer_ref.copy_from_mg(mgdof, v, lv);
+          v -= bv.block(b);
+          deallog << "Diff copy_from_mg b" << b << ": " << v.l2_norm() << std::endl;
+        }
+
+      // check copy_to_mg
+      // use un-initialized level-block vector to make sure that it will be
+      // set correctly in copy_to_mg().
+      MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> lbv2(0, tr.n_global_levels()-1);
+      for (unsigned int b=0; b < nb; ++b)
+        for (unsigned int i = 0; i < bv.block(b).local_size(); ++i)
+          bv.block(b).local_element(i) = (double)Testing::rand()/RAND_MAX;
+
+      transfer.copy_to_mg(mgdof, lbv2, bv);
+      // Also check that the block vector has its (global) size set on each level:
+      for (unsigned int l=lv.min_level(); l<=lv.max_level(); ++l)
+        {
+          unsigned int total_size = 0;
+          for (unsigned int b=0; b < nb; ++b)
+            total_size += lbv2[l].block(b).size();
+
+          AssertThrow (total_size == lbv2[l].size(),
+                       ExcDimensionMismatch(total_size,lbv2[l].size()));
+        }
+
+      // Finally check the difference:
+      for (unsigned int b=0; b < nb; ++b)
+        {
+          v = bv.block(b);
+          transfer_ref.copy_to_mg(mgdof, lv, v);
+          for (unsigned int l=lv.min_level(); l<=lv.max_level(); ++l)
+            {
+              lv[l] -= lbv2[l].block(b);
+              deallog << "Diff copy_to_mg   l" << l << ": " << lv[l].l2_norm() << std::endl;
+            }
+        }
+
+      deallog << std::endl;
+    }
+}
+
+
+int main(int argc, char **argv)
+{
+  // no threading in this test...
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+  mpi_initlog();
+
+  check<2,double>(1);
+  check<2,double>(3);
+  check<3,double>(1);
+  check<3,double>(3);
+  check<2,float> (2);
+  check<3,float> (2);
+}
diff --git a/tests/multigrid/transfer_05.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/multigrid/transfer_05.with_mpi=true.with_p4est=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..bace0b4
--- /dev/null
@@ -0,0 +1,367 @@
+
+DEAL::FE: FE_Q<2>(1)
+DEAL::no. cells: 88
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::
+DEAL::no. cells: 250
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l6: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l6: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l6: 0.00000
+DEAL::
+DEAL::no. cells: 510
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::
+DEAL::FE: FE_Q<2>(3)
+DEAL::no. cells: 34
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::
+DEAL::no. cells: 88
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::
+DEAL::no. cells: 141
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::
+DEAL::FE: FE_Q<3>(1)
+DEAL::no. cells: 15
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::
+DEAL::no. cells: 694
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::
+DEAL::no. cells: 1791
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::
+DEAL::FE: FE_Q<3>(3)
+DEAL::no. cells: 1
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::
+DEAL::no. cells: 15
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::
+DEAL::no. cells: 223
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::
+DEAL::FE: FE_Q<2>(2)
+DEAL::no. cells: 88
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::
+DEAL::no. cells: 250
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l6: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l6: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l6: 0.00000
+DEAL::
+DEAL::no. cells: 510
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::
+DEAL::FE: FE_Q<3>(2)
+DEAL::no. cells: 15
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::
+DEAL::no. cells: 694
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l5: 0.00000
+DEAL::
+DEAL::no. cells: 1791
+DEAL::Diff copy_from_mg b0: 0.00000
+DEAL::Diff copy_from_mg b1: 0.00000
+DEAL::Diff copy_from_mg b2: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::Diff copy_to_mg   l0: 0.00000
+DEAL::Diff copy_to_mg   l1: 0.00000
+DEAL::Diff copy_to_mg   l2: 0.00000
+DEAL::Diff copy_to_mg   l3: 0.00000
+DEAL::Diff copy_to_mg   l4: 0.00000
+DEAL::
diff --git a/tests/multigrid/transfer_matrix_free_02_block.cc b/tests/multigrid/transfer_matrix_free_02_block.cc
new file mode 100644 (file)
index 0000000..1b3915c
--- /dev/null
@@ -0,0 +1,166 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Check MGTransferBlockMatrixFree by comparison with MGTransferMatrixFree on a
+// series of adaptive meshes for FE_Q. This is, essentially, a
+// copy-paste of transfer_matrix_free_02.cc
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/la_parallel_block_vector.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/multigrid/mg_transfer.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+
+
+template <int dim, typename Number>
+void check(const unsigned int fe_degree)
+{
+  FE_Q<dim> fe(fe_degree);
+  deallog << "FE: " << fe.get_name() << std::endl;
+
+  // run a few different sizes...
+  unsigned int sizes [] = {1, 2, 3};
+  for (unsigned int cycle=0; cycle<sizeof(sizes)/sizeof(unsigned int); ++cycle)
+    {
+      unsigned int n_refinements = 0;
+      unsigned int n_subdiv = sizes[cycle];
+      if (n_subdiv > 1)
+        while (n_subdiv%2 == 0)
+          {
+            n_refinements += 1;
+            n_subdiv /= 2;
+          }
+      n_refinements += 3-dim;
+      if (fe_degree < 3)
+        n_refinements += 1;
+
+      parallel::distributed::Triangulation<dim>
+      tr(MPI_COMM_WORLD,
+         Triangulation<dim>::limit_level_difference_at_vertices,
+         parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+      GridGenerator::subdivided_hyper_cube(tr, n_subdiv);
+      tr.refine_global(n_refinements);
+
+      // adaptive refinement into a circle
+      for (typename Triangulation<dim>::active_cell_iterator cell=tr.begin_active(); cell != tr.end(); ++cell)
+        if (cell->is_locally_owned() &&
+            cell->center().norm() < 0.5)
+          cell->set_refine_flag();
+      tr.execute_coarsening_and_refinement();
+      for (typename Triangulation<dim>::active_cell_iterator cell=tr.begin_active(); cell != tr.end(); ++cell)
+        if (cell->is_locally_owned() &&
+            cell->center().norm() > 0.3 && cell->center().norm() < 0.4)
+          cell->set_refine_flag();
+      tr.execute_coarsening_and_refinement();
+      for (typename Triangulation<dim>::active_cell_iterator cell=tr.begin_active(); cell != tr.end(); ++cell)
+        if (cell->is_locally_owned() &&
+            cell->center().norm() > 0.33 && cell->center().norm() < 0.37)
+          cell->set_refine_flag();
+      tr.execute_coarsening_and_refinement();
+
+      deallog << "no. cells: " << tr.n_global_active_cells() << std::endl;
+
+      DoFHandler<dim> mgdof(tr);
+      mgdof.distribute_dofs(fe);
+      mgdof.distribute_mg_dofs(fe);
+
+      MGConstrainedDoFs mg_constrained_dofs;
+      ZeroFunction<dim> zero_function;
+      typename FunctionMap<dim>::type dirichlet_boundary;
+      dirichlet_boundary[0] = &zero_function;
+      mg_constrained_dofs.initialize(mgdof, dirichlet_boundary);
+
+      // build reference
+      MGTransferMatrixFree<dim, Number> transfer_ref(mg_constrained_dofs);
+      transfer_ref.build(mgdof);
+
+      // build matrix-free transfer
+      MGTransferBlockMatrixFree<dim, Number> transfer(mg_constrained_dofs);
+      transfer.build(mgdof);
+
+      const unsigned int nb = 3;
+      // check prolongation for all levels using random vector
+      for (unsigned int level=1; level<mgdof.get_triangulation().n_global_levels(); ++level)
+        {
+          LinearAlgebra::distributed::BlockVector<Number> v1(nb), v2(nb), v3(nb);
+          for (unsigned int b = 0; b < nb; ++b)
+            {
+              v1.block(b).reinit(mgdof.locally_owned_mg_dofs(level-1), MPI_COMM_WORLD);
+              v2.block(b).reinit(mgdof.locally_owned_mg_dofs(level), MPI_COMM_WORLD);
+              v3.block(b).reinit(mgdof.locally_owned_mg_dofs(level), MPI_COMM_WORLD);
+
+              for (unsigned int i=0; i<v1.block(b).local_size(); ++i)
+                v1.block(b).local_element(i) = (double)Testing::rand()/RAND_MAX;
+
+              transfer_ref.prolongate(level, v2.block(b), v1.block(b));
+            }
+
+          transfer.prolongate(level, v3, v1);
+          v3 -= v2;
+          deallog << "Diff prolongate   l" << level << ": " << v3.l2_norm() << std::endl;
+        }
+
+      // check restriction for all levels using random vector
+      for (unsigned int level=1; level<mgdof.get_triangulation().n_global_levels(); ++level)
+        {
+          LinearAlgebra::distributed::BlockVector<Number> v1(nb), v2(nb), v3(nb);
+          for (unsigned int b = 0; b < nb; ++b)
+            {
+              v1.block(b).reinit(mgdof.locally_owned_mg_dofs(level), MPI_COMM_WORLD);
+              v2.block(b).reinit(mgdof.locally_owned_mg_dofs(level-1), MPI_COMM_WORLD);
+              v3.block(b).reinit(mgdof.locally_owned_mg_dofs(level-1), MPI_COMM_WORLD);
+
+              for (unsigned int i=0; i<v1.block(b).local_size(); ++i)
+                v1.block(b).local_element(i) = (double)Testing::rand()/RAND_MAX;
+
+              transfer_ref.restrict_and_add(level, v2.block(b), v1.block(b));
+            }
+
+          transfer.restrict_and_add(level, v3, v1);
+          v3 -= v2;
+          deallog << "Diff restrict     l" << level << ": " << v3.l2_norm() << std::endl;
+
+          v2 = 1.;
+          v3 = 1.;
+          transfer.restrict_and_add(level, v2, v1);
+          for (unsigned int b = 0; b < nb; ++b)
+            transfer_ref.restrict_and_add(level, v3.block(b), v1.block(b));
+          v3 -= v2;
+          deallog << "Diff restrict add l" << level << ": " << v3.l2_norm() << std::endl;
+        }
+      deallog << std::endl;
+    }
+}
+
+
+int main(int argc, char **argv)
+{
+  // no threading in this test...
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+  mpi_initlog();
+
+  check<2,double>(1);
+  check<2,double>(3);
+  check<3,double>(1);
+  check<3,double>(3);
+  check<2,float> (2);
+  check<3,float> (2);
+}
diff --git a/tests/multigrid/transfer_matrix_free_02_block.with_mpi=true.with_p4est=true.with_trilinos=true.mpirun=5.output b/tests/multigrid/transfer_matrix_free_02_block.with_mpi=true.with_p4est=true.with_trilinos=true.mpirun=5.output
new file mode 100644 (file)
index 0000000..40a8160
--- /dev/null
@@ -0,0 +1,259 @@
+
+DEAL::FE: FE_Q<2>(1)
+DEAL::no. cells: 88
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff prolongate   l3: 0
+DEAL::Diff prolongate   l4: 0
+DEAL::Diff prolongate   l5: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::Diff restrict     l3: 0
+DEAL::Diff restrict add l3: 0
+DEAL::Diff restrict     l4: 0
+DEAL::Diff restrict add l4: 0
+DEAL::Diff restrict     l5: 0
+DEAL::Diff restrict add l5: 0
+DEAL::
+DEAL::no. cells: 250
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff prolongate   l3: 0
+DEAL::Diff prolongate   l4: 0
+DEAL::Diff prolongate   l5: 0
+DEAL::Diff prolongate   l6: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::Diff restrict     l3: 0
+DEAL::Diff restrict add l3: 0
+DEAL::Diff restrict     l4: 0
+DEAL::Diff restrict add l4: 0
+DEAL::Diff restrict     l5: 0
+DEAL::Diff restrict add l5: 0
+DEAL::Diff restrict     l6: 0
+DEAL::Diff restrict add l6: 0
+DEAL::
+DEAL::no. cells: 510
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff prolongate   l3: 0
+DEAL::Diff prolongate   l4: 0
+DEAL::Diff prolongate   l5: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::Diff restrict     l3: 0
+DEAL::Diff restrict add l3: 0
+DEAL::Diff restrict     l4: 0
+DEAL::Diff restrict add l4: 0
+DEAL::Diff restrict     l5: 0
+DEAL::Diff restrict add l5: 0
+DEAL::
+DEAL::FE: FE_Q<2>(3)
+DEAL::no. cells: 34
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff prolongate   l3: 0
+DEAL::Diff prolongate   l4: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::Diff restrict     l3: 0
+DEAL::Diff restrict add l3: 0
+DEAL::Diff restrict     l4: 0
+DEAL::Diff restrict add l4: 0
+DEAL::
+DEAL::no. cells: 88
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff prolongate   l3: 0
+DEAL::Diff prolongate   l4: 0
+DEAL::Diff prolongate   l5: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::Diff restrict     l3: 0
+DEAL::Diff restrict add l3: 0
+DEAL::Diff restrict     l4: 0
+DEAL::Diff restrict add l4: 0
+DEAL::Diff restrict     l5: 0
+DEAL::Diff restrict add l5: 0
+DEAL::
+DEAL::no. cells: 141
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff prolongate   l3: 0
+DEAL::Diff prolongate   l4: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::Diff restrict     l3: 0
+DEAL::Diff restrict add l3: 0
+DEAL::Diff restrict     l4: 0
+DEAL::Diff restrict add l4: 0
+DEAL::
+DEAL::FE: FE_Q<3>(1)
+DEAL::no. cells: 15
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::
+DEAL::no. cells: 694
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff prolongate   l3: 0
+DEAL::Diff prolongate   l4: 0
+DEAL::Diff prolongate   l5: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::Diff restrict     l3: 0
+DEAL::Diff restrict add l3: 0
+DEAL::Diff restrict     l4: 0
+DEAL::Diff restrict add l4: 0
+DEAL::Diff restrict     l5: 0
+DEAL::Diff restrict add l5: 0
+DEAL::
+DEAL::no. cells: 1791
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff prolongate   l3: 0
+DEAL::Diff prolongate   l4: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::Diff restrict     l3: 0
+DEAL::Diff restrict add l3: 0
+DEAL::Diff restrict     l4: 0
+DEAL::Diff restrict add l4: 0
+DEAL::
+DEAL::FE: FE_Q<3>(3)
+DEAL::no. cells: 1
+DEAL::
+DEAL::no. cells: 15
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::
+DEAL::no. cells: 223
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff prolongate   l3: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::Diff restrict     l3: 0
+DEAL::Diff restrict add l3: 0
+DEAL::
+DEAL::FE: FE_Q<2>(2)
+DEAL::no. cells: 88
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff prolongate   l3: 0
+DEAL::Diff prolongate   l4: 0
+DEAL::Diff prolongate   l5: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::Diff restrict     l3: 0
+DEAL::Diff restrict add l3: 0
+DEAL::Diff restrict     l4: 0
+DEAL::Diff restrict add l4: 0
+DEAL::Diff restrict     l5: 0
+DEAL::Diff restrict add l5: 0
+DEAL::
+DEAL::no. cells: 250
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff prolongate   l3: 0
+DEAL::Diff prolongate   l4: 0
+DEAL::Diff prolongate   l5: 0
+DEAL::Diff prolongate   l6: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::Diff restrict     l3: 0
+DEAL::Diff restrict add l3: 0
+DEAL::Diff restrict     l4: 0
+DEAL::Diff restrict add l4: 0
+DEAL::Diff restrict     l5: 0
+DEAL::Diff restrict add l5: 0
+DEAL::Diff restrict     l6: 0
+DEAL::Diff restrict add l6: 0
+DEAL::
+DEAL::no. cells: 510
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff prolongate   l3: 0
+DEAL::Diff prolongate   l4: 0
+DEAL::Diff prolongate   l5: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::Diff restrict     l3: 0
+DEAL::Diff restrict add l3: 0
+DEAL::Diff restrict     l4: 0
+DEAL::Diff restrict add l4: 0
+DEAL::Diff restrict     l5: 0
+DEAL::Diff restrict add l5: 0
+DEAL::
+DEAL::FE: FE_Q<3>(2)
+DEAL::no. cells: 15
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::
+DEAL::no. cells: 694
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff prolongate   l3: 0
+DEAL::Diff prolongate   l4: 0
+DEAL::Diff prolongate   l5: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::Diff restrict     l3: 0
+DEAL::Diff restrict add l3: 0
+DEAL::Diff restrict     l4: 0
+DEAL::Diff restrict add l4: 0
+DEAL::Diff restrict     l5: 0
+DEAL::Diff restrict add l5: 0
+DEAL::
+DEAL::no. cells: 1791
+DEAL::Diff prolongate   l1: 0
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff prolongate   l3: 0
+DEAL::Diff prolongate   l4: 0
+DEAL::Diff restrict     l1: 0
+DEAL::Diff restrict add l1: 0
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::Diff restrict     l3: 0
+DEAL::Diff restrict add l3: 0
+DEAL::Diff restrict     l4: 0
+DEAL::Diff restrict add l4: 0
+DEAL::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.