#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>
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);
}
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);
}
}
projects / dealii-svn.git / commitdiff