Make sure to clear all variables in MGTransferMatrixFree.

diff --git a/doc/news/changes/minor/20170222MartinKronbichler b/doc/news/changes/minor/20170222MartinKronbichler
new file mode 100644 (file)
index 0000000..6118425
--- /dev/null
+++ b/doc/news/changes/minor/20170222MartinKronbichler
@@ -0,0 +1,5 @@
+Fixed: The class MGTransferMatrixFree::clear() did not reset all member
+variables, sometimes treating Dirichlet boundary conditions incorrectly if a
+transfer object was reused. This is now fixed.
+<br>
+(Martin Kronbichler, 2017/02/22)

diff --git a/source/multigrid/mg_transfer_matrix_free.cc b/source/multigrid/mg_transfer_matrix_free.cc
index b2a7e720059d2e6e4f6eaefeada4d7b5cf4c7c0b..17e7b9f43492f2997c0f8704cae92f5ab55dfb95 100644 (file)
--- a/source/multigrid/mg_transfer_matrix_free.cc
+++ b/source/multigrid/mg_transfer_matrix_free.cc
@@ -83,6 +83,7 @@ void MGTransferMatrixFree<dim,Number>::clear ()
   n_child_cell_dofs = 0;
   level_dof_indices.clear();
   parent_child_connect.clear();
+  dirichlet_indices.clear();
   n_owned_level_cells.clear();
   prolongation_matrix_1d.clear();
   evaluation_data.clear();
@@ -101,7 +102,7 @@ void MGTransferMatrixFree<dim,Number>::build
   internal::MGTransfer::ElementInfo<Number> elem_info;
 
   internal::MGTransfer::setup_transfer<dim,Number>(mg_dof,
-                                                   (const MGConstrainedDoFs *)this->mg_constrained_dofs,
+                                                   this->mg_constrained_dofs,
                                                    elem_info,
                                                    level_dof_indices,
                                                    parent_child_connect,

diff --git a/tests/multigrid/transfer_matrix_free_09.cc b/tests/multigrid/transfer_matrix_free_09.cc
new file mode 100644 (file)
index 0000000..d5816e5
--- /dev/null
+++ b/tests/multigrid/transfer_matrix_free_09.cc
@@ -0,0 +1,114 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Check that MGTransferMatrixFree::clear works correctly by comparing a
+// cleared transfer system with a freshly constructed one
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/multigrid/mg_transfer.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+#include <deal.II/numerics/vector_tools.h>
+
+
+template <int dim, typename Number>
+void check(const FiniteElement<dim> &fe)
+{
+  deallog.threshold_double(std::max(1e2*(double)std::numeric_limits<Number>::epsilon(),
+                                    1e-11));
+  deallog << "FE: " << fe.get_name() << std::endl;
+
+  MGConstrainedDoFs mg_constrained_dofs;
+  MGTransferMatrixFree<dim,Number> transfer;
+  Triangulation<dim> tr(Triangulation<dim>::limit_level_difference_at_vertices);
+  GridGenerator::hyper_cube(tr);
+  tr.refine_global(6-dim);
+
+  // run a few different sizes...
+  for (unsigned int c=0; c<4; ++c)
+    {
+      for (typename Triangulation<dim>::active_cell_iterator cell=tr.begin_active();
+           cell != tr.end(); ++cell)
+        if ((cell->center().norm() < 0.5 && (cell->level() < 5 ||
+                                             cell->center().norm() > 0.45)
+             ||
+             (dim == 2 && cell->center().norm() > 1.2)))
+          cell->set_refine_flag();
+      tr.execute_coarsening_and_refinement ();
+
+      DoFHandler<dim> mgdof(tr);
+      mgdof.distribute_dofs(fe);
+      mgdof.distribute_mg_dofs(fe);
+
+      Tensor<1,dim> exponents_monomial;
+      for (unsigned int d=0; d<dim; ++d)
+        exponents_monomial[d] = 1;
+      LinearAlgebra::distributed::Vector<double> vref;
+      vref.reinit(mgdof.n_dofs());
+      VectorTools::interpolate(mgdof, Functions::Monomial<dim>(exponents_monomial), vref);
+
+      deallog << "no. cells: " << tr.n_global_active_cells() << std::endl;
+
+      std::set<types::boundary_id> dirichlet_boundary;
+      dirichlet_boundary.insert(0);
+      mg_constrained_dofs.initialize(mgdof);
+      mg_constrained_dofs.make_zero_boundary_constraints(mgdof, dirichlet_boundary);
+
+      // build matrix-free transfer
+      transfer.clear();
+      transfer.initialize_constraints(mg_constrained_dofs);
+      transfer.build(mgdof);
+      MGTransferMatrixFree<dim, Number> transfer_ref(mg_constrained_dofs);
+      transfer_ref.build(mgdof);
+      MGLevelObject<LinearAlgebra::distributed::Vector<Number> > vectors(0, tr.n_global_levels()-1);
+      transfer_ref.copy_to_mg(mgdof, vectors, vref);
+      for (unsigned int level=vectors.max_level(); level>0; --level)
+        {
+          LinearAlgebra::distributed::Vector<Number> vec2(vectors[level-1]);
+          transfer_ref.restrict_and_add(level, vectors[level-1], vectors[level]);
+          transfer.restrict_and_add(level, vec2, vectors[level]);
+          vec2 -= vectors[level-1];
+          deallog << "Error in restriction:  " << (double)vec2.linfty_norm() << std::endl;
+        }
+
+      for (unsigned int level=1; level<vectors.max_level(); ++level)
+        {
+          LinearAlgebra::distributed::Vector<Number> vec2(vectors[level+1]);
+          transfer_ref.prolongate(level+1, vectors[level+1], vectors[level]);
+          transfer.prolongate(level+1, vec2, vectors[level]);
+          vec2 -= vectors[level+1];
+          deallog << "Error in prolongation: " << (double)vec2.linfty_norm() << 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>(FE_DGQ<2>(2));
+  check<2,double>(FE_Q<2>(2));
+  check<3,double>(FE_Q<3>(1));
+  check<2,float>(FE_Q<2>(2));
+}

diff --git a/tests/multigrid/transfer_matrix_free_09.output b/tests/multigrid/transfer_matrix_free_09.output
new file mode 100644 (file)
index 0000000..37eaf90
--- /dev/null
+++ b/tests/multigrid/transfer_matrix_free_09.output
@@ -0,0 +1,205 @@
+
+DEAL::FE: FE_DGQ<2>(2)
+DEAL::no. cells: 448
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::no. cells: 766
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::no. cells: 1924
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::no. cells: 6082
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::FE: FE_Q<2>(2)
+DEAL::no. cells: 448
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::no. cells: 766
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::no. cells: 1924
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::no. cells: 6082
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::FE: FE_Q<3>(1)
+DEAL::no. cells: 757
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::no. cells: 2808
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::no. cells: 6847
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::no. cells: 33776
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::FE: FE_Q<2>(2)
+DEAL::no. cells: 448
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::no. cells: 766
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::no. cells: 1924
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::no. cells: 6082
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in restriction:  0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0
+DEAL::Error in prolongation: 0