MatrixFreeTools::compute_matrix() for PETScWrappers::MPI::SparseMatrix

diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h
index d3498fba298e215e6d512ecb2d87dbec813b93d8..7040ffffd83eca50ea22d266b357d19b912aec4c 100644 (file)
--- a/include/deal.II/matrix_free/tools.h
+++ b/include/deal.II/matrix_free/tools.h
@@ -826,7 +826,8 @@ namespace MatrixFreeTools
 
             for (unsigned int v = 0; v < n_filled_lanes; ++v)
               {
-                const auto cell_v = matrix_free.get_cell_iterator(cell, v);
+                const auto cell_v =
+                  matrix_free.get_cell_iterator(cell, v, dof_no);
 
                 if (matrix_free.get_mg_level() != numbers::invalid_unsigned_int)
                   cell_v->get_mg_dof_indices(dof_indices);
@@ -837,7 +838,6 @@ namespace MatrixFreeTools
                   dof_indices_mf[j] = dof_indices[lexicographic_numbering[j]];
 
                 constraints.distribute_local_to_global(matrices[v],
-                                                       dof_indices_mf,
                                                        dof_indices_mf,
                                                        dst);
               }
@@ -847,13 +847,6 @@ namespace MatrixFreeTools
       matrix);
 
     matrix.compress(VectorOperation::add);
-
-    for (auto &entry : matrix)
-      if (entry.row() == entry.column() && entry.value() == 0.0)
-        {
-          Assert(constraints.is_constrained(entry.row()), ExcNotImplemented());
-          entry.value() = 1.0;
-        }
   }
 
   template <typename CLASS,

diff --git a/tests/matrix_free/compute_diagonal_util.h b/tests/matrix_free/compute_diagonal_util.h
index 6d2a3c5ddda641a034667e40ac22c24a7978a620..6486dc4dc13fff44752702d645b612d9337b0391 100644 (file)
--- a/tests/matrix_free/compute_diagonal_util.h
+++ b/tests/matrix_free/compute_diagonal_util.h
@@ -188,6 +188,10 @@ public:
 
         for (; a1 != A1.end(); ++a1, ++a2, ++a_ref)
           {
+            if (a1->row() == a1->column() &&
+                constraints.is_constrained(a1->row()))
+              continue;
+
             Assert(std::abs(a1->value() - a_ref->value()) < 1e-6,
                    ExcNotImplemented());
             Assert(std::abs(a2->value() - a_ref->value()) < 1e-6,