Fix bug in distribute_local_to_global

If a local matrix with all diagonal elements equal to zero is distributed to a global matrix, the l1 norm of the local matrix divided by the size of the local matrix is now added to those diagonal elements of the global matrix which correspond to a constrained dof. In case the entire local matrix is zero, 1 is added. Previously zero was added for both cases, possibly resulting in singular global matrices. Additionally, a test for the patched version of distribute_local_to_global has been added.
Fixes #7658

diff --git a/doc/news/changes/minor/20190213SebastianStark b/doc/news/changes/minor/20190213SebastianStark
new file mode 100644 (file)
index 0000000..4471051
--- /dev/null
+++ b/doc/news/changes/minor/20190213SebastianStark
@@ -0,0 +1,3 @@
+Fixed: Handle the case properly that a local matrix with all diagonal elements equal to zero is distributed with AffineConstraints::distribute_local_to_global() to the global matrix while constraints apply.
+<br>
+(Sebastian Stark, 2019/02/13)

diff --git a/include/deal.II/lac/affine_constraints.templates.h b/include/deal.II/lac/affine_constraints.templates.h
index dc0df47d7f5a807374d93dfe560f4166e905f4b1..9eb936ec9f73ad81bb0b2e48dc35784b8e8ff3d4 100644 (file)
--- a/include/deal.II/lac/affine_constraints.templates.h
+++ b/include/deal.II/lac/affine_constraints.templates.h
@@ -3325,9 +3325,13 @@ namespace internals
   }
 
   // to make sure that the global matrix remains invertible, we need to do
-  // something with the diagonal elements. add the absolute value of the local
-  // matrix, so the resulting entry will always be positive and furthermore be
-  // in the same order of magnitude as the other elements of the matrix
+  // something with the diagonal elements. Add the average of the
+  // absolute values of the local matrix diagonals, so the resulting entry
+  // will always be positive and furthermore be in the same order of magnitude
+  // as the other elements of the matrix. If all local matrix diagonals are
+  // zero, add the l1 norm of the local matrix divided by the matrix size
+  // to the diagonal of the global matrix. If the entire local matrix is zero,
+  // add 1 to the diagonal of the global matrix.
   //
   // note that this also captures the special case that a dof is both
   // constrained and fixed (this can happen for hanging nodes in 3d that also
@@ -3355,6 +3359,16 @@ namespace internals
           average_diagonal += std::abs(local_matrix(i, i));
         average_diagonal /= static_cast<number>(local_matrix.m());
 
+        // handle the case that all diagonal elements are zero
+        if (average_diagonal == static_cast<number>(0.))
+          {
+            average_diagonal = static_cast<number>(local_matrix.l1_norm()) /
+                               static_cast<number>(local_matrix.m());
+            // if the entire matrix is zero, use 1. for the diagonal
+            if (average_diagonal == static_cast<number>(0.))
+              average_diagonal = static_cast<number>(1.);
+          }
+
         for (size_type i = 0; i < global_rows.n_constraints(); i++)
           {
             const size_type local_row  = global_rows.constraint_origin(i);

diff --git a/tests/lac/constraint_matrix_distribute_local_global.cc b/tests/lac/constraint_matrix_distribute_local_global.cc
new file mode 100644 (file)
index 0000000..c76d137
--- /dev/null
+++ b/tests/lac/constraint_matrix_distribute_local_global.cc
@@ -0,0 +1,78 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2005 - 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <stdio.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+// Tests distribute_local_to_global for
+// (1) the case that a local matrix
+//     with all diagonals equal to zero is
+//     distributed while a dof is
+//     constrained, and
+// (2) the case that all entries of the local matrix
+//     are zero while a dof is constrained.
+
+using namespace dealii;
+
+int
+main()
+{
+  initlog();
+
+  // set up constraint
+  AffineConstraints<double> constraints;
+  constraints.add_line(0);
+  constraints.close();
+
+  // global matrix and vector
+  FullMatrix<double> global_matrix(2);
+  Vector<double>     global_vector(2);
+
+  // first test: add matrix with all diagonals zero
+  FullMatrix<double> local_matrix(2);
+  local_matrix(1, 0) = local_matrix(0, 1) = 1.;
+  Vector<double> local_vector(2);
+  constraints.distribute_local_to_global(
+    local_matrix, local_vector, {0, 1}, global_matrix, global_vector, true);
+  // output result
+  for (unsigned int m = 0; m < global_matrix.m(); ++m)
+    {
+      for (unsigned int n = 0; n < global_matrix.n(); ++n)
+        deallog << global_matrix(m, n) << " ";
+      deallog << std::endl;
+    }
+  deallog << std::endl;
+
+  // second test: add matrix with all entries zero
+  global_matrix = 0.;
+  local_matrix  = 0.;
+  constraints.distribute_local_to_global(
+    local_matrix, local_vector, {0, 1}, global_matrix, global_vector, true);
+  // output result
+  for (unsigned int m = 0; m < global_matrix.m(); ++m)
+    {
+      for (unsigned int n = 0; n < global_matrix.n(); ++n)
+        deallog << global_matrix(m, n) << " ";
+      deallog << std::endl;
+    }
+  deallog << std::endl;
+}

diff --git a/tests/lac/constraint_matrix_distribute_local_global.output b/tests/lac/constraint_matrix_distribute_local_global.output
new file mode 100644 (file)
index 0000000..6e129c4
--- /dev/null
+++ b/tests/lac/constraint_matrix_distribute_local_global.output
@@ -0,0 +1,7 @@
+
+DEAL::0.500000 0.00000 
+DEAL::0.00000 0.00000 
+DEAL::
+DEAL::1.00000 0.00000 
+DEAL::0.00000 0.00000 
+DEAL::