From 79a1be0b7e1c67e844e4cd44a9e777efa1ec8887 Mon Sep 17 00:00:00 2001
From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Tue, 8 Jan 2013 21:19:57 +0000
Subject: [PATCH] Replace two untested and undocumented functions by a call to
 MatrixTools::apply_boundary_values.

git-svn-id: https://svn.dealii.org/trunk@27982 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/source/multigrid/mg_tools.cc | 405 +++------------------------
 1 file changed, 36 insertions(+), 369 deletions(-)

diff --git a/deal.II/source/multigrid/mg_tools.cc b/deal.II/source/multigrid/mg_tools.cc
index ce14f7e9fd..bd45b98853 100644
--- a/deal.II/source/multigrid/mg_tools.cc
+++ b/deal.II/source/multigrid/mg_tools.cc
@@ -34,6 +34,7 @@
 #include <deal.II/dofs/dof_tools.h>
 #include <deal.II/fe/fe.h>
 #include <deal.II/fe/component_mask.h>
+#include <deal.II/numerics/matrix_tools.h>
 
 #include <vector>
 #include <algorithm>
@@ -1625,131 +1626,23 @@ namespace MGTools
     const std::set<unsigned int> &boundary_dofs,
     SparseMatrix<number> &matrix,
     const bool preserve_symmetry,
-    const bool ignore_zeros)
+    const bool /*ignore_zeros*/)
   {
-    // if no boundary values are to be applied
-    // simply return
-    if (boundary_dofs.size() == 0)
-      return;
-
-
-    const unsigned int n_dofs = matrix.m();
-
-    // if a diagonal entry is zero
-    // later, then we use another
-    // number instead. take it to be
-    // the first nonzero diagonal
-    // element of the matrix, or 1 if
-    // there is no such thing
-    number first_nonzero_diagonal_entry = 1;
-    for (unsigned int i=0; i<n_dofs; ++i)
-      if (matrix.diag_element(i) != 0)
-        {
-          first_nonzero_diagonal_entry = matrix.diag_element(i);
-          break;
-        }
-
-
-    std::set<unsigned int>::const_iterator dof  = boundary_dofs.begin(),
-                                           endd = boundary_dofs.end();
-    const SparsityPattern    &sparsity    = matrix.get_sparsity_pattern();
-    const std::size_t  *sparsity_rowstart = sparsity.get_rowstart_indices();
-    const unsigned int *sparsity_colnums  = sparsity.get_column_numbers();
-    for (; dof != endd; ++dof)
-      {
-        Assert (*dof < n_dofs, ExcInternalError());
-
-        const unsigned int dof_number = *dof;
-        // for each boundary dof:
-
-        // set entries of this line
-        // to zero except for the diagonal
-        // entry. Note that the diagonal
-        // entry is always the first one
-        // for square matrices, i.e.
-        // we shall not set
-        // matrix.global_entry(
-        //     sparsity_rowstart[dof.first])
-        const unsigned int last = sparsity_rowstart[dof_number+1];
-        for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
-          matrix.global_entry(j) = 0.;
-
-
-        // set right hand side to
-        // wanted value: if main diagonal
-        // entry nonzero, don't touch it
-        // and scale rhs accordingly. If
-        // zero, take the first main
-        // diagonal entry we can find, or
-        // one if no nonzero main diagonal
-        // element exists. Normally, however,
-        // the main diagonal entry should
-        // not be zero.
-        //
-        // store the new rhs entry to make
-        // the gauss step more efficient
-        if (!ignore_zeros)
-          matrix.set (dof_number, dof_number,
-                      first_nonzero_diagonal_entry);
-        // if the user wants to have
-        // the symmetry of the matrix
-        // preserved, and if the
-        // sparsity pattern is
-        // symmetric, then do a Gauss
-        // elimination step with the
-        // present row
-        if (preserve_symmetry)
-          {
-            // we have to loop over all
-            // rows of the matrix which
-            // have a nonzero entry in
-            // the column which we work
-            // in presently. if the
-            // sparsity pattern is
-            // symmetric, then we can
-            // get the positions of
-            // these rows cheaply by
-            // looking at the nonzero
-            // column numbers of the
-            // present row. we need not
-            // look at the first entry,
-            // since that is the
-            // diagonal element and
-            // thus the present row
-            for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
-              {
-                const unsigned int row = sparsity_colnums[j];
-
-                // find the position of
-                // element
-                // (row,dof_number)
-                const unsigned int *
-                p = Utilities::lower_bound(&sparsity_colnums[sparsity_rowstart[row]+1],
-                                           &sparsity_colnums[sparsity_rowstart[row+1]],
-                                           dof_number);
-
-                // check whether this line has
-                // an entry in the regarding column
-                // (check for ==dof_number and
-                // != next_row, since if
-                // row==dof_number-1, *p is a
-                // past-the-end pointer but points
-                // to dof_number anyway...)
-                //
-                // there should be such an entry!
-                Assert ((*p == dof_number) &&
-                        (p != &sparsity_colnums[sparsity_rowstart[row+1]]),
-                        ExcInternalError());
-
-                const unsigned int global_entry
-                  = (p - &sparsity_colnums[sparsity_rowstart[0]]);
-
-                // correct right hand side
-                // set matrix entry to zero
-                matrix.global_entry(global_entry) = 0.;
-              }
-          }
-      }
+    // this function is not documented and not tested in the testsuite
+    // so it isn't quite clear what it's supposed to do. it also isn't
+    // used anywhere else in the library. in avoiding the use of
+    // a deprecated function, I therefore threw away the original function
+    // and replaced it by the following, which I believe should work
+
+    std::map<unsigned int, double> boundary_values;
+    for (std::set<unsigned int>::const_iterator p=boundary_dofs.begin();
+        p != boundary_dofs.end(); ++p)
+      boundary_values[*p] = 0;
+
+    Vector<number> dummy(matrix.m());
+    MatrixTools::apply_boundary_values (boundary_values,
+                                        matrix, dummy, dummy,
+                                        preserve_symmetry);
   }
 
 
@@ -1773,251 +1666,25 @@ namespace MGTools
       Assert (matrix.block(i,i).get_sparsity_pattern().optimize_diagonal(),
               SparsityPattern::ExcDiagonalNotOptimized());
 
-
-    // if no boundary values are to be applied
-    // simply return
-    if (boundary_dofs.size() == 0)
-      return;
-
-
-    const unsigned int n_dofs = matrix.m();
-
-    // if a diagonal entry is zero
-    // later, then we use another
-    // number instead. take it to be
-    // the first nonzero diagonal
-    // element of the matrix, or 1 if
-    // there is no such thing
-    number first_nonzero_diagonal_entry = 0;
-    for (unsigned int diag_block=0; diag_block<blocks; ++diag_block)
-      {
-        for (unsigned int i=0; i<matrix.block(diag_block,diag_block).n(); ++i)
-          if (matrix.block(diag_block,diag_block).diag_element(i) != 0)
-            {
-              first_nonzero_diagonal_entry
-                = matrix.block(diag_block,diag_block).diag_element(i);
-              break;
-            }
-        // check whether we have found
-        // something in the present
-        // block
-        if (first_nonzero_diagonal_entry != 0)
-          break;
-      }
-    // nothing found on all diagonal
-    // blocks? if so, use 1.0 instead
-    if (first_nonzero_diagonal_entry == 0)
-      first_nonzero_diagonal_entry = 1;
-
-
-    std::set<unsigned int>::const_iterator dof  = boundary_dofs.begin(),
-                                           endd = boundary_dofs.end();
-    const BlockSparsityPattern &
-    sparsity_pattern = matrix.get_sparsity_pattern();
-
-    // pointer to the mapping between
-    // global and block indices. since
-    // the row and column mappings are
-    // equal, store a pointer on only
-    // one of them
-    const BlockIndices &
-    index_mapping = sparsity_pattern.get_column_indices();
-
-    // now loop over all boundary dofs
-    for (; dof != endd; ++dof)
-      {
-        Assert (*dof < n_dofs, ExcInternalError());
-
-        // get global index and index
-        // in the block in which this
-        // dof is located
-        const unsigned int dof_number = *dof;
-        const std::pair<unsigned int,unsigned int>
-        block_index = index_mapping.global_to_local (dof_number);
-
-        // for each boundary dof:
-
-        // set entries of this line
-        // to zero except for the diagonal
-        // entry. Note that the diagonal
-        // entry is always the first one
-        // for square matrices, i.e.
-        // we shall not set
-        // matrix.global_entry(
-        //     sparsity_rowstart[dof.first])
-        // of the diagonal block
-        for (unsigned int block_col=0; block_col<blocks; ++block_col)
-          {
-            const SparsityPattern &
-            local_sparsity = sparsity_pattern.block(block_index.first,
-                                                    block_col);
-
-            // find first and last
-            // entry in the present row
-            // of the present
-            // block. exclude the main
-            // diagonal element, which
-            // is the diagonal element
-            // of a diagonal block,
-            // which must be a square
-            // matrix so the diagonal
-            // element is the first of
-            // this row.
-            const unsigned int
-            last  = local_sparsity.get_rowstart_indices()[block_index.second+1],
-            first = (block_col == block_index.first ?
-                     local_sparsity.get_rowstart_indices()[block_index.second]+1 :
-                     local_sparsity.get_rowstart_indices()[block_index.second]);
-
-            for (unsigned int j=first; j<last; ++j)
-              matrix.block(block_index.first,block_col).global_entry(j) = 0.;
-          }
-
-        matrix.block(block_index.first, block_index.first)
-        .diag_element(block_index.second)
-          = first_nonzero_diagonal_entry;
-
-        // if the user wants to have
-        // the symmetry of the matrix
-        // preserved, and if the
-        // sparsity pattern is
-        // symmetric, then do a Gauss
-        // elimination step with the
-        // present row. this is a
-        // little more complicated for
-        // block matrices.
-        if (preserve_symmetry)
-          {
-            // we have to loop over all
-            // rows of the matrix which
-            // have a nonzero entry in
-            // the column which we work
-            // in presently. if the
-            // sparsity pattern is
-            // symmetric, then we can
-            // get the positions of
-            // these rows cheaply by
-            // looking at the nonzero
-            // column numbers of the
-            // present row.
-            //
-            // note that if we check
-            // whether row @p{row} in
-            // block (r,c) is non-zero,
-            // then we have to check
-            // for the existence of
-            // column @p{row} in block
-            // (c,r), i.e. of the
-            // transpose block
-            for (unsigned int block_row=0; block_row<blocks; ++block_row)
-              {
-                // get pointers to the
-                // sparsity patterns of
-                // this block and of
-                // the transpose one
-                const SparsityPattern &this_sparsity
-                  = sparsity_pattern.block (block_row, block_index.first);
-                const SparsityPattern &transpose_sparsity
-                  = sparsity_pattern.block (block_index.first, block_row);
-
-                // traverse the row of
-                // the transpose block
-                // to find the
-                // interesting rows in
-                // the present block.
-                // don't use the
-                // diagonal element of
-                // the diagonal block
-                const unsigned int
-                first = (block_index.first == block_row ?
-                         transpose_sparsity.get_rowstart_indices()[block_index.second]+1 :
-                         transpose_sparsity.get_rowstart_indices()[block_index.second]),
-                        last  = transpose_sparsity.get_rowstart_indices()[block_index.second+1];
-
-                for (unsigned int j=first; j<last; ++j)
-                  {
-                    // get the number
-                    // of the column in
-                    // this row in
-                    // which a nonzero
-                    // entry is. this
-                    // is also the row
-                    // of the transpose
-                    // block which has
-                    // an entry in the
-                    // interesting row
-                    const unsigned int row = transpose_sparsity.get_column_numbers()[j];
-
-                    // find the
-                    // position of
-                    // element
-                    // (row,dof_number)
-                    // in this block
-                    // (not in the
-                    // transpose
-                    // one). note that
-                    // we have to take
-                    // care of special
-                    // cases with
-                    // square
-                    // sub-matrices
-                    const unsigned int *p = 0;
-                    if (this_sparsity.n_rows() == this_sparsity.n_cols())
-                      {
-                        if (this_sparsity.get_column_numbers()
-                            [this_sparsity.get_rowstart_indices()[row]]
-                            ==
-                            block_index.second)
-                          p = &this_sparsity.get_column_numbers()
-                              [this_sparsity.get_rowstart_indices()[row]];
-                        else
-                          p = Utilities::lower_bound(&this_sparsity.get_column_numbers()
-                                                     [this_sparsity.get_rowstart_indices()[row]+1],
-                                                     &this_sparsity.get_column_numbers()
-                                                     [this_sparsity.get_rowstart_indices()[row+1]],
-                                                     block_index.second);
-                      }
-                    else
-                      p = Utilities::lower_bound(&this_sparsity.get_column_numbers()
-                                                 [this_sparsity.get_rowstart_indices()[row]],
-                                                 &this_sparsity.get_column_numbers()
-                                                 [this_sparsity.get_rowstart_indices()[row+1]],
-                                                 block_index.second);
-
-                    // check whether this line has
-                    // an entry in the regarding column
-                    // (check for ==dof_number and
-                    // != next_row, since if
-                    // row==dof_number-1, *p is a
-                    // past-the-end pointer but points
-                    // to dof_number anyway...)
-                    //
-                    // there should be
-                    // such an entry!
-                    // note, however,
-                    // that this
-                    // assertion will
-                    // fail sometimes
-                    // if the sparsity
-                    // pattern is not
-                    // symmetric!
-                    Assert ((*p == block_index.second) &&
-                            (p != &this_sparsity.get_column_numbers()
-                             [this_sparsity.get_rowstart_indices()[row+1]]),
-                            ExcInternalError());
-
-                    const unsigned int global_entry
-                      = (p
-                         -
-                         &this_sparsity.get_column_numbers()
-                         [this_sparsity.get_rowstart_indices()[0]]);
-
-                    // set matrix entry to zero
-                    matrix.block(block_row,block_index.first).global_entry(global_entry) = 0.;
-                  }
-              }
-          }
-      }
+    // this function is not documented and not tested in the testsuite
+    // so it isn't quite clear what it's supposed to do. it also isn't
+    // used anywhere else in the library. in avoiding the use of
+    // a deprecated function, I therefore threw away the original function
+    // and replaced it by the following, which I believe should work
+
+    std::map<unsigned int, double> boundary_values;
+    for (std::set<unsigned int>::const_iterator p=boundary_dofs.begin();
+        p != boundary_dofs.end(); ++p)
+      boundary_values[*p] = 0;
+
+    BlockVector<number> dummy(matrix.n_block_rows());
+    for (unsigned int i=0; i<matrix.n_block_rows(); ++i)
+      dummy.block(i).reinit (matrix.block(i,i).m());
+    dummy.collect_sizes();
+
+    MatrixTools::apply_boundary_values (boundary_values,
+                                        matrix, dummy, dummy,
+                                        preserve_symmetry);
   }
 }
 
-- 
2.39.5