From f1e01d0845b0ed1257b5b4aca88503470a45e570 Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 7 Jun 2004 16:06:56 +0000 Subject: [PATCH] Improve the computation of new diagonal entries when eliminating a constrained DoF. git-svn-id: https://svn.dealii.org/trunk@9395 0785d39b-7218-0410-832d-ea1e28bc413d --- .../numerics/matrices.all_dimensions.cc | 51 ++++++++++++++++++- 1 file changed, 49 insertions(+), 2 deletions(-) diff --git a/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc b/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc index 01f141fadf..ef786a0040 100644 --- a/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc +++ b/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc @@ -669,6 +669,27 @@ local_apply_boundary_values (const std::map &boundary_value // otherwise traverse all the dofs used in // the local matrices and vectors and see // what's there to do + + // if we need to treat an entry, then we + // set the diagonal entry to its absolute + // value. if it is zero, we used to set it + // to one, which is a really terrible + // choice that can lead to hours of + // searching for bugs in programs (I + // experienced this :-( ) if the matrix + // entries are otherwise very large. this + // is so since iterative solvers would + // simply not correct boundary nodes for + // their correct values since the residual + // contributions of their rows of the + // linear system is almost zero if the + // diagonal entry is one. thus, set it to + // the average absolute value of the + // nonzero diagonal elements. + // + // we only compute this value lazily the + // first time we need it. + double average_diagonal = 0; const unsigned int n_local_dofs = local_dof_indices.size(); for (unsigned int i=0; i &boundary_value // replace diagonal entry by its // absolute value to make sure that // everything remains positive, or - // by one if zero + // by the average diagonal value if + // zero if (local_matrix(i,i) == 0.) - local_matrix(i,i) = 1.; + { + // if average diagonal hasn't + // yet been computed, do so now + if (average_diagonal == 0.) + { + unsigned int nonzero_diagonals = 0; + for (unsigned int k=0; k