From d42fc27dc1c223718812005a8ba05825d97a73a8 Mon Sep 17 00:00:00 2001 From: kanschat Date: Wed, 29 Sep 2010 00:41:25 +0000 Subject: [PATCH] implementation of MGMatrixBlockVector git-svn-id: https://svn.dealii.org/trunk@22185 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/matrix_block.h | 260 +++++++++++++++++-------- 1 file changed, 180 insertions(+), 80 deletions(-) diff --git a/deal.II/lac/include/lac/matrix_block.h b/deal.II/lac/include/lac/matrix_block.h index e44ddc07eb..260d493855 100644 --- a/deal.II/lac/include/lac/matrix_block.h +++ b/deal.II/lac/include/lac/matrix_block.h @@ -27,6 +27,19 @@ DEAL_II_NAMESPACE_OPEN +template class MatrixBlock; + +namespace internal +{ + template + void + reinit(MatrixBlock& v, const BlockSparsityPattern& p); + + template + void + reinit(MatrixBlock >&, const BlockSparsityPattern&); +} + /** * A wrapper around a matrix object, storing the coordinates in a * block matrix as well. @@ -72,6 +85,10 @@ DEAL_II_NAMESPACE_OPEN * Here, we have not gained very much, except that we do not need to * set up empty blocks in the block system. * + * @note This class expects, that the row and column BlockIndices + * objects for the system are equal. If they are not, some functions + * will throw ExcNotImplemented. + * * @todo Example for the product preconditioner of the pressure Schur * complement. * @@ -82,6 +99,7 @@ DEAL_II_NAMESPACE_OPEN */ template class MatrixBlock + : public Subscriptor { public: /** @@ -101,9 +119,22 @@ class MatrixBlock * initializing the matrix. */ - MatrixBlock(unsigned int i, unsigned int j, - const BlockIndices* block_indices = 0); + MatrixBlock(unsigned int i, unsigned int j); + /** + * Reinitialize the matrix for a + * new BlockSparsityPattern. This + * adujusts the #matrix as well + * as the #block_indices. + * + * @note The row and column block + * structure of the sparsity + * pattern must be equal. + */ + void reinit(const BlockSparsityPattern& sparsity); + + operator MATRIX() const; + /** * Add value to the * element (i,j). Throws @@ -327,7 +358,7 @@ class MatrixBlock * The matrix itself */ MATRIX matrix; - + private: /** * The BlockIndices of the whole * system. Using row() and @@ -336,7 +367,9 @@ class MatrixBlock * row and column degrees of * freedom for this block. */ - SmartPointer > block_indices; + BlockIndices block_indices; + + friend void internal::reinit<>(MatrixBlock&, const BlockSparsityPattern&); }; @@ -358,9 +391,7 @@ class MatrixBlockVector : public NamedData(row,column) * in the block system. */ - void add(unsigned int row, unsigned int column, - const std::string& name, - const BlockIndices* block_indices); + void add(unsigned int row, unsigned int column, const std::string& name); /** * For matrices using a @@ -371,6 +402,25 @@ class MatrixBlockVector : public NamedData -class MGMatrixBlockVector : public NamedData > > > +class MGMatrixBlockVector + : public NamedData > > > { public: /** @@ -419,15 +470,10 @@ class MGMatrixBlockVector : public NamedData& sparsity); + + /** + * Clears the object. + * + * Since often only clearing of + * the individual matrices is + * desired, but not removing the + * blocks themselves, there is an + * optional argument. If the + * argument is missing or @p + * false, all matrices will be + * mepty, but the size of this + * object and the block positions + * will not change. If @p + * really_clean is @p true, then + * the object will contain no + * blocks at the end. + */ + void clear (bool really_clean = false); + /** * Access a constant reference to * the block at position i. @@ -452,15 +518,6 @@ class MGMatrixBlockVector : public NamedDatai for read and write - * access. - */ - MATRIX& matrix(unsigned int i); - - }; @@ -470,41 +527,24 @@ namespace internal { template void - reinit(MatrixBlockVector&, const BlockSparsityPattern&) + reinit(MatrixBlock& v, const BlockSparsityPattern& p) { - Assert(false, ExcNotImplemented()); + Assert(p.get_row_indices() == p.get_column_indices(), ExcNotImplemented()); + v.block_indices = p.get_row_indices(); } template void - reinit(MatrixBlockVector >& v, const BlockSparsityPattern& p) + reinit(MatrixBlock >& v, const BlockSparsityPattern& p) { - for (unsigned int i=0;imatrix.reinit(p.block(v(i)->row, v(i)->column)); + Assert(p.get_row_indices() == p.get_column_indices(), ExcNotImplemented()); + v.block_indices = p.get_row_indices(); + v.matrix.reinit(p.block(v.row, v.column)); } - - - template - void - reinit(MGMatrixBlockVector&, const MGLevelObject&) - { - Assert(false, ExcNotImplemented()); - } - - - template - void - reinit(MGMatrixBlockVector >& v, - const MGLevelObject& sparsity) - { - for (unsigned int i=0;imatrix[level].reinit(sparsity[level].block(v(i)->row, v(i)->column)); - } - } + template inline MatrixBlock::MatrixBlock() @@ -518,6 +558,7 @@ template inline MatrixBlock::MatrixBlock(const MatrixBlock& M) : + Subscriptor(), row(M.row), column(M.column), matrix(M.matrix), @@ -527,13 +568,29 @@ MatrixBlock::MatrixBlock(const MatrixBlock& M) template inline -MatrixBlock::MatrixBlock(unsigned int i, unsigned int j, - const BlockIndices* block_indices) +MatrixBlock::MatrixBlock(unsigned int i, unsigned int j) : - row(i), column(j), block_indices(block_indices) + row(i), column(j) {} +template +inline +void +MatrixBlock::reinit(const BlockSparsityPattern& sparsity) +{ + internal::reinit(*this, sparsity); +} + + +template +inline +MatrixBlock::operator MATRIX() const +{ + return matrix; +} + + template inline void MatrixBlock::add ( @@ -541,12 +598,12 @@ MatrixBlock::add ( const unsigned int gj, const typename MATRIX::value_type value) { - Assert(block_indices != 0, ExcNotInitialized()); + Assert(block_indices.size() != 0, ExcNotInitialized()); const std::pair bi - = block_indices->global_to_local(gi); + = block_indices.global_to_local(gi); const std::pair bj - = block_indices->global_to_local(gj); + = block_indices.global_to_local(gj); Assert (bi.first == row, ExcBlockIndexMismatch(bi.first, row)); Assert (bj.first == column, ExcBlockIndexMismatch(bj.first, column)); @@ -564,7 +621,7 @@ MatrixBlock::add (const std::vector& row_indices, const FullMatrix& values, const bool elide_zero_values) { - Assert(block_indices != 0, ExcNotInitialized()); + Assert(block_indices.size() != 0, ExcNotInitialized()); AssertDimension (row_indices.size(), values.m()); AssertDimension (col_indices.size(), values.n()); @@ -586,9 +643,9 @@ MatrixBlock::add (const unsigned int b_row, const bool, const bool) { - Assert(block_indices != 0, ExcNotInitialized()); + Assert(block_indices.size() != 0, ExcNotInitialized()); const std::pair bi - = block_indices->global_to_local(b_row); + = block_indices.global_to_local(b_row); // In debug mode, we check whether // all indices are in the correct @@ -604,7 +661,7 @@ MatrixBlock::add (const unsigned int b_row, for (unsigned int j=0;j bj - = block_indices->global_to_local(col_indices[j]); + = block_indices.global_to_local(col_indices[j]); Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column)); matrix.add(bi.second, bj.second, values[j]); @@ -621,7 +678,7 @@ MatrixBlock::add (const std::vector &indices, const FullMatrix &values, const bool elide_zero_values) { - Assert(block_indices != 0, ExcNotInitialized()); + Assert(block_indices.size() != 0, ExcNotInitialized()); AssertDimension (indices.size(), values.m()); Assert (values.n() == values.m(), ExcNotQuadratic()); @@ -642,7 +699,7 @@ MatrixBlock::add (const unsigned int row, const std::vector &values, const bool elide_zero_values) { - Assert(block_indices != 0, ExcNotInitialized()); + Assert(block_indices.size() != 0, ExcNotInitialized()); AssertDimension (col_indices.size(), values.size()); add (row, col_indices.size(), &col_indices[0], &values[0], elide_zero_values); @@ -706,10 +763,9 @@ template inline void MatrixBlockVector::add( unsigned int row, unsigned int column, - const std::string& name, - const BlockIndices* block_indices) + const std::string& name) { - std_cxx1x::shared_ptr > p(new MatrixBlock(row, column, block_indices)); + std_cxx1x::shared_ptr > p(new MatrixBlock(row, column)); NamedData > >::add(p, name); } @@ -718,7 +774,26 @@ template inline void MatrixBlockVector::reinit(const BlockSparsityPattern& sparsity) { - internal::reinit(*this, sparsity); + for (unsigned int i=0;isize();++i) + { + block(i).reinit(sparsity); + } +} + + +template +inline void +MatrixBlockVector::clear(bool really_clean) +{ + if (really_clean) + { + Assert(false, ExcNotImplemented()); + } + else + { + for (unsigned int i=0;isize();++i) + matrix(i).clear(); + } } @@ -754,24 +829,14 @@ template inline void MGMatrixBlockVector::add( unsigned int row, unsigned int column, - const std::string& name, - const BlockIndices* block_indices) + const std::string& name) { std_cxx1x::shared_ptr > > - p(new MGLevelObject >(row, column, block_indices)); + p(new MGLevelObject >(row, column)); NamedData > >::add(p, name); } -template -inline void -MGMatrixBlockVector::reinit(const MGLevelObject& sparsity) -{ - internal::reinit(*this, sparsity); -} - - - template inline const MGLevelObject >& MGMatrixBlockVector::block(unsigned int i) const @@ -789,13 +854,48 @@ MGMatrixBlockVector::block(unsigned int i) template -inline MATRIX& -MGMatrixBlockVector::matrix(unsigned int i) +inline void +MGMatrixBlockVector::reinit(const MGLevelObject& sparsity) { - return (*this)(i)->matrix; + for (unsigned int i=0;isize();++i) + { + MGLevelObject >& o = block(i); + const unsigned int row = o[o.min_level()].row; + const unsigned int col = o[o.min_level()].column; + + o.resize(sparsity.min_level(), sparsity.max_level()); + for (unsigned int level = o.min_level();level <= o.max_level();++level) + { + o[level].row = row; + o[level].column = col; + internal::reinit(o[level], sparsity[level]); + } + } } +template +inline void +MGMatrixBlockVector::clear(bool really_clean) +{ + if (really_clean) + { + Assert(false, ExcNotImplemented()); + } + else + { + for (unsigned int i=0;isize();++i) + { + MGLevelObject >& o = block(i); + for (unsigned int level = o.min_level();level <= o.max_level();++level) + o[level].matrix.clear(); + } + } +} + + + + DEAL_II_NAMESPACE_CLOSE -- 2.39.5