]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add tests
authorDaniel Arndt <d.arndt@math.uni-goettingen.de>
Sun, 1 Oct 2017 18:14:50 +0000 (20:14 +0200)
committerDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Mon, 2 Oct 2017 08:17:34 +0000 (10:17 +0200)
tests/matrix_free/parallel_multigrid_adaptive_07.cc [new file with mode: 0644]
tests/matrix_free/parallel_multigrid_adaptive_07.with_mpi=true.with_p4est=true.mpirun=3.output [new file with mode: 0644]
tests/multigrid/transfer_06.cc [new file with mode: 0644]
tests/multigrid/transfer_06.with_mpi=true.with_p4est=true.mpirun=3.output [new file with mode: 0644]
tests/multigrid/transfer_matrix_free_02_block_02.cc [new file with mode: 0644]
tests/multigrid/transfer_matrix_free_02_block_02.with_mpi=true.with_p4est=true.with_trilinos=true.mpirun=5.output [new file with mode: 0644]

diff --git a/tests/matrix_free/parallel_multigrid_adaptive_07.cc b/tests/matrix_free/parallel_multigrid_adaptive_07.cc
new file mode 100644 (file)
index 0000000..b737326
--- /dev/null
@@ -0,0 +1,419 @@
+// ---------------------------------------------------------------------
+//
+// 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();
+  }
+
+}
diff --git a/tests/matrix_free/parallel_multigrid_adaptive_07.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/matrix_free/parallel_multigrid_adaptive_07.with_mpi=true.with_p4est=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..158a6bc
--- /dev/null
@@ -0,0 +1,108 @@
+
+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
diff --git a/tests/multigrid/transfer_06.cc b/tests/multigrid/transfer_06.cc
new file mode 100644 (file)
index 0000000..751b5b0
--- /dev/null
@@ -0,0 +1,218 @@
+// ---------------------------------------------------------------------
+//
+// 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);
+}
diff --git a/tests/multigrid/transfer_06.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/multigrid/transfer_06.with_mpi=true.with_p4est=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..0a70122
--- /dev/null
@@ -0,0 +1,265 @@
+
+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::
diff --git a/tests/multigrid/transfer_matrix_free_02_block_02.cc b/tests/multigrid/transfer_matrix_free_02_block_02.cc
new file mode 100644 (file)
index 0000000..6ab03ee
--- /dev/null
@@ -0,0 +1,181 @@
+// ---------------------------------------------------------------------
+//
+// 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);
+}
diff --git a/tests/multigrid/transfer_matrix_free_02_block_02.with_mpi=true.with_p4est=true.with_trilinos=true.mpirun=5.output b/tests/multigrid/transfer_matrix_free_02_block_02.with_mpi=true.with_p4est=true.with_trilinos=true.mpirun=5.output
new file mode 100644 (file)
index 0000000..62ed8e7
--- /dev/null
@@ -0,0 +1,265 @@
+
+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::

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.