--- /dev/null
+//---------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------------------------------------------------------
+
+
+#include <base/logstream.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
+#include <lac/block_vector.h>
+#include <multigrid/mg_tools.h>
+#include <multigrid/mg_base.h>
+#include <multigrid/mg_level_object.h>
+
+#include <vector>
+#include <algorithm>
+#include <numeric>
+
+
+template <typename number>
+void
+MGTools::apply_boundary_values (
+ const std::vector<unsigned int> &boundary_values,
+ SparseMatrix<number>& matrix,
+ const bool preserve_symmetry)
+{
+ // if no boundary values are to be applied
+ // simply return
+ if (boundary_values.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::vector<unsigned int>::const_iterator dof = boundary_values.begin(),
+ endd = boundary_values.end();
+ const SparsityPattern &sparsity = matrix.get_sparsity_pattern();
+ const unsigned int *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
+ 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 = std::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.;
+ }
+ }
+ }
+}
+
+
+
+template <typename number>
+void
+MGTools::apply_boundary_values (
+ const std::vector<unsigned int>& boundary_values,
+ BlockSparseMatrix<number>& matrix,
+ const bool preserve_symmetry)
+{
+ const unsigned int blocks = matrix.n_block_rows();
+
+ Assert (matrix.n_block_rows() == matrix.n_block_cols(),
+ ExcNotQuadratic());
+ Assert (matrix.get_sparsity_pattern().get_row_indices() ==
+ matrix.get_sparsity_pattern().get_column_indices(),
+ ExcNotQuadratic());
+
+ for (unsigned int i=0; i<blocks; ++i)
+ Assert (matrix.block(i,i).get_sparsity_pattern().optimize_diagonal(),
+ SparsityPattern::ExcDiagonalNotOptimized());
+
+
+ // if no boundary values are to be applied
+ // simply return
+ if (boundary_values.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::vector<unsigned int>::const_iterator dof = boundary_values.begin(),
+ endd = boundary_values.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 = std::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 = std::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.;
+ }
+ }
+ }
+ }
+}
+
+
+template void MGTools::apply_boundary_values (
+ const std::vector<unsigned int>& boundary_values,
+ SparseMatrix<float>& matrix,
+ const bool preserve_symmetry);
+template void MGTools::apply_boundary_values (
+ const std::vector<unsigned int>& boundary_values,
+ SparseMatrix<double>& matrix,
+ const bool preserve_symmetry);
+template void MGTools::apply_boundary_values (
+ const std::vector<unsigned int>& boundary_values,
+ BlockSparseMatrix<float>& matrix,
+ const bool preserve_symmetry);
+template void MGTools::apply_boundary_values (
+ const std::vector<unsigned int>& boundary_values,
+ BlockSparseMatrix<double>& matrix,
+ const bool preserve_symmetry);