From: kronbichler Date: Tue, 2 Sep 2008 14:27:26 +0000 (+0000) Subject: Added apply_boundary_conditions() function to Trilinos Block matrices and vectors. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c970fe9f4cc175834a508f2c7da5c35e4a5c1451;p=dealii-svn.git Added apply_boundary_conditions() function to Trilinos Block matrices and vectors. git-svn-id: https://svn.dealii.org/trunk@16707 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h index 476afa305d..9b8db14cba 100644 --- a/deal.II/deal.II/include/numerics/matrices.h +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -65,6 +65,8 @@ namespace TrilinosWrappers { class SparseMatrix; class Vector; + class BlockSparseMatrix; + class BlockVector; } #endif @@ -1137,10 +1139,6 @@ class MatrixTools : public MatrixCreator * is that it doesn't handle the case * where the matrix is distributed across * an MPI system. - * - * This function is used in - * @ref step_17 "step-17" and - * @ref step_18 "step-18". */ static void apply_boundary_values (const std::map &boundary_values, @@ -1148,6 +1146,18 @@ class MatrixTools : public MatrixCreator 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 &boundary_values, + TrilinosWrappers::BlockSparseMatrix &matrix, + TrilinosWrappers::BlockVector &solution, + TrilinosWrappers::BlockVector &right_hand_side, + const bool eliminate_columns = true); #endif /** diff --git a/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc b/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc index 3034f81236..9fcf748c6c 100644 --- a/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc +++ b/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc @@ -29,6 +29,8 @@ #ifdef DEAL_II_USE_TRILINOS # include # include +# include +# include #endif #include @@ -700,7 +702,7 @@ namespace TrilinosWrappers Assert (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.m(), solution.size())); // if no boundary values are to be applied // simply return @@ -713,7 +715,6 @@ namespace TrilinosWrappers ExcInternalError()); Assert (local_range == solution.local_range(), ExcInternalError()); - // we have to read and write from this // matrix (in this order). this will only @@ -786,7 +787,7 @@ namespace TrilinosWrappers // now also set appropriate values for // the rhs for (unsigned int i=0; i &boundary_values, right_hand_side, eliminate_columns); } +void +MatrixTools:: +apply_boundary_values (const std::map &boundary_values, + TrilinosWrappers::BlockSparseMatrix &matrix, + TrilinosWrappers::BlockVector &solution, + TrilinosWrappers::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 > block_boundary_values(n_blocks); + { + int offset = 0, block=0; + for (std::map::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 (index,dof->second)); + } + } + + // Now call the non-block variants on + // the diagonal subblocks and the + // solution/rhs. + for (unsigned int block=0; block local_range + = matrix.block(block_m,0).local_range(); + + std::vector constrained_rows; + for (std::map::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_nblock_indices = v.get_block_indices(); if (this->components.size() != this->n_blocks()) this->components.resize(this->n_blocks()); for (unsigned int i=0;in_blocks();++i) - block(i).reinit(v.block(i)); + block(i).reinit(v.block(i), fast); } diff --git a/deal.II/lac/source/trilinos_block_sparse_matrix.cc b/deal.II/lac/source/trilinos_block_sparse_matrix.cc index 04f43fe0f5..90731a9cc4 100644 --- a/deal.II/lac/source/trilinos_block_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_block_sparse_matrix.cc @@ -156,16 +156,31 @@ namespace TrilinosWrappers this->block(r,c).compress(); } + + void BlockSparseMatrix::collect_sizes () { compress(); BaseClass::collect_sizes (); } - + + + + unsigned int + BlockSparseMatrix::n_nonzero_elements () const + { + unsigned int n_nonzero = 0; + for (unsigned int rows = 0; rowsn_block_rows(); ++rows) + for (unsigned int cols = 0; colsn_block_cols(); ++cols) + n_nonzero += this->block(rows,cols).n_nonzero_elements(); + + return n_nonzero; + } + } - + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index f83268fe4e..b48139bc9f 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -492,10 +492,10 @@ namespace TrilinosWrappers int diag_index = (int)(diag_find - col_indices); for (int j=0; jFilled()) @@ -750,9 +750,9 @@ namespace TrilinosWrappers { Assert (&src != &dst, ExcSourceEqualsDestination()); - Assert (row_map.SameAs(dst.map), + Assert (row_map.SameAs(src.map), ExcMessage ("Row map of matrix does not fit with vector map!")); - Assert (col_map.SameAs(src.map), + Assert (col_map.SameAs(dst.map), ExcMessage ("Column map of matrix does not fit with vector map!")); if (!matrix->Filled())