MatrixFreeTools::compute_matrix(): specialize assembly of diagonal blocks

diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h
index 92e3f83166ea88537f3292653e6c56f80f6a83e0..1afa1983af997b224c56fa12f7230f3027a29e8c 100644 (file)
--- a/include/deal.II/matrix_free/tools.h
+++ b/include/deal.II/matrix_free/tools.h
@@ -2221,11 +2221,20 @@ namespace MatrixFreeTools
 
                   for (unsigned int v = 0; v < n_filled_lanes; ++v)
                     for (unsigned int bi = 0; bi < n_blocks; ++bi)
-                      constraints.distribute_local_to_global(
-                        matrices[bi][bj][v],
-                        dof_indices_mf[bi][v],
-                        dof_indices_mf[bj][v],
-                        matrix);
+                      if (bi == bj)
+                        // specialization for blocks on the diagonal
+                        // to writing into diagonal elements of the
+                        // matrix if the corresponding degree of freedom
+                        // is constrained, see also the documentation
+                        // of AffineConstraints::distribute_local_to_global()
+                        constraints.distribute_local_to_global(
+                          matrices[bi][bi][v], dof_indices_mf[bi][v], matrix);
+                      else
+                        constraints.distribute_local_to_global(
+                          matrices[bi][bj][v],
+                          dof_indices_mf[bi][v],
+                          dof_indices_mf[bj][v],
+                          matrix);
                 }
             }
         };