* @ingroup numerics
* @author Wolfgang Bangerth, 1998, 2000, 2004, 2005
*/
-class MatrixTools
+namespace MatrixTools
{
-// using namespace MatrixCreator
- public:
- /**
- * Apply dirichlet boundary conditions
- * to the system matrix and vectors
- * as described in the general
- * documentation.
- */
- template <typename number>
- static void
- apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- SparseMatrix<number> &matrix,
- Vector<number> &solution,
- Vector<number> &right_hand_side,
- const bool eliminate_columns = true);
+ /**
+ * Import namespace MatrixCreator for
+ * backward compatibility with older
+ * versions of deal.II in which these
+ * namespaces were classes and class
+ * MatrixTools was publicly derived from
+ * class MatrixCreator.
+ */
+ using namespace MatrixCreator;
- /**
- * Apply dirichlet boundary
- * conditions to the system
- * matrix and vectors as
- * described in the general
- * documentation. This function
- * works for block sparse
- * matrices and block vectors
- */
- template <typename number>
- static void
- apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- BlockSparseMatrix<number> &matrix,
- BlockVector<number> &solution,
- BlockVector<number> &right_hand_side,
- const bool eliminate_columns = true);
+ /**
+ * Apply dirichlet boundary conditions
+ * to the system matrix and vectors
+ * as described in the general
+ * documentation.
+ */
+ template <typename number>
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ SparseMatrix<number> &matrix,
+ Vector<number> &solution,
+ Vector<number> &right_hand_side,
+ const bool eliminate_columns = true);
+
+ /**
+ * Apply dirichlet boundary
+ * conditions to the system
+ * matrix and vectors as
+ * described in the general
+ * documentation. This function
+ * works for block sparse
+ * matrices and block vectors
+ */
+ template <typename number>
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ BlockSparseMatrix<number> &matrix,
+ BlockVector<number> &solution,
+ BlockVector<number> &right_hand_side,
+ const bool eliminate_columns = true);
#ifdef DEAL_II_USE_PETSC
- /**
- * Apply dirichlet boundary conditions to
- * the system matrix and vectors as
- * described in the general
- * documentation. This function works on
- * the classes that are used to wrap
- * PETSc objects.
- *
- * Note that this function is not very
- * efficient: it needs to alternatingly
- * read and write into the matrix, a
- * situation that PETSc does not handle
- * too well. In addition, we only get rid
- * of rows corresponding to boundary
- * nodes, but the corresponding case of
- * deleting the respective columns
- * (i.e. if @p eliminate_columns is @p
- * true) is not presently implemented,
- * and probably will never because it is
- * too expensive without direct access to
- * the PETSc data structures. (This leads
- * to the situation where the action
- * indicates by the default value of the
- * last argument is actually not
- * implemented; that argument has
- * <code>true</code> as its default value
- * to stay consistent with the other
- * functions of same name in this class.)
- * A third reason against this function
- * is that it doesn't handle the case
- * where the matrix is distributed across
- * an MPI system.
- *
- * This function is used in
- * step-17 and
- * step-18.
- */
- static void
- apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- PETScWrappers::SparseMatrix &matrix,
- PETScWrappers::Vector &solution,
- PETScWrappers::Vector &right_hand_side,
- const bool eliminate_columns = true);
+ /**
+ * Apply dirichlet boundary conditions to
+ * the system matrix and vectors as
+ * described in the general
+ * documentation. This function works on
+ * the classes that are used to wrap
+ * PETSc objects.
+ *
+ * Note that this function is not very
+ * efficient: it needs to alternatingly
+ * read and write into the matrix, a
+ * situation that PETSc does not handle
+ * too well. In addition, we only get rid
+ * of rows corresponding to boundary
+ * nodes, but the corresponding case of
+ * deleting the respective columns
+ * (i.e. if @p eliminate_columns is @p
+ * true) is not presently implemented,
+ * and probably will never because it is
+ * too expensive without direct access to
+ * the PETSc data structures. (This leads
+ * to the situation where the action
+ * indicates by the default value of the
+ * last argument is actually not
+ * implemented; that argument has
+ * <code>true</code> as its default value
+ * to stay consistent with the other
+ * functions of same name in this class.)
+ * A third reason against this function
+ * is that it doesn't handle the case
+ * where the matrix is distributed across
+ * an MPI system.
+ *
+ * This function is used in
+ * step-17 and
+ * step-18.
+ */
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ PETScWrappers::SparseMatrix &matrix,
+ PETScWrappers::Vector &solution,
+ PETScWrappers::Vector &right_hand_side,
+ const bool eliminate_columns = true);
- /**
- * Same function, but for parallel PETSc
- * matrices.
- */
- static void
- apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- PETScWrappers::MPI::SparseMatrix &matrix,
- PETScWrappers::MPI::Vector &solution,
- PETScWrappers::MPI::Vector &right_hand_side,
- const bool eliminate_columns = true);
+ /**
+ * Same function, but for parallel PETSc
+ * matrices.
+ */
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ PETScWrappers::MPI::SparseMatrix &matrix,
+ PETScWrappers::MPI::Vector &solution,
+ PETScWrappers::MPI::Vector &right_hand_side,
+ const bool eliminate_columns = true);
- /**
- * Same function, but for
- * parallel PETSc matrices. Note
- * that this function only
- * operates on the local range of
- * the parallel matrix, i.e. it
- * only eliminates rows
- * corresponding to degrees of
- * freedom for which the row is
- * stored on the present
- * processor. All other boundary
- * nodes are ignored, and it
- * doesn't matter whether they
- * are present in the first
- * argument to this function or
- * not. A consequence of this,
- * however, is that this function
- * has to be called from all
- * processors that participate in
- * sharing the contents of the
- * given matrices and vectors. It
- * is also implied that the local
- * range for all objects passed
- * to this function is the same.
- */
- static void
- apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- PETScWrappers::MPI::SparseMatrix &matrix,
- PETScWrappers::Vector &solution,
- PETScWrappers::MPI::Vector &right_hand_side,
- const bool eliminate_columns = true);
+ /**
+ * Same function, but for
+ * parallel PETSc matrices. Note
+ * that this function only
+ * operates on the local range of
+ * the parallel matrix, i.e. it
+ * only eliminates rows
+ * corresponding to degrees of
+ * freedom for which the row is
+ * stored on the present
+ * processor. All other boundary
+ * nodes are ignored, and it
+ * doesn't matter whether they
+ * are present in the first
+ * argument to this function or
+ * not. A consequence of this,
+ * however, is that this function
+ * has to be called from all
+ * processors that participate in
+ * sharing the contents of the
+ * given matrices and vectors. It
+ * is also implied that the local
+ * range for all objects passed
+ * to this function is the same.
+ */
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ PETScWrappers::MPI::SparseMatrix &matrix,
+ PETScWrappers::Vector &solution,
+ PETScWrappers::MPI::Vector &right_hand_side,
+ const bool eliminate_columns = true);
- /**
- * Same as above but for BlockSparseMatrix.
- */
- static void
- apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- PETScWrappers::MPI::BlockSparseMatrix &matrix,
- PETScWrappers::MPI::BlockVector &solution,
- PETScWrappers::MPI::BlockVector &right_hand_side,
- const bool eliminate_columns = true);
+ /**
+ * Same as above but for BlockSparseMatrix.
+ */
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ PETScWrappers::MPI::BlockSparseMatrix &matrix,
+ PETScWrappers::MPI::BlockVector &solution,
+ PETScWrappers::MPI::BlockVector &right_hand_side,
+ const bool eliminate_columns = true);
#endif
#ifdef DEAL_II_USE_TRILINOS
- /**
- * Apply dirichlet boundary
- * conditions to the system matrix
- * and vectors as described in the
- * general documentation. This
- * function works on the classes
- * that are used to wrap Trilinos
- * objects.
- *
- * Note that this function is not
- * very efficient: it needs to
- * alternatingly read and write
- * into the matrix, a situation
- * that Trilinos does not handle
- * too well. In addition, we only
- * get rid of rows corresponding to
- * boundary nodes, but the
- * corresponding case of deleting
- * the respective columns (i.e. if
- * @p eliminate_columns is @p true)
- * is not presently implemented,
- * and probably will never because
- * it is too expensive without
- * direct access to the Trilinos
- * data structures. (This leads to
- * the situation where the action
- * indicates by the default value
- * of the last argument is actually
- * not implemented; that argument
- * has <code>true</code> as its
- * default value to stay consistent
- * with the other functions of same
- * name in this class.) A third
- * reason against this function is
- * that it doesn't handle the case
- * where the matrix is distributed
- * across an MPI system.
- */
- static void
- apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- TrilinosWrappers::SparseMatrix &matrix,
- TrilinosWrappers::Vector &solution,
- TrilinosWrappers::Vector &right_hand_side,
- const bool eliminate_columns = true);
+ /**
+ * Apply dirichlet boundary
+ * conditions to the system matrix
+ * and vectors as described in the
+ * general documentation. This
+ * function works on the classes
+ * that are used to wrap Trilinos
+ * objects.
+ *
+ * Note that this function is not
+ * very efficient: it needs to
+ * alternatingly read and write
+ * into the matrix, a situation
+ * that Trilinos does not handle
+ * too well. In addition, we only
+ * get rid of rows corresponding to
+ * boundary nodes, but the
+ * corresponding case of deleting
+ * the respective columns (i.e. if
+ * @p eliminate_columns is @p true)
+ * is not presently implemented,
+ * and probably will never because
+ * it is too expensive without
+ * direct access to the Trilinos
+ * data structures. (This leads to
+ * the situation where the action
+ * indicates by the default value
+ * of the last argument is actually
+ * not implemented; that argument
+ * has <code>true</code> as its
+ * default value to stay consistent
+ * with the other functions of same
+ * name in this class.) A third
+ * reason against this function is
+ * that it doesn't handle the case
+ * where the matrix is distributed
+ * across an MPI system.
+ */
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ TrilinosWrappers::SparseMatrix &matrix,
+ TrilinosWrappers::Vector &solution,
+ TrilinosWrappers::Vector &right_hand_side,
+ const bool eliminate_columns = true);
- /**
- * This function does the same as
- * the one above, except now
- * working on block structures.
- */
- static void
- apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- TrilinosWrappers::BlockSparseMatrix &matrix,
- TrilinosWrappers::BlockVector &solution,
- TrilinosWrappers::BlockVector &right_hand_side,
- const bool eliminate_columns = true);
+ /**
+ * This function does the same as
+ * the one above, except now
+ * working on block structures.
+ */
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ TrilinosWrappers::BlockSparseMatrix &matrix,
+ TrilinosWrappers::BlockVector &solution,
+ TrilinosWrappers::BlockVector &right_hand_side,
+ const bool eliminate_columns = true);
- /**
- * Apply dirichlet boundary
- * conditions to the system matrix
- * and vectors as described in the
- * general documentation. This
- * function works on the classes
- * that are used to wrap Trilinos
- * objects.
- *
- * Note that this function is not
- * very efficient: it needs to
- * alternatingly read and write
- * into the matrix, a situation
- * that Trilinos does not handle
- * too well. In addition, we only
- * get rid of rows corresponding to
- * boundary nodes, but the
- * corresponding case of deleting
- * the respective columns (i.e. if
- * @p eliminate_columns is @p true)
- * is not presently implemented,
- * and probably will never because
- * it is too expensive without
- * direct access to the Trilinos
- * data structures. (This leads to
- * the situation where the action
- * indicates by the default value
- * of the last argument is actually
- * not implemented; that argument
- * has <code>true</code> as its
- * default value to stay consistent
- * with the other functions of same
- * name in this class.) This
- * function does work on MPI vector
- * types.
- */
- static void
- apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- TrilinosWrappers::SparseMatrix &matrix,
- TrilinosWrappers::MPI::Vector &solution,
- TrilinosWrappers::MPI::Vector &right_hand_side,
- const bool eliminate_columns = true);
+ /**
+ * Apply dirichlet boundary
+ * conditions to the system matrix
+ * and vectors as described in the
+ * general documentation. This
+ * function works on the classes
+ * that are used to wrap Trilinos
+ * objects.
+ *
+ * Note that this function is not
+ * very efficient: it needs to
+ * alternatingly read and write
+ * into the matrix, a situation
+ * that Trilinos does not handle
+ * too well. In addition, we only
+ * get rid of rows corresponding to
+ * boundary nodes, but the
+ * corresponding case of deleting
+ * the respective columns (i.e. if
+ * @p eliminate_columns is @p true)
+ * is not presently implemented,
+ * and probably will never because
+ * it is too expensive without
+ * direct access to the Trilinos
+ * data structures. (This leads to
+ * the situation where the action
+ * indicates by the default value
+ * of the last argument is actually
+ * not implemented; that argument
+ * has <code>true</code> as its
+ * default value to stay consistent
+ * with the other functions of same
+ * name in this class.) This
+ * function does work on MPI vector
+ * types.
+ */
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ TrilinosWrappers::SparseMatrix &matrix,
+ TrilinosWrappers::MPI::Vector &solution,
+ TrilinosWrappers::MPI::Vector &right_hand_side,
+ const bool eliminate_columns = true);
- /**
- * This function does the same as
- * the one above, except now working
- * on block structures.
- */
- static void
- apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- TrilinosWrappers::BlockSparseMatrix &matrix,
- TrilinosWrappers::MPI::BlockVector &solution,
- TrilinosWrappers::MPI::BlockVector &right_hand_side,
- const bool eliminate_columns = true);
+ /**
+ * This function does the same as
+ * the one above, except now working
+ * on block structures.
+ */
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ TrilinosWrappers::BlockSparseMatrix &matrix,
+ TrilinosWrappers::MPI::BlockVector &solution,
+ TrilinosWrappers::MPI::BlockVector &right_hand_side,
+ const bool eliminate_columns = true);
#endif
- /**
- * Rather than applying boundary
- * values to the global matrix
- * and vector after creating the
- * global matrix, this function
- * does so during assembly, by
- * modifying the local matrix and
- * vector contributions. If you
- * call this function on all
- * local contributions, the
- * resulting matrix will have the
- * same entries, and the final
- * call to
- * apply_boundary_values() on the
- * global system will not be
- * necessary.
- *
- * Since this function does not
- * have to work on the
- * complicated data structures of
- * sparse matrices, it is
- * relatively cheap. It may
- * therefore be a win if you have
- * many fixed degrees of freedom
- * (e.g. boundary nodes), or if
- * access to the sparse matrix is
- * expensive (e.g. for block
- * sparse matrices, or for PETSc
- * or trilinos
- * matrices). However, it doesn't
- * work as expected if there are
- * also hanging nodes to be
- * considered. More caveats are
- * listed in the general
- * documentation of this class.
- */
- static void
- local_apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- const std::vector<unsigned int> &local_dof_indices,
- FullMatrix<double> &local_matrix,
- Vector<double> &local_rhs,
- const bool eliminate_columns);
+ /**
+ * Rather than applying boundary
+ * values to the global matrix
+ * and vector after creating the
+ * global matrix, this function
+ * does so during assembly, by
+ * modifying the local matrix and
+ * vector contributions. If you
+ * call this function on all
+ * local contributions, the
+ * resulting matrix will have the
+ * same entries, and the final
+ * call to
+ * apply_boundary_values() on the
+ * global system will not be
+ * necessary.
+ *
+ * Since this function does not
+ * have to work on the
+ * complicated data structures of
+ * sparse matrices, it is
+ * relatively cheap. It may
+ * therefore be a win if you have
+ * many fixed degrees of freedom
+ * (e.g. boundary nodes), or if
+ * access to the sparse matrix is
+ * expensive (e.g. for block
+ * sparse matrices, or for PETSc
+ * or trilinos
+ * matrices). However, it doesn't
+ * work as expected if there are
+ * also hanging nodes to be
+ * considered. More caveats are
+ * listed in the general
+ * documentation of this class.
+ */
+ void
+ local_apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ const std::vector<unsigned int> &local_dof_indices,
+ FullMatrix<double> &local_matrix,
+ Vector<double> &local_rhs,
+ const bool eliminate_columns);
- /**
- * Exception
- */
- DeclException0 (ExcBlocksDontMatch);
-};
+ /**
+ * Exception
+ */
+ DeclException0 (ExcBlocksDontMatch);
+}
} // namespace MatrixCreator
-
-//TODO:[WB] I don't think that the optimized storage of diagonals is needed (GK)
-template <typename number>
-void
-MatrixTools::apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- SparseMatrix<number> &matrix,
- Vector<number> &solution,
- Vector<number> &right_hand_side,
- const bool eliminate_columns)
+namespace MatrixTools
{
- // Require that diagonals are first
- // in each row
- Assert (matrix.get_sparsity_pattern().optimize_diagonal(),
- typename SparsityPattern::ExcDiagonalNotOptimized());
- Assert (matrix.n() == right_hand_side.size(),
- ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
- Assert (matrix.n() == solution.size(),
- ExcDimensionMismatch(matrix.n(), solution.size()));
- // 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;
- }
+//TODO:[WB] I don't think that the optimized storage of diagonals is needed (GK)
+ template <typename number>
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ SparseMatrix<number> &matrix,
+ Vector<number> &solution,
+ Vector<number> &right_hand_side,
+ const bool eliminate_columns)
+ {
+ // Require that diagonals are first
+ // in each row
+ Assert (matrix.get_sparsity_pattern().optimize_diagonal(),
+ typename SparsityPattern::ExcDiagonalNotOptimized());
+ Assert (matrix.n() == right_hand_side.size(),
+ ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
+ Assert (matrix.n() == solution.size(),
+ ExcDimensionMismatch(matrix.n(), solution.size()));
+ // if no boundary values are to be applied
+ // simply return
+ if (boundary_values.size() == 0)
+ return;
- std::map<unsigned int,double>::const_iterator dof = boundary_values.begin(),
- endd = boundary_values.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->first < n_dofs, ExcInternalError());
-
- const unsigned int dof_number = dof->first;
- // 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
- number new_rhs;
- if (matrix.diag_element(dof_number) != 0.0)
- {
- new_rhs = dof->second * matrix.diag_element(dof_number);
- right_hand_side(dof_number) = new_rhs;
- }
- else
- {
- matrix.set (dof_number, dof_number,
- first_nonzero_diagonal_entry);
- new_rhs = dof->second * first_nonzero_diagonal_entry;
- right_hand_side(dof_number) = new_rhs;
- }
+ const unsigned int n_dofs = matrix.m();
- // 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 (eliminate_columns)
+ // 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)
{
- // store the only nonzero entry
- // of this line for the Gauss
- // elimination step
- const number diagonal_entry = matrix.diag_element(dof_number);
-
- // 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
- right_hand_side(row) -= matrix.global_entry(global_entry) /
- diagonal_entry * new_rhs;
-
- // set matrix entry to zero
- matrix.global_entry(global_entry) = 0.;
- }
+ first_nonzero_diagonal_entry = matrix.diag_element(i);
+ break;
}
- // preset solution vector
- solution(dof_number) = dof->second;
- }
-}
-
-
-template <typename number>
-void
-MatrixTools::apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- BlockSparseMatrix<number> &matrix,
- BlockVector<number> &solution,
- BlockVector<number> &right_hand_side,
- const bool eliminate_columns)
-{
- const unsigned int blocks = matrix.n_block_rows();
-
- Assert (matrix.n() == right_hand_side.size(),
- ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
- Assert (matrix.n() == solution.size(),
- ExcDimensionMismatch(matrix.n(), solution.size()));
- 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());
- Assert (matrix.get_sparsity_pattern().get_column_indices() ==
- solution.get_block_indices (),
- ExcBlocksDontMatch ());
- Assert (matrix.get_sparsity_pattern().get_row_indices() ==
- right_hand_side.get_block_indices (),
- ExcBlocksDontMatch ());
-
- 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)
+ std::map<unsigned int,double>::const_iterator dof = boundary_values.begin(),
+ endd = boundary_values.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->first < n_dofs, ExcInternalError());
+
+ const unsigned int dof_number = dof->first;
+ // 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
+ number new_rhs;
+ if (matrix.diag_element(dof_number) != 0.0)
{
- first_nonzero_diagonal_entry
- = matrix.block(diag_block,diag_block).diag_element(i);
- break;
+ new_rhs = dof->second * matrix.diag_element(dof_number);
+ right_hand_side(dof_number) = new_rhs;
+ }
+ else
+ {
+ matrix.set (dof_number, dof_number,
+ first_nonzero_diagonal_entry);
+ new_rhs = dof->second * first_nonzero_diagonal_entry;
+ right_hand_side(dof_number) = new_rhs;
}
- // 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::map<unsigned int,double>::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->first < n_dofs, ExcInternalError());
-
- // get global index and index
- // in the block in which this
- // dof is located
- const unsigned int dof_number = dof->first;
- 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.;
- }
-
- // 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
- number new_rhs;
- if (matrix.block(block_index.first, block_index.first)
- .diag_element(block_index.second) != 0.0)
- new_rhs = dof->second *
- matrix.block(block_index.first, block_index.first)
- .diag_element(block_index.second);
- else
- {
- matrix.block(block_index.first, block_index.first)
- .diag_element(block_index.second)
- = first_nonzero_diagonal_entry;
- new_rhs = dof->second * first_nonzero_diagonal_entry;
- }
- right_hand_side.block(block_index.first)(block_index.second)
- = new_rhs;
-
-
- // 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 (eliminate_columns)
- {
- // store the only nonzero entry
- // of this line for the Gauss
- // elimination step
- const number diagonal_entry
- = matrix.block(block_index.first,block_index.first)
- .diag_element(block_index.second);
-
- // 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]]);
-
- // correct right hand side
- right_hand_side.block(block_row)(row)
- -= matrix.block(block_row,block_index.first).global_entry(global_entry) /
- diagonal_entry * new_rhs;
-
- // set matrix entry to zero
- matrix.block(block_row,block_index.first).global_entry(global_entry) = 0.;
- }
- }
- }
- // preset solution vector
- solution.block(block_index.first)(block_index.second) = dof->second;
- }
-}
+ // 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 (eliminate_columns)
+ {
+ // store the only nonzero entry
+ // of this line for the Gauss
+ // elimination step
+ const number diagonal_entry = matrix.diag_element(dof_number);
+
+ // 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
+ right_hand_side(row) -= matrix.global_entry(global_entry) /
+ diagonal_entry * new_rhs;
+
+ // set matrix entry to zero
+ matrix.global_entry(global_entry) = 0.;
+ }
+ }
+ // preset solution vector
+ solution(dof_number) = dof->second;
+ }
+ }
-#ifdef DEAL_II_USE_PETSC
-namespace PETScWrappers
-{
- template <typename PETScMatrix, typename PETScVector>
+ template <typename number>
void
apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- PETScMatrix &matrix,
- PETScVector &solution,
- PETScVector &right_hand_side,
- const bool eliminate_columns)
+ BlockSparseMatrix<number> &matrix,
+ BlockVector<number> &solution,
+ BlockVector<number> &right_hand_side,
+ const bool eliminate_columns)
{
- Assert (eliminate_columns == false, ExcNotImplemented());
+ const unsigned int blocks = matrix.n_block_rows();
Assert (matrix.n() == right_hand_side.size(),
- ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
+ ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
Assert (matrix.n() == solution.size(),
- ExcDimensionMismatch(matrix.n(), solution.size()));
+ ExcDimensionMismatch(matrix.n(), solution.size()));
+ 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());
+ Assert (matrix.get_sparsity_pattern().get_column_indices() ==
+ solution.get_block_indices (),
+ ExcBlocksDontMatch ());
+ Assert (matrix.get_sparsity_pattern().get_row_indices() ==
+ right_hand_side.get_block_indices (),
+ ExcBlocksDontMatch ());
+
+ 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 no boundary values are to be applied
+ // simply return
if (boundary_values.size() == 0)
return;
- const std::pair<unsigned int, unsigned int> local_range
- = matrix.local_range();
- Assert (local_range == right_hand_side.local_range(),
- ExcInternalError());
- Assert (local_range == solution.local_range(),
- ExcInternalError());
+ const unsigned int n_dofs = matrix.m();
- // we have to read and write from this
- // matrix (in this order). this will only
- // work if we compress the matrix first,
- // done here
- matrix.compress ();
+ // 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::map<unsigned int,double>::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->first < n_dofs, ExcInternalError());
+
+ // get global index and index
+ // in the block in which this
+ // dof is located
+ const unsigned int dof_number = dof->first;
+ 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.;
+ }
- // determine the first nonzero diagonal
- // entry from within the part of the
- // matrix that we can see. if we can't
- // find such an entry, take one
- PetscScalar average_nonzero_diagonal_entry = 1;
- for (unsigned int i=local_range.first; i<local_range.second; ++i)
- if (matrix.diag_element(i) != 0)
- {
- average_nonzero_diagonal_entry = std::fabs(matrix.diag_element(i));
- break;
- }
-
- // figure out which rows of the matrix we
- // have to eliminate on this processor
- std::vector<unsigned int> constrained_rows;
- for (std::map<unsigned int,double>::const_iterator
- dof = boundary_values.begin();
- dof != boundary_values.end();
- ++dof)
- if ((dof->first >= local_range.first) &&
- (dof->first < local_range.second))
- constrained_rows.push_back (dof->first);
-
- // then eliminate these rows and set
- // their diagonal entry to what we have
- // determined above. note that for petsc
- // matrices interleaving read with write
- // operations is very expensive. thus, we
- // here always replace the diagonal
- // element, rather than first checking
- // whether it is nonzero and in that case
- // preserving it. this is different from
- // the case of deal.II sparse matrices
- // treated in the other functions.
- matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
-
- // the next thing is to set right hand
- // side to the wanted value. there's one
- // drawback: if we write to individual
- // vector elements, then we have to do
- // that on all processors. however, some
- // processors may not need to set
- // anything because their chunk of
- // matrix/rhs do not contain any boundary
- // nodes. therefore, rather than using
- // individual calls, we use one call for
- // all elements, thereby making sure that
- // all processors call this function,
- // even if some only have an empty set of
- // elements to set
- right_hand_side.compress ();
- solution.compress ();
-
- std::vector<unsigned int> indices;
- std::vector<PetscScalar> solution_values;
- for (std::map<unsigned int,double>::const_iterator
- dof = boundary_values.begin();
- dof != boundary_values.end();
- ++dof)
- if ((dof->first >= local_range.first) &&
- (dof->first < local_range.second))
- {
- indices.push_back (dof->first);
- solution_values.push_back (dof->second);
- }
- solution.set (indices, solution_values);
-
- // now also set appropriate values for
- // the rhs
- for (unsigned int i=0; i<solution_values.size(); ++i)
- solution_values[i] *= average_nonzero_diagonal_entry;
-
- right_hand_side.set (indices, solution_values);
-
- // clean up
- matrix.compress ();
- solution.compress ();
- right_hand_side.compress ();
- }
-}
+ // 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
+ number new_rhs;
+ if (matrix.block(block_index.first, block_index.first)
+ .diag_element(block_index.second) != 0.0)
+ new_rhs = dof->second *
+ matrix.block(block_index.first, block_index.first)
+ .diag_element(block_index.second);
+ else
+ {
+ matrix.block(block_index.first, block_index.first)
+ .diag_element(block_index.second)
+ = first_nonzero_diagonal_entry;
+ new_rhs = dof->second * first_nonzero_diagonal_entry;
+ }
+ right_hand_side.block(block_index.first)(block_index.second)
+ = new_rhs;
+
+
+ // 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 (eliminate_columns)
+ {
+ // store the only nonzero entry
+ // of this line for the Gauss
+ // elimination step
+ const number diagonal_entry
+ = matrix.block(block_index.first,block_index.first)
+ .diag_element(block_index.second);
+
+ // 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]]);
-void
-MatrixTools::
-apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- PETScWrappers::SparseMatrix &matrix,
- PETScWrappers::Vector &solution,
- PETScWrappers::Vector &right_hand_side,
- const bool eliminate_columns)
-{
- // simply redirect to the generic function
- // used for both petsc matrix types
- PETScWrappers::apply_boundary_values (boundary_values, matrix, solution,
- right_hand_side, eliminate_columns);
-}
+ // correct right hand side
+ right_hand_side.block(block_row)(row)
+ -= matrix.block(block_row,block_index.first).global_entry(global_entry) /
+ diagonal_entry * new_rhs;
+ // set matrix entry to zero
+ matrix.block(block_row,block_index.first).global_entry(global_entry) = 0.;
+ }
+ }
+ }
+ // preset solution vector
+ solution.block(block_index.first)(block_index.second) = dof->second;
+ }
+ }
-void
-MatrixTools::
-apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- PETScWrappers::MPI::SparseMatrix &matrix,
- PETScWrappers::MPI::Vector &solution,
- PETScWrappers::MPI::Vector &right_hand_side,
- const bool eliminate_columns)
-{
- // simply redirect to the generic function
- // used for both petsc matrix types
- PETScWrappers::apply_boundary_values (boundary_values, matrix, solution,
- right_hand_side, eliminate_columns);
- // compress the matrix once we're done
- matrix.compress ();
-}
+#ifdef DEAL_II_USE_PETSC
-void
-MatrixTools::
-apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- PETScWrappers::MPI::BlockSparseMatrix &matrix,
- PETScWrappers::MPI::BlockVector &solution,
- PETScWrappers::MPI::BlockVector &right_hand_side,
- const bool eliminate_columns)
-{
- Assert (matrix.n() == right_hand_side.size(),
- ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
- Assert (matrix.n() == solution.size(),
- ExcDimensionMismatch(matrix.n(), solution.size()));
- Assert (matrix.n_block_rows() == matrix.n_block_cols(),
- ExcNotQuadratic());
-
- const unsigned int n_blocks = matrix.n_block_rows();
-
- matrix.compress();
-
- // We need to find the subdivision
- // into blocks for the boundary values.
- // To this end, generate a vector of
- // maps with the respective indices.
- std::vector<std::map<unsigned int,double> > block_boundary_values(n_blocks);
+ namespace internal
{
- int offset = 0, block=0;
- for (std::map<unsigned int,double>::const_iterator
- dof = boundary_values.begin();
- dof != boundary_values.end();
- ++dof)
+ namespace PETScWrappers
+ {
+ template <typename PETScMatrix, typename PETScVector>
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ PETScMatrix &matrix,
+ PETScVector &solution,
+ PETScVector &right_hand_side,
+ const bool eliminate_columns)
{
- if (dof->first >= matrix.block(block,0).m() + offset)
- {
- offset += matrix.block(block,0).m();
- block++;
- }
- const unsigned int index = dof->first - offset;
- block_boundary_values[block].insert(std::pair<unsigned int, double> (index,dof->second));
- }
- }
+ Assert (eliminate_columns == false, ExcNotImplemented());
- // Now call the non-block variants on
- // the diagonal subblocks and the
- // solution/rhs.
- for (unsigned int block=0; block<n_blocks; ++block)
- PETScWrappers::apply_boundary_values(block_boundary_values[block],
- matrix.block(block,block),
- solution.block(block),
- right_hand_side.block(block),
- eliminate_columns);
-
- // Finally, we need to do something
- // about the off-diagonal matrices. This
- // is luckily not difficult. Just clear
- // the whole row.
- for (unsigned int block_m=0; block_m<n_blocks; ++block_m)
- {
- const std::pair<unsigned int, unsigned int> local_range
- = matrix.block(block_m,0).local_range();
+ Assert (matrix.n() == right_hand_side.size(),
+ ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
+ Assert (matrix.n() == solution.size(),
+ ExcDimensionMismatch(matrix.n(), solution.size()));
- std::vector<unsigned int> constrained_rows;
- for (std::map<unsigned int,double>::const_iterator
- dof = block_boundary_values[block_m].begin();
- dof != block_boundary_values[block_m].end();
- ++dof)
- if ((dof->first >= local_range.first) &&
- (dof->first < local_range.second))
- constrained_rows.push_back (dof->first);
-
- for (unsigned int block_n=0; block_n<n_blocks; ++block_n)
- if (block_m != block_n)
- matrix.block(block_m,block_n).clear_rows(constrained_rows);
- }
-}
+ // if no boundary values are to be applied
+ // simply return
+ if (boundary_values.size() == 0)
+ return;
-#endif
+ const std::pair<unsigned int, unsigned int> local_range
+ = matrix.local_range();
+ Assert (local_range == right_hand_side.local_range(),
+ ExcInternalError());
+ Assert (local_range == solution.local_range(),
+ ExcInternalError());
+
+
+ // we have to read and write from this
+ // matrix (in this order). this will only
+ // work if we compress the matrix first,
+ // done here
+ matrix.compress ();
+
+ // determine the first nonzero diagonal
+ // entry from within the part of the
+ // matrix that we can see. if we can't
+ // find such an entry, take one
+ PetscScalar average_nonzero_diagonal_entry = 1;
+ for (unsigned int i=local_range.first; i<local_range.second; ++i)
+ if (matrix.diag_element(i) != 0)
+ {
+ average_nonzero_diagonal_entry = std::fabs(matrix.diag_element(i));
+ break;
+ }
+
+ // figure out which rows of the matrix we
+ // have to eliminate on this processor
+ std::vector<unsigned int> constrained_rows;
+ for (std::map<unsigned int,double>::const_iterator
+ dof = boundary_values.begin();
+ dof != boundary_values.end();
+ ++dof)
+ if ((dof->first >= local_range.first) &&
+ (dof->first < local_range.second))
+ constrained_rows.push_back (dof->first);
+ // then eliminate these rows and set
+ // their diagonal entry to what we have
+ // determined above. note that for petsc
+ // matrices interleaving read with write
+ // operations is very expensive. thus, we
+ // here always replace the diagonal
+ // element, rather than first checking
+ // whether it is nonzero and in that case
+ // preserving it. this is different from
+ // the case of deal.II sparse matrices
+ // treated in the other functions.
+ matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
+
+ // the next thing is to set right hand
+ // side to the wanted value. there's one
+ // drawback: if we write to individual
+ // vector elements, then we have to do
+ // that on all processors. however, some
+ // processors may not need to set
+ // anything because their chunk of
+ // matrix/rhs do not contain any boundary
+ // nodes. therefore, rather than using
+ // individual calls, we use one call for
+ // all elements, thereby making sure that
+ // all processors call this function,
+ // even if some only have an empty set of
+ // elements to set
+ right_hand_side.compress ();
+ solution.compress ();
+
+ std::vector<unsigned int> indices;
+ std::vector<PetscScalar> solution_values;
+ for (std::map<unsigned int,double>::const_iterator
+ dof = boundary_values.begin();
+ dof != boundary_values.end();
+ ++dof)
+ if ((dof->first >= local_range.first) &&
+ (dof->first < local_range.second))
+ {
+ indices.push_back (dof->first);
+ solution_values.push_back (dof->second);
+ }
+ solution.set (indices, solution_values);
+
+ // now also set appropriate values for
+ // the rhs
+ for (unsigned int i=0; i<solution_values.size(); ++i)
+ solution_values[i] *= average_nonzero_diagonal_entry;
+
+ right_hand_side.set (indices, solution_values);
+
+ // clean up
+ matrix.compress ();
+ solution.compress ();
+ right_hand_side.compress ();
+ }
+ }
+ }
-#ifdef DEAL_II_USE_TRILINOS
-namespace TrilinosWrappers
-{
- template <typename TrilinosMatrix, typename TrilinosVector>
void
apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- TrilinosMatrix &matrix,
- TrilinosVector &solution,
- TrilinosVector &right_hand_side,
- const bool eliminate_columns)
+ PETScWrappers::SparseMatrix &matrix,
+ PETScWrappers::Vector &solution,
+ PETScWrappers::Vector &right_hand_side,
+ const bool eliminate_columns)
{
- Assert (eliminate_columns == false, ExcNotImplemented());
+ // simply redirect to the generic function
+ // used for both petsc matrix types
+ internal::PETScWrappers::apply_boundary_values (boundary_values, matrix, solution,
+ right_hand_side, eliminate_columns);
+ }
- Assert (matrix.n() == right_hand_side.size(),
- ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
- Assert (matrix.n() == solution.size(),
- ExcDimensionMismatch(matrix.m(), solution.size()));
- // if no boundary values are to be applied
- // simply return
- if (boundary_values.size() == 0)
- return;
- const std::pair<unsigned int, unsigned int> local_range
- = matrix.local_range();
- Assert (local_range == right_hand_side.local_range(),
- ExcInternalError());
- Assert (local_range == solution.local_range(),
- ExcInternalError());
-
- // we have to read and write from this
- // matrix (in this order). this will only
- // work if we compress the matrix first,
- // done here
- matrix.compress ();
+ void
- // determine the first nonzero diagonal
- // entry from within the part of the
- // matrix that we can see. if we can't
- // find such an entry, take one
- TrilinosScalar average_nonzero_diagonal_entry = 1;
- for (unsigned int i=local_range.first; i<local_range.second; ++i)
- if (matrix.diag_element(i) != 0)
- {
- average_nonzero_diagonal_entry = std::fabs(matrix.diag_element(i));
- break;
- }
-
- // figure out which rows of the matrix we
- // have to eliminate on this processor
- std::vector<unsigned int> constrained_rows;
- for (std::map<unsigned int,double>::const_iterator
- dof = boundary_values.begin();
- dof != boundary_values.end();
- ++dof)
- if ((dof->first >= local_range.first) &&
- (dof->first < local_range.second))
- constrained_rows.push_back (dof->first);
-
- // then eliminate these rows and
- // set their diagonal entry to
- // what we have determined
- // above. if the value already is
- // nonzero, it will be preserved,
- // in accordance with the basic
- // matrix classes in deal.II.
- matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
-
- // the next thing is to set right
- // hand side to the wanted
- // value. there's one drawback:
- // if we write to individual
- // vector elements, then we have
- // to do that on all
- // processors. however, some
- // processors may not need to set
- // anything because their chunk
- // of matrix/rhs do not contain
- // any boundary nodes. therefore,
- // rather than using individual
- // calls, we use one call for all
- // elements, thereby making sure
- // that all processors call this
- // function, even if some only
- // have an empty set of elements
- // to set
- right_hand_side.compress ();
- solution.compress ();
-
- std::vector<unsigned int> indices;
- std::vector<TrilinosScalar> solution_values;
- for (std::map<unsigned int,double>::const_iterator
- dof = boundary_values.begin();
- dof != boundary_values.end();
- ++dof)
- if ((dof->first >= local_range.first) &&
- (dof->first < local_range.second))
- {
- indices.push_back (dof->first);
- solution_values.push_back (dof->second);
- }
- solution.set (indices, solution_values);
-
- // now also set appropriate
- // values for the rhs
- for (unsigned int i=0; i<solution_values.size(); ++i)
- solution_values[i] *= matrix.diag_element(indices[i]);
-
- right_hand_side.set (indices, solution_values);
-
- // clean up
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ PETScWrappers::MPI::SparseMatrix &matrix,
+ PETScWrappers::MPI::Vector &solution,
+ PETScWrappers::MPI::Vector &right_hand_side,
+ const bool eliminate_columns)
+ {
+ // simply redirect to the generic function
+ // used for both petsc matrix types
+ internal::PETScWrappers::apply_boundary_values (boundary_values, matrix, solution,
+ right_hand_side, eliminate_columns);
+
+ // compress the matrix once we're done
matrix.compress ();
- solution.compress ();
- right_hand_side.compress ();
}
-
- template <typename TrilinosMatrix, typename TrilinosBlockVector>
void
- apply_block_boundary_values (const std::map<unsigned int,double> &boundary_values,
- TrilinosMatrix &matrix,
- TrilinosBlockVector &solution,
- TrilinosBlockVector &right_hand_side,
- const bool eliminate_columns)
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ PETScWrappers::MPI::BlockSparseMatrix &matrix,
+ PETScWrappers::MPI::BlockVector &solution,
+ PETScWrappers::MPI::BlockVector &right_hand_side,
+ const bool eliminate_columns)
{
- Assert (eliminate_columns == false, ExcNotImplemented());
-
Assert (matrix.n() == right_hand_side.size(),
ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
Assert (matrix.n() == solution.size(),
matrix.compress();
- // We need to find the subdivision
- // into blocks for the boundary values.
- // To this end, generate a vector of
- // maps with the respective indices.
+ // We need to find the subdivision
+ // into blocks for the boundary values.
+ // To this end, generate a vector of
+ // maps with the respective indices.
std::vector<std::map<unsigned int,double> > block_boundary_values(n_blocks);
{
int offset = 0, block=0;
block++;
}
const unsigned int index = dof->first - offset;
- block_boundary_values[block].insert(
- std::pair<unsigned int, double> (index,dof->second));
+ block_boundary_values[block].insert(std::pair<unsigned int, double> (index,dof->second));
}
}
- // Now call the non-block variants on
- // the diagonal subblocks and the
- // solution/rhs.
+ // Now call the non-block variants on
+ // the diagonal subblocks and the
+ // solution/rhs.
for (unsigned int block=0; block<n_blocks; ++block)
- TrilinosWrappers::apply_boundary_values(block_boundary_values[block],
- matrix.block(block,block),
- solution.block(block),
- right_hand_side.block(block),
- eliminate_columns);
-
- // Finally, we need to do something
- // about the off-diagonal matrices. This
- // is luckily not difficult. Just clear
- // the whole row.
+ PETScWrappers::apply_boundary_values(block_boundary_values[block],
+ matrix.block(block,block),
+ solution.block(block),
+ right_hand_side.block(block),
+ eliminate_columns);
+
+ // Finally, we need to do something
+ // about the off-diagonal matrices. This
+ // is luckily not difficult. Just clear
+ // the whole row.
for (unsigned int block_m=0; block_m<n_blocks; ++block_m)
{
const std::pair<unsigned int, unsigned int> local_range
}
}
-}
+#endif
-void
-MatrixTools::
-apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- TrilinosWrappers::SparseMatrix &matrix,
- TrilinosWrappers::Vector &solution,
- TrilinosWrappers::Vector &right_hand_side,
- const bool eliminate_columns)
-{
- // simply redirect to the generic function
- // used for both trilinos matrix types
- TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution,
- right_hand_side, eliminate_columns);
-}
+#ifdef DEAL_II_USE_TRILINOS
+ namespace internal
+ {
+ namespace TrilinosWrappers
+ {
+ template <typename TrilinosMatrix, typename TrilinosVector>
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ TrilinosMatrix &matrix,
+ TrilinosVector &solution,
+ TrilinosVector &right_hand_side,
+ const bool eliminate_columns)
+ {
+ Assert (eliminate_columns == false, ExcNotImplemented());
+ Assert (matrix.n() == right_hand_side.size(),
+ ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
+ Assert (matrix.n() == solution.size(),
+ ExcDimensionMismatch(matrix.m(), solution.size()));
-void
-MatrixTools::
-apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- TrilinosWrappers::SparseMatrix &matrix,
- TrilinosWrappers::MPI::Vector &solution,
- TrilinosWrappers::MPI::Vector &right_hand_side,
- const bool eliminate_columns)
-{
- // simply redirect to the generic function
- // used for both trilinos matrix types
- TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution,
- right_hand_side, eliminate_columns);
-}
+ // if no boundary values are to be applied
+ // simply return
+ if (boundary_values.size() == 0)
+ return;
+ const std::pair<unsigned int, unsigned int> local_range
+ = matrix.local_range();
+ Assert (local_range == right_hand_side.local_range(),
+ ExcInternalError());
+ Assert (local_range == solution.local_range(),
+ ExcInternalError());
+
+ // we have to read and write from this
+ // matrix (in this order). this will only
+ // work if we compress the matrix first,
+ // done here
+ matrix.compress ();
+
+ // determine the first nonzero diagonal
+ // entry from within the part of the
+ // matrix that we can see. if we can't
+ // find such an entry, take one
+ TrilinosScalar average_nonzero_diagonal_entry = 1;
+ for (unsigned int i=local_range.first; i<local_range.second; ++i)
+ if (matrix.diag_element(i) != 0)
+ {
+ average_nonzero_diagonal_entry = std::fabs(matrix.diag_element(i));
+ break;
+ }
+ // figure out which rows of the matrix we
+ // have to eliminate on this processor
+ std::vector<unsigned int> constrained_rows;
+ for (std::map<unsigned int,double>::const_iterator
+ dof = boundary_values.begin();
+ dof != boundary_values.end();
+ ++dof)
+ if ((dof->first >= local_range.first) &&
+ (dof->first < local_range.second))
+ constrained_rows.push_back (dof->first);
-void
-MatrixTools::
-apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- TrilinosWrappers::BlockSparseMatrix &matrix,
- TrilinosWrappers::BlockVector &solution,
- TrilinosWrappers::BlockVector &right_hand_side,
- const bool eliminate_columns)
-{
- TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix,
- solution, right_hand_side,
- eliminate_columns);
-}
+ // then eliminate these rows and
+ // set their diagonal entry to
+ // what we have determined
+ // above. if the value already is
+ // nonzero, it will be preserved,
+ // in accordance with the basic
+ // matrix classes in deal.II.
+ matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
+
+ // the next thing is to set right
+ // hand side to the wanted
+ // value. there's one drawback:
+ // if we write to individual
+ // vector elements, then we have
+ // to do that on all
+ // processors. however, some
+ // processors may not need to set
+ // anything because their chunk
+ // of matrix/rhs do not contain
+ // any boundary nodes. therefore,
+ // rather than using individual
+ // calls, we use one call for all
+ // elements, thereby making sure
+ // that all processors call this
+ // function, even if some only
+ // have an empty set of elements
+ // to set
+ right_hand_side.compress ();
+ solution.compress ();
+
+ std::vector<unsigned int> indices;
+ std::vector<TrilinosScalar> solution_values;
+ for (std::map<unsigned int,double>::const_iterator
+ dof = boundary_values.begin();
+ dof != boundary_values.end();
+ ++dof)
+ if ((dof->first >= local_range.first) &&
+ (dof->first < local_range.second))
+ {
+ indices.push_back (dof->first);
+ solution_values.push_back (dof->second);
+ }
+ solution.set (indices, solution_values);
+ // now also set appropriate
+ // values for the rhs
+ for (unsigned int i=0; i<solution_values.size(); ++i)
+ solution_values[i] *= matrix.diag_element(indices[i]);
+ right_hand_side.set (indices, solution_values);
-void
-MatrixTools::
-apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- TrilinosWrappers::BlockSparseMatrix &matrix,
- TrilinosWrappers::MPI::BlockVector &solution,
- TrilinosWrappers::MPI::BlockVector &right_hand_side,
- const bool eliminate_columns)
-{
- TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix,
- solution, right_hand_side,
- eliminate_columns);
-}
+ // clean up
+ matrix.compress ();
+ solution.compress ();
+ right_hand_side.compress ();
+ }
-#endif
+ template <typename TrilinosMatrix, typename TrilinosBlockVector>
+ void
+ apply_block_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ TrilinosMatrix &matrix,
+ TrilinosBlockVector &solution,
+ TrilinosBlockVector &right_hand_side,
+ const bool eliminate_columns)
+ {
+ Assert (eliminate_columns == false, ExcNotImplemented());
-void
-MatrixTools::
-local_apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- const std::vector<unsigned int> &local_dof_indices,
- FullMatrix<double> &local_matrix,
- Vector<double> &local_rhs,
- const bool eliminate_columns)
-{
- Assert (local_dof_indices.size() == local_matrix.m(),
- ExcDimensionMismatch(local_dof_indices.size(),
- local_matrix.m()));
- Assert (local_dof_indices.size() == local_matrix.n(),
- ExcDimensionMismatch(local_dof_indices.size(),
- local_matrix.n()));
- Assert (local_dof_indices.size() == local_rhs.size(),
- ExcDimensionMismatch(local_dof_indices.size(),
- local_rhs.size()));
-
- // if there is nothing to do, then exit
- // right away
- if (boundary_values.size() == 0)
- return;
-
- // 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<n_local_dofs; ++i)
- {
- const std::map<unsigned int, double>::const_iterator
- boundary_value = boundary_values.find (local_dof_indices[i]);
- if (boundary_value != boundary_values.end())
- {
- // remove this row, except for the
- // diagonal element
- for (unsigned j=0; j<n_local_dofs; ++j)
- if (i != j)
- local_matrix(i,j) = 0;
-
- // replace diagonal entry by its
- // absolute value to make sure that
- // everything remains positive, or
- // by the average diagonal value if
- // zero
- if (local_matrix(i,i) == 0.)
- {
- // 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<n_local_dofs; ++k)
- if (local_matrix(k,k) != 0.)
- {
- average_diagonal += std::fabs(local_matrix(k,k));
- ++nonzero_diagonals;
- }
- if (nonzero_diagonals != 0)
- average_diagonal /= nonzero_diagonals;
- else
- average_diagonal = 0;
- }
-
- // only if all diagonal entries
- // are zero, then resort to the
- // last measure: choose one
- if (average_diagonal == 0.)
- average_diagonal = 1.;
-
- local_matrix(i,i) = average_diagonal;
- }
- else
- local_matrix(i,i) = std::fabs(local_matrix(i,i));
-
- // and replace rhs entry by correct
- // value
- local_rhs(i) = local_matrix(i,i) * boundary_value->second;
-
- // finally do the elimination step
- // if requested
- if (eliminate_columns == true)
- {
- for (unsigned int row=0; row<n_local_dofs; ++row)
- if (row != i)
- {
- local_rhs(row) -= local_matrix(row,i) * boundary_value->second;
- local_matrix(row,i) = 0;
- }
- }
- }
+ Assert (matrix.n() == right_hand_side.size(),
+ ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
+ Assert (matrix.n() == solution.size(),
+ ExcDimensionMismatch(matrix.n(), solution.size()));
+ Assert (matrix.n_block_rows() == matrix.n_block_cols(),
+ ExcNotQuadratic());
+
+ const unsigned int n_blocks = matrix.n_block_rows();
+
+ matrix.compress();
+
+ // We need to find the subdivision
+ // into blocks for the boundary values.
+ // To this end, generate a vector of
+ // maps with the respective indices.
+ std::vector<std::map<unsigned int,double> > block_boundary_values(n_blocks);
+ {
+ int offset = 0, block=0;
+ for (std::map<unsigned int,double>::const_iterator
+ dof = boundary_values.begin();
+ dof != boundary_values.end();
+ ++dof)
+ {
+ if (dof->first >= matrix.block(block,0).m() + offset)
+ {
+ offset += matrix.block(block,0).m();
+ block++;
+ }
+ const unsigned int index = dof->first - offset;
+ block_boundary_values[block].insert(
+ std::pair<unsigned int, double> (index,dof->second));
+ }
+ }
+
+ // Now call the non-block variants on
+ // the diagonal subblocks and the
+ // solution/rhs.
+ for (unsigned int block=0; block<n_blocks; ++block)
+ TrilinosWrappers::apply_boundary_values(block_boundary_values[block],
+ matrix.block(block,block),
+ solution.block(block),
+ right_hand_side.block(block),
+ eliminate_columns);
+
+ // Finally, we need to do something
+ // about the off-diagonal matrices. This
+ // is luckily not difficult. Just clear
+ // the whole row.
+ for (unsigned int block_m=0; block_m<n_blocks; ++block_m)
+ {
+ const std::pair<unsigned int, unsigned int> local_range
+ = matrix.block(block_m,0).local_range();
+
+ std::vector<unsigned int> constrained_rows;
+ for (std::map<unsigned int,double>::const_iterator
+ dof = block_boundary_values[block_m].begin();
+ dof != block_boundary_values[block_m].end();
+ ++dof)
+ if ((dof->first >= local_range.first) &&
+ (dof->first < local_range.second))
+ constrained_rows.push_back (dof->first);
+
+ for (unsigned int block_n=0; block_n<n_blocks; ++block_n)
+ if (block_m != block_n)
+ matrix.block(block_m,block_n).clear_rows(constrained_rows);
+ }
+ }
}
+ }
+
+
+
+
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ TrilinosWrappers::SparseMatrix &matrix,
+ TrilinosWrappers::Vector &solution,
+ TrilinosWrappers::Vector &right_hand_side,
+ const bool eliminate_columns)
+ {
+ // simply redirect to the generic function
+ // used for both trilinos matrix types
+ internal::TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution,
+ right_hand_side, eliminate_columns);
+ }
+
+
+
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ TrilinosWrappers::SparseMatrix &matrix,
+ TrilinosWrappers::MPI::Vector &solution,
+ TrilinosWrappers::MPI::Vector &right_hand_side,
+ const bool eliminate_columns)
+ {
+ // simply redirect to the generic function
+ // used for both trilinos matrix types
+ internal::TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution,
+ right_hand_side, eliminate_columns);
+ }
+
+
+
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ TrilinosWrappers::BlockSparseMatrix &matrix,
+ TrilinosWrappers::BlockVector &solution,
+ TrilinosWrappers::BlockVector &right_hand_side,
+ const bool eliminate_columns)
+ {
+ internal::TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix,
+ solution, right_hand_side,
+ eliminate_columns);
+ }
+
+
+
+ void
+ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ TrilinosWrappers::BlockSparseMatrix &matrix,
+ TrilinosWrappers::MPI::BlockVector &solution,
+ TrilinosWrappers::MPI::BlockVector &right_hand_side,
+ const bool eliminate_columns)
+ {
+ internal::TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix,
+ solution, right_hand_side,
+ eliminate_columns);
+ }
+
+#endif
+
+
+
+ void
+ local_apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ const std::vector<unsigned int> &local_dof_indices,
+ FullMatrix<double> &local_matrix,
+ Vector<double> &local_rhs,
+ const bool eliminate_columns)
+ {
+ Assert (local_dof_indices.size() == local_matrix.m(),
+ ExcDimensionMismatch(local_dof_indices.size(),
+ local_matrix.m()));
+ Assert (local_dof_indices.size() == local_matrix.n(),
+ ExcDimensionMismatch(local_dof_indices.size(),
+ local_matrix.n()));
+ Assert (local_dof_indices.size() == local_rhs.size(),
+ ExcDimensionMismatch(local_dof_indices.size(),
+ local_rhs.size()));
+
+ // if there is nothing to do, then exit
+ // right away
+ if (boundary_values.size() == 0)
+ return;
+
+ // 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<n_local_dofs; ++i)
+ {
+ const std::map<unsigned int, double>::const_iterator
+ boundary_value = boundary_values.find (local_dof_indices[i]);
+ if (boundary_value != boundary_values.end())
+ {
+ // remove this row, except for the
+ // diagonal element
+ for (unsigned j=0; j<n_local_dofs; ++j)
+ if (i != j)
+ local_matrix(i,j) = 0;
+
+ // replace diagonal entry by its
+ // absolute value to make sure that
+ // everything remains positive, or
+ // by the average diagonal value if
+ // zero
+ if (local_matrix(i,i) == 0.)
+ {
+ // 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<n_local_dofs; ++k)
+ if (local_matrix(k,k) != 0.)
+ {
+ average_diagonal += std::fabs(local_matrix(k,k));
+ ++nonzero_diagonals;
+ }
+ if (nonzero_diagonals != 0)
+ average_diagonal /= nonzero_diagonals;
+ else
+ average_diagonal = 0;
+ }
+
+ // only if all diagonal entries
+ // are zero, then resort to the
+ // last measure: choose one
+ if (average_diagonal == 0.)
+ average_diagonal = 1.;
+
+ local_matrix(i,i) = average_diagonal;
+ }
+ else
+ local_matrix(i,i) = std::fabs(local_matrix(i,i));
+
+ // and replace rhs entry by correct
+ // value
+ local_rhs(i) = local_matrix(i,i) * boundary_value->second;
+
+ // finally do the elimination step
+ // if requested
+ if (eliminate_columns == true)
+ {
+ for (unsigned int row=0; row<n_local_dofs; ++row)
+ if (row != i)
+ {
+ local_rhs(row) -= local_matrix(row,i) * boundary_value->second;
+ local_matrix(row,i) = 0;
+ }
+ }
+ }
+ }
+ }
}
// explicit instantiations
#include "matrices.inst"
-template
-void
-MatrixTools::apply_boundary_values<double> (const std::map<unsigned int,double> &boundary_values,
- SparseMatrix<double> &matrix,
- Vector<double> &solution,
- Vector<double> &right_hand_side,
- const bool eliminate_columns);
-template
-void
-MatrixTools::apply_boundary_values<float> (const std::map<unsigned int,double> &boundary_values,
- SparseMatrix<float> &matrix,
- Vector<float> &solution,
- Vector<float> &right_hand_side,
- const bool eliminate_columns);
-
-template
-void
-MatrixTools::apply_boundary_values<double> (const std::map<unsigned int,double> &boundary_values,
- BlockSparseMatrix<double> &matrix,
- BlockVector<double> &solution,
- BlockVector<double> &right_hand_side,
- const bool eliminate_columns);
-template
-void
-MatrixTools::apply_boundary_values<float> (const std::map<unsigned int,double> &boundary_values,
- BlockSparseMatrix<float> &matrix,
- BlockVector<float> &solution,
- BlockVector<float> &right_hand_side,
- const bool eliminate_columns);
+namespace MatrixTools
+{
+ template
+ void
+ apply_boundary_values<double> (const std::map<unsigned int,double> &boundary_values,
+ SparseMatrix<double> &matrix,
+ Vector<double> &solution,
+ Vector<double> &right_hand_side,
+ const bool eliminate_columns);
+ template
+ void
+ apply_boundary_values<float> (const std::map<unsigned int,double> &boundary_values,
+ SparseMatrix<float> &matrix,
+ Vector<float> &solution,
+ Vector<float> &right_hand_side,
+ const bool eliminate_columns);
+ template
+ void
+ apply_boundary_values<double> (const std::map<unsigned int,double> &boundary_values,
+ BlockSparseMatrix<double> &matrix,
+ BlockVector<double> &solution,
+ BlockVector<double> &right_hand_side,
+ const bool eliminate_columns);
+ template
+ void
+ apply_boundary_values<float> (const std::map<unsigned int,double> &boundary_values,
+ BlockSparseMatrix<float> &matrix,
+ BlockVector<float> &solution,
+ BlockVector<float> &right_hand_side,
+ const bool eliminate_columns);
+}
DEAL_II_NAMESPACE_CLOSE
projects / dealii-svn.git / commitdiff