--- /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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// This test is similar to parallel_multigrid_adaptive_06 but we also test
+// for different polynomial degree in different blocks.
+// We expect to have the same iteration numbers as in
+// parallel_multigrid_adaptive_06 with repsect to the highest polynomial
+// degree used.
+
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/std_cxx14/algorithm.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_1, int fe_degree_2, int n_q_points_1d, 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)
+ {
+ laplace1.initialize(data, std::vector<unsigned int>(1,0));
+ laplace2.initialize(data, std::vector<unsigned int>(1,1));
+ }
+
+ void initialize (std::shared_ptr<const MatrixFree<dim,value_type>> data,
+ const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
+ const unsigned int level)
+ {
+ laplace1.initialize(data, mg_constrained_dofs[0], level, std::vector<unsigned int>(1,0));
+ laplace2.initialize(data, mg_constrained_dofs[1], level, std::vector<unsigned int>(1,1));
+ }
+
+ void vmult_interface_down (BlockVectorType &dst, const BlockVectorType &src) const
+ {
+ laplace1.vmult_interface_down(dst.block(0), src.block(0));
+ laplace2.vmult_interface_down(dst.block(1), src.block(1));
+ }
+
+ void vmult_interface_up (BlockVectorType &dst, const BlockVectorType &src) const
+ {
+ laplace1.vmult_interface_up(dst.block(0), src.block(0));
+ laplace2.vmult_interface_up(dst.block(1), src.block(1));
+ }
+
+ void vmult (BlockVectorType &dst, const BlockVectorType &src) const
+ {
+ laplace1.vmult(dst.block(0), src.block(0));
+ laplace2.vmult(dst.block(1), src.block(1));
+ }
+
+ void Tvmult (BlockVectorType &dst, const BlockVectorType &src) const
+ {
+ laplace1.Tvmult(dst.block(0), src.block(0));
+ laplace2.Tvmult(dst.block(1), src.block(1));
+ }
+
+ void vmult_add (BlockVectorType &dst, const BlockVectorType &src) const
+ {
+ laplace1.vmult_add(dst.block(0), src.block(0));
+ laplace2.vmult_add(dst.block(1), src.block(1));
+ }
+
+ void Tvmult_add (BlockVectorType &dst, const BlockVectorType &src) const
+ {
+ laplace1.Tvmult_add(dst.block(0), src.block(0));
+ laplace2.Tvmult_add(dst.block(1), src.block(1));
+ }
+
+ void precondition_Jacobi (BlockVectorType &dst, const BlockVectorType &src, const value_type omega) const
+ {
+ laplace1.precondition_Jacobi(dst.block(0), src.block(0), omega);
+ laplace2.precondition_Jacobi(dst.block(1), src.block(1), omega);
+ }
+
+ void compute_diagonal ()
+ {
+ laplace1.compute_diagonal ();
+ laplace2.compute_diagonal ();
+ }
+
+
+ virtual void clear()
+ {
+ laplace1.clear();
+ laplace2.clear();
+ }
+
+private:
+
+ MatrixFreeOperators::LaplaceOperator<dim, fe_degree_1, n_q_points_1d, 1, typename BlockVectorType::BlockType> laplace1;
+ MatrixFreeOperators::LaplaceOperator<dim, fe_degree_2, n_q_points_1d, 1, typename BlockVectorType::BlockType> laplace2;
+
+};
+
+
+
+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_1, int fe_degree_2, int n_q_points_1d, typename number>
+void do_test (const std::vector<const DoFHandler<dim>*> &dof)
+{
+ const unsigned int nb = 2;
+ if (types_are_equal<number,float>::value == true)
+ {
+ deallog.push("float");
+ }
+ else
+ {
+ }
+
+ for (unsigned int i=0; i<dof.size(); ++i)
+ {
+ deallog << "Testing " << dof[i]->get_fe().get_name();
+ deallog << std::endl;
+ deallog << "Number of degrees of freedom: " << dof[i]->n_dofs() << std::endl;
+ }
+
+ std::vector<IndexSet> locally_relevant_dofs(dof.size());
+ for (unsigned int i=0; i<dof.size(); ++i)
+ DoFTools::extract_locally_relevant_dofs(*dof[i], locally_relevant_dofs[i]);
+
+ // Dirichlet BC
+ ZeroFunction<dim> zero_function;
+ typename FunctionMap<dim>::type dirichlet_boundary;
+ dirichlet_boundary[0] = &zero_function;
+
+ // fine-level constraints
+ std::vector<ConstraintMatrix> constraints(dof.size());
+ std::vector<const ConstraintMatrix *> constraints_ptrs(dof.size());
+ for (unsigned int i=0; i<dof.size(); ++i)
+ {
+ constraints[i].reinit(locally_relevant_dofs[i]);
+ DoFTools::make_hanging_node_constraints(*dof[i], constraints[i]);
+ VectorTools::interpolate_boundary_values(*dof[i], dirichlet_boundary,
+ constraints[i]);
+ constraints[i].close();
+ constraints_ptrs[i] = &constraints[i];
+ }
+ QGauss<1> quad(n_q_points_1d);
+ constexpr int max_degree = std_cxx14::max(fe_degree_1, fe_degree_2);
+ MappingQ<dim> mapping(max_degree);
+
+ 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;
+
+ std::shared_ptr<MatrixFree<dim, double> >fine_level_data (new MatrixFree<dim, double>());
+ fine_level_data->reinit (mapping, dof, constraints_ptrs, quad, fine_level_additional_data);
+
+ BlockLaplace<dim,fe_degree_1,fe_degree_2,n_q_points_1d,LinearAlgebra::distributed::BlockVector<number>> fine_matrix;
+
+ 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),b);
+ fine_level_data->initialize_dof_vector(sol.block(b),b);
+ }
+
+ in.collect_sizes();
+ sol.collect_sizes();
+
+ // set constant rhs vector
+ {
+ for (unsigned int b = 0; b < nb; ++b)
+ {
+ // this is to make it consistent with parallel_multigrid_adaptive.cc
+ ConstraintMatrix hanging_node_constraints;
+
+ hanging_node_constraints.reinit(locally_relevant_dofs[b]);
+ DoFTools::make_hanging_node_constraints(*dof[b], hanging_node_constraints);
+ hanging_node_constraints.close();
+
+ for (unsigned int i = 0; i < in.block(b).local_size(); ++i)
+ if (!hanging_node_constraints.is_constrained(in.block(b).get_partitioner()->local_to_global(i)))
+ in.block(b).local_element(i) = 1.;
+ }
+ }
+
+ // level constraints:
+ std::vector<MGConstrainedDoFs> mg_constrained_dofs(dof.size());
+ for (unsigned int i=0; i<dof.size(); ++i)
+ mg_constrained_dofs[i].initialize(*dof[i], dirichlet_boundary);
+
+ // set up multigrid in analogy to step-37
+ typedef BlockLaplace<dim,fe_degree_1,fe_degree_2,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[0]->get_triangulation().n_global_levels()-1);
+ mg_level_data.resize(0, dof[0]->get_triangulation().n_global_levels()-1);
+ for (unsigned int level = 0; level<dof[0]->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;
+
+ std::vector<ConstraintMatrix> level_constraints(dof.size());
+ std::vector<const ConstraintMatrix * > level_constraints_ptrs(dof.size());
+ for (unsigned int i=0; i<dof.size(); ++i)
+ {
+ IndexSet relevant_dofs;
+ DoFTools::extract_locally_relevant_level_dofs(*dof[i], level, relevant_dofs);
+ level_constraints[i].reinit(relevant_dofs);
+ level_constraints[i].add_lines(mg_constrained_dofs[i].get_boundary_indices(level));
+ level_constraints[i].close();
+ level_constraints_ptrs[i] = &level_constraints[i];
+ }
+
+ mg_level_data[level].reinit (mapping, dof, level_constraints_ptrs, quad,
+ 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[0]->get_triangulation().n_global_levels()-1);
+ for (unsigned int level=0; level<dof[0]->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[0]->get_triangulation().n_global_levels(); ++level)
+ mg_matrices[level].clear();
+}
+
+
+
+template <int dim, int fe_degree_1, int fe_degree_2>
+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);
+ constexpr int max_degree = std_cxx14::max(fe_degree_1, fe_degree_2);
+ const unsigned int n_runs = max_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_1 (fe_degree_1);
+ DoFHandler<dim> dof_1 (tria);
+ dof_1.distribute_dofs(fe_1);
+ dof_1.distribute_mg_dofs();
+
+ FE_Q<dim> fe_2 (fe_degree_2);
+ DoFHandler<dim> dof_2 (tria);
+ dof_2.distribute_dofs(fe_2);
+ dof_2.distribute_mg_dofs();
+
+ std::vector<const DoFHandler<dim, dim> *> dh_ptrs {&dof_1, &dof_2};
+
+ constexpr int n_q_points_1d = std_cxx14::max(fe_degree_1, fe_degree_2)+1;
+ do_test<dim, fe_degree_1, fe_degree_2, n_q_points_1d, double> (dh_ptrs);
+ }
+}
+
+
+
+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);
+ }
+
+ // 2 blocks
+ {
+ deallog.push("2d");
+ test<2,1,1>();
+ deallog << std::endl;
+ test<2,1,3>();
+ deallog << std::endl;
+ test<2,3,1>();
+ deallog.pop();
+ deallog << std::endl;
+ deallog.push("3d");
+ test<3,1,1>();
+ deallog << std::endl;
+ test<3,1,2>();
+ deallog << std::endl;
+ test<3,2,1>();
+ deallog.pop();
+ }
+
+}
--- /dev/null
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 507
+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::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::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::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::
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 507
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 4259
+DEAL:2d:cg::Starting value 67.90
+DEAL:2d:cg::Convergence step 9 value 1.556e-06
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 881
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 7457
+DEAL:2d:cg::Starting value 87.62
+DEAL:2d:cg::Convergence step 9 value 8.294e-06
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 2147
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 18535
+DEAL:2d:cg::Starting value 137.9
+DEAL:2d:cg::Convergence step 9 value 1.073e-05
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 4259
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 507
+DEAL:2d:cg::Starting value 67.90
+DEAL:2d:cg::Convergence step 9 value 1.556e-06
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 7457
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 881
+DEAL:2d:cg::Starting value 87.62
+DEAL:2d:cg::Convergence step 9 value 8.294e-06
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 18535
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 2147
+DEAL:2d:cg::Starting value 137.9
+DEAL:2d:cg::Convergence step 9 value 1.073e-05
+DEAL::
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 1103
+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::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::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::
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 1103
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Number of degrees of freedom: 7494
+DEAL:3d:cg::Starting value 88.58
+DEAL:3d:cg::Convergence step 9 value 3.537e-06
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 3694
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Number of degrees of freedom: 26548
+DEAL:3d:cg::Starting value 162.2
+DEAL:3d:cg::Convergence step 10 value 2.508e-06
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Number of degrees of freedom: 7494
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 1103
+DEAL:3d:cg::Starting value 88.58
+DEAL:3d:cg::Convergence step 9 value 3.537e-06
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Number of degrees of freedom: 26548
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 3694
+DEAL:3d:cg::Starting value 162.2
+DEAL:3d:cg::Convergence step 10 value 2.508e-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. This is similar to transfer_05 but using
+// a separate DoFHandler for each block
+
+#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)
+{
+ const unsigned int nb = 2;
+
+ FE_Q<dim> fe_1(fe_degree+1);
+ FE_Q<dim> fe_2(fe_degree);
+
+ deallog << "FE: " << fe_1.get_name() << std::endl;
+ deallog << "FE: " << fe_2.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_1(tr);
+ mgdof_1.distribute_dofs(fe_1);
+ mgdof_1.distribute_mg_dofs(fe_1);
+
+ DoFHandler<dim> mgdof_2(tr);
+ mgdof_2.distribute_dofs(fe_2);
+ mgdof_2.distribute_mg_dofs(fe_2);
+
+ const std::vector<const DoFHandler<dim>*> mgdof_ptr {&mgdof_1, &mgdof_2};
+
+ std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(2);
+ ZeroFunction<dim> zero_function;
+ typename FunctionMap<dim>::type dirichlet_boundary;
+ dirichlet_boundary[0] = &zero_function;
+ for (unsigned int i=0; i< mgdof_ptr.size(); ++i)
+ mg_constrained_dofs_vector[i].initialize(*mgdof_ptr[i], dirichlet_boundary);
+
+ // build reference non-block
+ std::vector<MGTransferMatrixFree<dim, Number> > transfer_ref;
+ for (unsigned int b=0; b<nb; ++b)
+ {
+ transfer_ref.emplace_back(mg_constrained_dofs_vector[b]);
+ transfer_ref[b].build(*mgdof_ptr[b]);
+ }
+
+ // build matrix-free block transfer
+ MGTransferBlockMatrixFree<dim, Number> transfer(mg_constrained_dofs_vector);
+ transfer.build(mgdof_ptr);
+
+ 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
+ for (unsigned int b = 0; b < nb; ++b)
+ bv.block(b).reinit(mgdof_ptr[b]->locally_owned_dofs(), MPI_COMM_WORLD);
+
+ for (unsigned int l=lbv.min_level(); l<=lbv.max_level(); ++l)
+ {
+ lbv[l].reinit(nb);
+ for (unsigned int b = 0; b < nb; ++b)
+ lbv[l].block(b).reinit(mgdof_ptr[b]->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_ptr, bv, lbv);
+ for (unsigned int b=0; b < nb; ++b)
+ {
+ v.reinit(mgdof_ptr[b]->locally_owned_dofs(), MPI_COMM_WORLD);
+
+ for (unsigned int l=lv.min_level(); l<=lv.max_level(); ++l)
+ {
+ lv[l].reinit(mgdof_ptr[b]->locally_owned_mg_dofs(l), MPI_COMM_WORLD);
+ lv[l] = lbv[l].block(b);
+ }
+
+ transfer_ref[b].copy_from_mg(*mgdof_ptr[b], 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_ptr, 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.reinit(mgdof_ptr[b]->locally_owned_dofs(), MPI_COMM_WORLD);
+
+ v = bv.block(b);
+ transfer_ref[b].copy_to_mg(*mgdof_ptr[b], 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>(2)
+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_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_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_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>(4)
+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_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_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_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>(2)
+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_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_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_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>(4)
+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_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_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_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>(3)
+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_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_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_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>(3)
+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_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_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_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 similar to transfer_matrix_free_02.cc
+// but using separate DoFHandlers for each block.
+
+#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)
+{
+ const unsigned int nb = 2;
+
+ FE_Q<dim> fe_1(fe_degree+1);
+ FE_Q<dim> fe_2(fe_degree);
+
+ deallog << "FE: " << fe_1.get_name() << std::endl;
+ deallog << "FE: " << fe_2.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_1(tr);
+ mgdof_1.distribute_dofs(fe_1);
+ mgdof_1.distribute_mg_dofs(fe_1);
+
+ DoFHandler<dim> mgdof_2(tr);
+ mgdof_2.distribute_dofs(fe_2);
+ mgdof_2.distribute_mg_dofs(fe_2);
+
+ const std::vector<const DoFHandler<dim>*> mgdof_ptr {&mgdof_1, &mgdof_2};
+
+ std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(2);
+ ZeroFunction<dim> zero_function;
+ typename FunctionMap<dim>::type dirichlet_boundary;
+ dirichlet_boundary[0] = &zero_function;
+ for (unsigned int i=0; i< mgdof_ptr.size(); ++i)
+ mg_constrained_dofs_vector[i].initialize(*mgdof_ptr[i], dirichlet_boundary);
+
+ // build reference
+ std::vector<MGTransferMatrixFree<dim, Number> > transfer_ref;
+ for (unsigned int b=0; b<nb; ++b)
+ {
+ transfer_ref.emplace_back(mg_constrained_dofs_vector[b]);
+ transfer_ref[b].build(*mgdof_ptr[b]);
+ }
+
+ // build matrix-free transfer
+ MGTransferBlockMatrixFree<dim, Number> transfer(mg_constrained_dofs_vector);
+ transfer.build(mgdof_ptr);
+
+ // check prolongation for all levels using random vector
+ for (unsigned int level=1; level<mgdof_ptr[0]->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_ptr[b]->locally_owned_mg_dofs(level-1), MPI_COMM_WORLD);
+ v2.block(b).reinit(mgdof_ptr[b]->locally_owned_mg_dofs(level), MPI_COMM_WORLD);
+ v3.block(b).reinit(mgdof_ptr[b]->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[b].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_ptr[0]->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_ptr[b]->locally_owned_mg_dofs(level), MPI_COMM_WORLD);
+ v2.block(b).reinit(mgdof_ptr[b]->locally_owned_mg_dofs(level-1), MPI_COMM_WORLD);
+ v3.block(b).reinit(mgdof_ptr[b]->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[b].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[b].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>(2)
+DEAL::FE: FE_Q<2>(1)
+DEAL::no. cells: 88
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff prolongate l3: 0.00000
+DEAL::Diff prolongate l4: 0.00000
+DEAL::Diff prolongate l5: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::Diff restrict l3: 0.00000
+DEAL::Diff restrict add l3: 0.00000
+DEAL::Diff restrict l4: 0.00000
+DEAL::Diff restrict add l4: 0.00000
+DEAL::Diff restrict l5: 0.00000
+DEAL::Diff restrict add l5: 0.00000
+DEAL::
+DEAL::no. cells: 250
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff prolongate l3: 0.00000
+DEAL::Diff prolongate l4: 0.00000
+DEAL::Diff prolongate l5: 0.00000
+DEAL::Diff prolongate l6: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::Diff restrict l3: 0.00000
+DEAL::Diff restrict add l3: 0.00000
+DEAL::Diff restrict l4: 0.00000
+DEAL::Diff restrict add l4: 0.00000
+DEAL::Diff restrict l5: 0.00000
+DEAL::Diff restrict add l5: 0.00000
+DEAL::Diff restrict l6: 0.00000
+DEAL::Diff restrict add l6: 0.00000
+DEAL::
+DEAL::no. cells: 510
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff prolongate l3: 0.00000
+DEAL::Diff prolongate l4: 0.00000
+DEAL::Diff prolongate l5: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::Diff restrict l3: 0.00000
+DEAL::Diff restrict add l3: 0.00000
+DEAL::Diff restrict l4: 0.00000
+DEAL::Diff restrict add l4: 0.00000
+DEAL::Diff restrict l5: 0.00000
+DEAL::Diff restrict add l5: 0.00000
+DEAL::
+DEAL::FE: FE_Q<2>(4)
+DEAL::FE: FE_Q<2>(3)
+DEAL::no. cells: 34
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff prolongate l3: 0.00000
+DEAL::Diff prolongate l4: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::Diff restrict l3: 0.00000
+DEAL::Diff restrict add l3: 0.00000
+DEAL::Diff restrict l4: 0.00000
+DEAL::Diff restrict add l4: 0.00000
+DEAL::
+DEAL::no. cells: 88
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff prolongate l3: 0.00000
+DEAL::Diff prolongate l4: 0.00000
+DEAL::Diff prolongate l5: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::Diff restrict l3: 0.00000
+DEAL::Diff restrict add l3: 0.00000
+DEAL::Diff restrict l4: 0.00000
+DEAL::Diff restrict add l4: 0.00000
+DEAL::Diff restrict l5: 0.00000
+DEAL::Diff restrict add l5: 0.00000
+DEAL::
+DEAL::no. cells: 141
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff prolongate l3: 0.00000
+DEAL::Diff prolongate l4: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::Diff restrict l3: 0.00000
+DEAL::Diff restrict add l3: 0.00000
+DEAL::Diff restrict l4: 0.00000
+DEAL::Diff restrict add l4: 0.00000
+DEAL::
+DEAL::FE: FE_Q<3>(2)
+DEAL::FE: FE_Q<3>(1)
+DEAL::no. cells: 15
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::
+DEAL::no. cells: 694
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff prolongate l3: 0.00000
+DEAL::Diff prolongate l4: 0.00000
+DEAL::Diff prolongate l5: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::Diff restrict l3: 0.00000
+DEAL::Diff restrict add l3: 0.00000
+DEAL::Diff restrict l4: 0.00000
+DEAL::Diff restrict add l4: 0.00000
+DEAL::Diff restrict l5: 0.00000
+DEAL::Diff restrict add l5: 0.00000
+DEAL::
+DEAL::no. cells: 1791
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff prolongate l3: 0.00000
+DEAL::Diff prolongate l4: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::Diff restrict l3: 0.00000
+DEAL::Diff restrict add l3: 0.00000
+DEAL::Diff restrict l4: 0.00000
+DEAL::Diff restrict add l4: 0.00000
+DEAL::
+DEAL::FE: FE_Q<3>(4)
+DEAL::FE: FE_Q<3>(3)
+DEAL::no. cells: 1
+DEAL::
+DEAL::no. cells: 15
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::
+DEAL::no. cells: 223
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff prolongate l3: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::Diff restrict l3: 0.00000
+DEAL::Diff restrict add l3: 0.00000
+DEAL::
+DEAL::FE: FE_Q<2>(3)
+DEAL::FE: FE_Q<2>(2)
+DEAL::no. cells: 88
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff prolongate l3: 0.00000
+DEAL::Diff prolongate l4: 0.00000
+DEAL::Diff prolongate l5: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::Diff restrict l3: 0.00000
+DEAL::Diff restrict add l3: 0.00000
+DEAL::Diff restrict l4: 0.00000
+DEAL::Diff restrict add l4: 0.00000
+DEAL::Diff restrict l5: 0.00000
+DEAL::Diff restrict add l5: 0.00000
+DEAL::
+DEAL::no. cells: 250
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff prolongate l3: 0.00000
+DEAL::Diff prolongate l4: 0.00000
+DEAL::Diff prolongate l5: 0.00000
+DEAL::Diff prolongate l6: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::Diff restrict l3: 0.00000
+DEAL::Diff restrict add l3: 0.00000
+DEAL::Diff restrict l4: 0.00000
+DEAL::Diff restrict add l4: 0.00000
+DEAL::Diff restrict l5: 0.00000
+DEAL::Diff restrict add l5: 0.00000
+DEAL::Diff restrict l6: 0.00000
+DEAL::Diff restrict add l6: 0.00000
+DEAL::
+DEAL::no. cells: 510
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff prolongate l3: 0.00000
+DEAL::Diff prolongate l4: 0.00000
+DEAL::Diff prolongate l5: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::Diff restrict l3: 0.00000
+DEAL::Diff restrict add l3: 0.00000
+DEAL::Diff restrict l4: 0.00000
+DEAL::Diff restrict add l4: 0.00000
+DEAL::Diff restrict l5: 0.00000
+DEAL::Diff restrict add l5: 0.00000
+DEAL::
+DEAL::FE: FE_Q<3>(3)
+DEAL::FE: FE_Q<3>(2)
+DEAL::no. cells: 15
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::
+DEAL::no. cells: 694
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff prolongate l3: 0.00000
+DEAL::Diff prolongate l4: 0.00000
+DEAL::Diff prolongate l5: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::Diff restrict l3: 0.00000
+DEAL::Diff restrict add l3: 0.00000
+DEAL::Diff restrict l4: 0.00000
+DEAL::Diff restrict add l4: 0.00000
+DEAL::Diff restrict l5: 0.00000
+DEAL::Diff restrict add l5: 0.00000
+DEAL::
+DEAL::no. cells: 1791
+DEAL::Diff prolongate l1: 0.00000
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff prolongate l3: 0.00000
+DEAL::Diff prolongate l4: 0.00000
+DEAL::Diff restrict l1: 0.00000
+DEAL::Diff restrict add l1: 0.00000
+DEAL::Diff restrict l2: 0.00000
+DEAL::Diff restrict add l2: 0.00000
+DEAL::Diff restrict l3: 0.00000
+DEAL::Diff restrict add l3: 0.00000
+DEAL::Diff restrict l4: 0.00000
+DEAL::Diff restrict add l4: 0.00000
+DEAL::