--- /dev/null
+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)
#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>
};
+/**
+ * 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
}
+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"
for (deal_II_dimension : DIMENSIONS; S1 : REAL_SCALARS)
{
template class MGTransferMatrixFree< deal_II_dimension, S1 >;
+ template class MGTransferBlockMatrixFree< deal_II_dimension, S1 >;
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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();
+ }
+
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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();
+ }
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+}
--- /dev/null
+
+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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+}
--- /dev/null
+
+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::