Fix zeroing vector without any non-constrained dofs.

diff --git a/include/deal.II/matrix_free/dof_info.templates.h b/include/deal.II/matrix_free/dof_info.templates.h
index 58ceba98791a0fa94990d20f189e1b22eb84b315..b46c0e91e4390f0bec5dff409374bf1606a2e507 100644 (file)
--- a/include/deal.II/matrix_free/dof_info.templates.h
+++ b/include/deal.II/matrix_free/dof_info.templates.h
@@ -1197,6 +1197,10 @@ namespace internal
                       }
                   }
           }
+      // ensure that all indices are touched at least during the last round
+      for (auto &i : touched_by)
+        if (i == numbers::invalid_unsigned_int)
+          i = task_info.cell_partition_data.back() - 1;
 
       vector_zero_range_list_index.resize(
         1 + task_info

diff --git a/tests/matrix_free/matrix_vector_24.cc b/tests/matrix_free/matrix_vector_24.cc
new file mode 100644 (file)
index 0000000..44275e5
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_24.cc
@@ -0,0 +1,194 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2013 - 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// this function tests the correctness of the implementation of matrix free
+// matrix-vector products for a special case of linear elements where all DoFs
+// are subject to constraints
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+#include "matrix_vector_mf.h"
+
+
+
+template <int dim,
+          int fe_degree,
+          typename Number,
+          typename VectorType = Vector<Number>,
+          int n_q_points_1d   = fe_degree + 1>
+class MatrixFreeVariant
+{
+public:
+  typedef VectorizedArray<Number> vector_t;
+
+  MatrixFreeVariant(const MatrixFree<dim, Number> &data_in)
+    : data(data_in)
+  {}
+
+  void
+  vmult(VectorType &dst, const VectorType &src) const
+  {
+    const std::function<
+      void(const MatrixFree<dim, typename VectorType::value_type> &,
+           VectorType &,
+           const VectorType &,
+           const std::pair<unsigned int, unsigned int> &)>
+      wrap = helmholtz_operator<dim, fe_degree, VectorType, n_q_points_1d>;
+    data.cell_loop(wrap, dst, src, true);
+    for (auto i : data.get_constrained_dofs())
+      dst(i) += src(i);
+  }
+
+private:
+  const MatrixFree<dim, Number> &data;
+};
+
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+
+  FE_Q<dim>       fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+  AffineConstraints<double> constraints;
+  VectorTools::interpolate_boundary_values(dof,
+                                           0,
+                                           ZeroFunction<dim>(),
+                                           constraints);
+  constraints.close();
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+  using number = double;
+
+  MatrixFree<dim, number> mf_data;
+  {
+    const QGauss<1>                                  quad(fe_degree + 1);
+    typename MatrixFree<dim, number>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+    mf_data.reinit(dof, constraints, quad, data);
+  }
+
+  MatrixFreeVariant<dim,
+                    fe_degree,
+                    number,
+                    LinearAlgebra::distributed::Vector<number>,
+                    fe_degree + 1>
+                                             mf(mf_data);
+  LinearAlgebra::distributed::Vector<number> in(dof.n_dofs()),
+    out(dof.n_dofs());
+  LinearAlgebra::distributed::Vector<number> out_dist(in);
+
+  for (unsigned int i = 0; i < dof.n_dofs(); ++i)
+    {
+      const double entry = random_value<double>();
+      in(i)              = entry;
+    }
+
+  out = in;
+  mf.vmult(out, in);
+
+  // assemble sparse matrix with (\nabla v, \nabla u) + (v, 10 * u)
+  SparsityPattern sparsity;
+  {
+    DynamicSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern(dof, csp, constraints, true);
+    sparsity.copy_from(csp);
+  }
+  SparseMatrix<double> sparse_matrix(sparsity);
+  {
+    QGauss<dim> quadrature_formula(fe_degree + 1);
+
+    FEValues<dim> fe_values(dof.get_fe(),
+                            quadrature_formula,
+                            update_values | update_gradients |
+                              update_JxW_values);
+
+    const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
+    const unsigned int n_q_points    = quadrature_formula.size();
+
+    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
+                                                   endc = dof.end();
+    for (; cell != endc; ++cell)
+      {
+        cell_matrix = 0;
+        fe_values.reinit(cell);
+
+        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+          for (unsigned int i = 0; i < dofs_per_cell; ++i)
+            {
+              for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                cell_matrix(i, j) += ((fe_values.shape_grad(i, q_point) *
+                                         fe_values.shape_grad(j, q_point) +
+                                       10. * fe_values.shape_value(i, q_point) *
+                                         fe_values.shape_value(j, q_point)) *
+                                      fe_values.JxW(q_point));
+            }
+
+        cell->get_dof_indices(local_dof_indices);
+        constraints.distribute_local_to_global(cell_matrix,
+                                               local_dof_indices,
+                                               sparse_matrix);
+      }
+  }
+  // set matrix entries to constrained rows to 1 for consistency with
+  // matrix-free
+  for (unsigned int i = 0; i < dof.n_dofs(); ++i)
+    if (constraints.is_constrained(i))
+      sparse_matrix.set(i, i, 1);
+
+  sparse_matrix.vmult(out_dist, in);
+  out -= out_dist;
+  const double diff_norm = out.linfty_norm() / out_dist.linfty_norm();
+
+  deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(3);
+
+  test<2, 1>();
+  test<3, 1>();
+}

diff --git a/tests/matrix_free/matrix_vector_24.output b/tests/matrix_free/matrix_vector_24.output
new file mode 100644 (file)
index 0000000..92d8da9
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_24.output
@@ -0,0 +1,7 @@
+
+DEAL::Testing FE_Q<2>(1)
+DEAL::Norm of difference: 0.00
+DEAL::
+DEAL::Testing FE_Q<3>(1)
+DEAL::Norm of difference: 0.00
+DEAL::