From 89f9d2766e23c57a8780eacf37252a679a09eb49 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Sat, 20 Mar 2010 03:35:39 +0000 Subject: [PATCH] add add to MatrixBlock git-svn-id: https://svn.dealii.org/trunk@20862 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/block_indices.h | 5 +- .../include/lac/constraint_matrix.templates.h | 14 +- deal.II/lac/include/lac/matrix_block.h | 270 +++++++++++++++++- 3 files changed, 272 insertions(+), 17 deletions(-) diff --git a/deal.II/lac/include/lac/block_indices.h b/deal.II/lac/include/lac/block_indices.h index 1c228b84fa..0ee28c3227 100644 --- a/deal.II/lac/include/lac/block_indices.h +++ b/deal.II/lac/include/lac/block_indices.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -15,6 +15,7 @@ #include +#include #include #include @@ -32,7 +33,7 @@ DEAL_II_NAMESPACE_OPEN * @ingroup data * @author Wolfgang Bangerth, Guido Kanschat, 2000, 2007 */ -class BlockIndices +class BlockIndices : public Subscriptor { public: diff --git a/deal.II/lac/include/lac/constraint_matrix.templates.h b/deal.II/lac/include/lac/constraint_matrix.templates.h index 6df7772742..5f20f9465f 100644 --- a/deal.II/lac/include/lac/constraint_matrix.templates.h +++ b/deal.II/lac/include/lac/constraint_matrix.templates.h @@ -2022,13 +2022,13 @@ ConstraintMatrix:: // standard (non-block) matrices template void -ConstraintMatrix:: -distribute_local_to_global (const FullMatrix &local_matrix, - const Vector &local_vector, - const std::vector &local_dof_indices, - MatrixType &global_matrix, - VectorType &global_vector, - internal::bool2type) const +ConstraintMatrix::distribute_local_to_global ( + const FullMatrix &local_matrix, + const Vector &local_vector, + const std::vector &local_dof_indices, + MatrixType &global_matrix, + VectorType &global_vector, + internal::bool2type) const { // check whether we work on real vectors // or we just used a dummy when calling diff --git a/deal.II/lac/include/lac/matrix_block.h b/deal.II/lac/include/lac/matrix_block.h index 93d26832fa..03b9a64260 100644 --- a/deal.II/lac/include/lac/matrix_block.h +++ b/deal.II/lac/include/lac/matrix_block.h @@ -15,6 +15,9 @@ #include #include +#include +#include +#include #include @@ -51,8 +54,9 @@ DEAL_II_NAMESPACE_OPEN * @author Guido Kanschat, 2006 */ template -struct MatrixBlock +class MatrixBlock { + public: /** * Constructor rendering an * uninitialized object. @@ -70,7 +74,124 @@ struct MatrixBlock * initializing the matrix. */ - MatrixBlock(unsigned int i, unsigned int j); + MatrixBlock(unsigned int i, unsigned int j, + const BlockIndices* block_indices = 0); + + /** + * Add value to the + * element (i,j). Throws + * an error if the entry does not + * exist or if value is + * not a finite number. Still, it + * is allowed to store zero + * values in non-existent fields. + */ + void add (const unsigned int i, + const unsigned int j, + const typename MATRIX::value_type value); + + /** + * Add all elements given in a + * FullMatrix into sparse + * matrix locations given by + * indices. In other words, + * this function adds the elements in + * full_matrix to the + * respective entries in calling + * matrix, using the local-to-global + * indexing specified by + * indices for both the rows + * and the columns of the + * matrix. This function assumes a + * quadratic sparse matrix and a + * quadratic full_matrix, the usual + * situation in FE calculations. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be added anyway or + * these should be filtered away and + * only non-zero data is added. The + * default value is true, + * i.e., zero values won't be added + * into the matrix. + */ + template + void add (const std::vector& indices, + const FullMatrix& full_matrix, + const bool elide_zero_values = true); + + /** + * Same function as before, but now + * including the possibility to use + * rectangular full_matrices and + * different local-to-global indexing + * on rows and columns, respectively. + */ + template + void add (const std::vector& row_indices, + const std::vector& col_indices, + const FullMatrix& full_matrix, + const bool elide_zero_values = true); + + /** + * Set several elements in the + * specified row of the matrix with + * column indices as given by + * col_indices to the + * respective value. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be added anyway or + * these should be filtered away and + * only non-zero data is added. The + * default value is true, + * i.e., zero values won't be added + * into the matrix. + */ + template + void add (const unsigned int row, + const std::vector& col_indices, + const std::vector& values, + const bool elide_zero_values = true); + + /** + * Add an array of values given by + * values in the given + * global matrix row at columns + * specified by col_indices in the + * sparse matrix. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be added anyway or + * these should be filtered away and + * only non-zero data is added. The + * default value is true, + * i.e., zero values won't be added + * into the matrix. + */ + template + void add (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices, + const number *values, + const bool elide_zero_values = true, + const bool col_indices_are_sorted = false); + + /** + * The block number computed from + * an index by using + * #block_indices does not match + * the block coordinates stored + * in this object. + */ + DeclException2(ExcBlockIndexMismatch, unsigned int, unsigned int, + << "Block index " << arg1 << " does not match " << arg2); /** * Row coordinate. This is the @@ -90,6 +211,16 @@ struct MatrixBlock * The matrix itself */ MATRIX matrix; + + /** + * The BlockIndices of the whole + * system. Using row() and + * column(), this allows us to + * find the index of the first + * row and column degrees of + * freedom for this block. + */ + SmartPointer > block_indices; }; @@ -109,7 +240,9 @@ class MatrixBlockVector : public NamedData * position * in the block system. */ - void add(unsigned int row, unsigned int column, const std::string& name); + void add(unsigned int row, unsigned int column, + const std::string& name, + const BlockIndices* block_indices); /** * Access a constant reference to * the block at position i. @@ -129,6 +262,7 @@ class MatrixBlockVector : public NamedData //----------------------------------------------------------------------// template +inline MatrixBlock::MatrixBlock() : row(deal_II_numbers::invalid_unsigned_int), @@ -137,27 +271,147 @@ MatrixBlock::MatrixBlock() template +inline MatrixBlock::MatrixBlock(const MatrixBlock& M) : row(M.row), column(M.column), - matrix(M.matrix) + matrix(M.matrix), + block_indices(M.block_indices) {} template -MatrixBlock::MatrixBlock(unsigned int i, unsigned int j) +inline +MatrixBlock::MatrixBlock(unsigned int i, unsigned int j, + const BlockIndices* block_indices) : - row(i), column(j) + row(i), column(j), block_indices(block_indices) {} + +template +inline void +MatrixBlock::add ( + const unsigned int gi, + const unsigned int gj, + const typename MATRIX::value_type value) +{ + Assert(block_indices != 0, ExcNotInitialized()); + + const std::pair bi + = block_indices->global_to_local(gi); + const std::pair bj + = block_indices->global_to_local(gj); + + Assert (bi.first == row, ExcBlockIndexMismatch(bi.first, row)); + Assert (bj.first == column, ExcBlockIndexMismatch(bj.first, column)); + + matrix.add(bi.second, bj.second, value); +} + + +template +template +inline +void +MatrixBlock::add (const std::vector& row_indices, + const std::vector& col_indices, + const FullMatrix& values, + const bool elide_zero_values) +{ + AssertDimension (row_indices.size(), values.m()); + AssertDimension (col_indices.size(), values.n()); + + for (unsigned int i=0; i +template +inline +void +MatrixBlock::add (const std::vector &indices, + const FullMatrix &values, + const bool elide_zero_values) +{ + AssertDimension (indices.size(), values.m()); + Assert (values.n() == values.m(), ExcNotQuadratic()); + + for (unsigned int i=0; i +template +inline +void +MatrixBlock::add (const unsigned int row, + const std::vector &col_indices, + const std::vector &values, + const bool elide_zero_values) +{ + AssertDimension (col_indices.size(), values.size()); + add (row, col_indices.size(), &col_indices[0], &values[0], + elide_zero_values); +} + + +template +template +inline +void +MatrixBlock::add (const unsigned int b_row, + const unsigned int n_cols, + const unsigned int *col_indices, + const number *values, + const bool, + const bool) +{ + const std::pair bi + = block_indices->global_to_local(b_row); + + // In debug mode, we check whether + // all indices are in the correct + // block. + + // Actually, for the time being, we + // leave it at this. While it may + // not be the most efficient way, + // it is at least thread safe. +//#ifdef DEBUG + Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row)); + + for (unsigned int j=0;j bj + = block_indices->global_to_local(col_indices[j]); + Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column)); + + matrix.add(bi.second, bj.second, values[j]); + } + +//#endif + +} + + //----------------------------------------------------------------------// template inline void -MatrixBlockVector::add(unsigned int row, unsigned int column, const std::string& name) +MatrixBlockVector::add( + unsigned int row, unsigned int column, + const std::string& name, + const BlockIndices* block_indices) { - boost::shared_ptr > p(new MatrixBlock(row, column)); + boost::shared_ptr > p(new MatrixBlock(row, column, block_indices)); NamedData > >::add(p, name); } -- 2.39.5