From 7af998f7d482c5f6994def49951900c6306cf0b8 Mon Sep 17 00:00:00 2001
From: David Wells <drwells@email.unc.edu>
Date: Thu, 29 Jul 2021 17:40:18 -0400
Subject: [PATCH] Make MassOperator's diagonals true diagonals, not lumped
 ones.

The previous implementation used row sums (i.e., mass lumping) instead of
computing the diagonal. This should be changed for two reasons:
1. It's not really the diagonal (and its inconsistent with LaplaceOperator).
2. Lumped mass matrices are only positive for either low-order elements (where
   the basis functions are all positive) or very well-behaved elements (like
   FE_Q). In particular, the values on the diagonal are either zero or negative
   for TET10 if we use lumping, which isn't going to work.
---
 .../incompatibilities/20210729DavidWells      |  4 ++
 include/deal.II/matrix_free/operators.h       | 57 ++++++++++++++++---
 2 files changed, 53 insertions(+), 8 deletions(-)
 create mode 100644 doc/news/changes/incompatibilities/20210729DavidWells

diff --git a/doc/news/changes/incompatibilities/20210729DavidWells b/doc/news/changes/incompatibilities/20210729DavidWells
new file mode 100644
index 0000000000..e4ddc17dfe
--- /dev/null
+++ b/doc/news/changes/incompatibilities/20210729DavidWells
@@ -0,0 +1,4 @@
+Changed: MatrixFreeOperators::MassOperator now computes its true diagonal (and
+not a lumped mass matrix) as its diagonal.
+<br>
+(David Wells, 2021/07/29)
diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h
index 8465ea3415..3608a37ffe 100644
--- a/include/deal.II/matrix_free/operators.h
+++ b/include/deal.II/matrix_free/operators.h
@@ -29,6 +29,7 @@
 
 #include <deal.II/matrix_free/fe_evaluation.h>
 #include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/tools.h>
 
 #include <deal.II/multigrid/mg_constrained_dofs.h>
 
@@ -40,7 +41,7 @@ namespace MatrixFreeOperators
 {
   namespace BlockHelper
   {
-    // workaroud for unifying non-block vector and block vector implementations
+    // workaround for unifying non-block vector and block vector implementations
     // a non-block vector has one block and the only subblock is the vector
     // itself
     template <typename VectorType>
@@ -357,6 +358,11 @@ namespace MatrixFreeOperators
      * A derived class needs to implement this function and resize and fill
      * the protected member inverse_diagonal_entries and/or diagonal_entries
      * accordingly.
+     *
+     * @note Since the diagonal is frequently used as a smoother or
+     * preconditioner, entries corresponding to constrained DoFs are set to 1
+     * (instead of the correct value of 0) so that the diagonal matrix is
+     * invertible.
      */
     virtual void
     compute_diagonal() = 0;
@@ -1793,6 +1799,9 @@ namespace MatrixFreeOperators
       typename Base<dim, VectorType, VectorizedArrayType>::value_type;
     Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
            ExcNotInitialized());
+    Assert(this->selected_rows == this->selected_columns,
+           ExcMessage("This function is only implemented for square (not "
+                      "rectangular) operators."));
 
     this->inverse_diagonal_entries =
       std::make_shared<DiagonalMatrix<VectorType>>();
@@ -1802,17 +1811,49 @@ namespace MatrixFreeOperators
     VectorType &diagonal_vector = this->diagonal_entries->get_vector();
     this->initialize_dof_vector(inverse_diagonal_vector);
     this->initialize_dof_vector(diagonal_vector);
-    inverse_diagonal_vector = Number(1.);
-    apply_add(diagonal_vector, inverse_diagonal_vector);
 
+    // Set up the action of the mass matrix in a way that's compatible with
+    // MatrixFreeTools::compute_diagonal:
+    auto diagonal_evaluation = [](auto &integrator) {
+      integrator.evaluate(EvaluationFlags::values);
+      for (unsigned int q = 0; q < integrator.n_q_points; ++q)
+        integrator.submit_value(integrator.get_value(q), q);
+      integrator.integrate(EvaluationFlags::values);
+    };
+
+    std::function<void(
+      FEEvaluation<
+        dim,
+        fe_degree,
+        n_q_points_1d,
+        n_components,
+        typename Base<dim, VectorType, VectorizedArrayType>::value_type,
+        VectorizedArrayType> &)>
+      diagonal_evaluation_f(diagonal_evaluation);
+
+    Assert(this->selected_rows.size() > 0, ExcInternalError());
+    for (unsigned int block_n = 0; block_n < this->selected_rows.size();
+         ++block_n)
+      MatrixFreeTools::compute_diagonal(*this->data,
+                                        BlockHelper::subblock(diagonal_vector,
+                                                              block_n),
+                                        diagonal_evaluation_f,
+                                        this->selected_rows[block_n]);
+
+    // Constrained entries will create zeros on the main diagonal, which we
+    // don't want
     this->set_constrained_entries_to_one(diagonal_vector);
+
     inverse_diagonal_vector = diagonal_vector;
 
-    const unsigned int locally_owned_size =
-      inverse_diagonal_vector.locally_owned_size();
-    for (unsigned int i = 0; i < locally_owned_size; ++i)
-      inverse_diagonal_vector.local_element(i) =
-        Number(1.) / inverse_diagonal_vector.local_element(i);
+    for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
+         ++i)
+      {
+        Assert(diagonal_vector.local_element(i) > Number(0),
+               ExcInternalError());
+        inverse_diagonal_vector.local_element(i) =
+          1. / inverse_diagonal_vector.local_element(i);
+      }
 
     inverse_diagonal_vector.update_ghost_values();
     diagonal_vector.update_ghost_values();
-- 
2.39.5