#include <base/subscriptor.h>
class SparsityPattern;
-template <int rows, int columns> class BlockSparsityPattern;
+class BlockSparsityPattern;
template <typename number> class SparseMatrix;
template <typename number, int rows, int columns> class BlockSparseMatrix;
class BlockIndices;
* condenses square block sparsity
* patterns.
*/
- template <int blocks>
- void condense (BlockSparsityPattern<blocks,blocks> &sparsity) const;
+ void condense (BlockSparsityPattern &sparsity) const;
/**
* Condense a given matrix. The associated
-template <int blocks>
-void ConstraintMatrix::condense (BlockSparsityPattern<blocks,blocks> &sparsity) const
-{
- Assert (sorted == true, ExcMatrixNotClosed());
- Assert (sparsity.is_compressed() == false, ExcMatrixIsClosed());
- Assert (sparsity.n_rows() == sparsity.n_cols(),
- ExcMatrixNotSquare());
- Assert (sparsity.get_column_indices() == sparsity.get_row_indices(),
- ExcMatrixNotSquare());
-
- const BlockIndices &
- index_mapping = sparsity.get_column_indices();
-
- // store for each index whether it
- // must be distributed or not. If entry
- // is -1, no distribution is necessary.
- // otherwise, the number states which
- // line in the constraint matrix handles
- // this index
- vector<int> distribute (sparsity.n_rows(), -1);
-
- for (unsigned int c=0; c<lines.size(); ++c)
- distribute[lines[c].line] = static_cast<signed int>(c);
-
- const unsigned int n_rows = sparsity.n_rows();
- for (unsigned int row=0; row<n_rows; ++row)
- {
- // get index of this row
- // within the blocks
- const pair<unsigned int,unsigned int>
- block_index = index_mapping.global_to_local(row);
- const unsigned int block_row = block_index.first;
-
- if (distribute[row] == -1)
- // regular line. loop over
- // all columns and see
- // whether this column must
- // be distributed
- {
-
- // to loop over all entries
- // in this row, we have to
- // loop over all blocks in
- // this blockrow and the
- // corresponding row
- // therein
- for (unsigned int block_col=0; block_col<blocks; ++block_col)
- {
- const SparsityPattern &
- block_sparsity = sparsity.block(block_row, block_col);
-
- const unsigned int
- first = block_sparsity.get_rowstart_indices()[block_index.second],
- last = block_sparsity.get_rowstart_indices()[block_index.second+1];
- for (unsigned int j=first; j<last; ++j)
- // end of row reached?
- if (block_sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
- {
- // nothing more
- // to do
- break;
- }
- else
- {
- const unsigned int global_col
- = index_mapping.local_to_global(block_col,
- block_sparsity.get_column_numbers()[j]);
-
- if (distribute[global_col] != -1)
- // distribute entry at regular
- // row @p{row} and irregular column
- // global_col
- {
- for (unsigned int q=0;
- q!=lines[distribute[global_col]]
- .entries.size(); ++q)
- sparsity.add (row,
- lines[distribute[global_col]].entries[q].first);
- };
- };
- };
- }
- else
- {
- // row must be
- // distributed. split the
- // whole row into the
- // chunks defined by the
- // blocks
- for (unsigned int block_col=0; block_col<blocks; ++block_col)
- {
- const SparsityPattern &
- block_sparsity = sparsity.block(block_row,block_col);
-
- const unsigned int
- first = block_sparsity.get_rowstart_indices()[block_index.second],
- last = block_sparsity.get_rowstart_indices()[block_index.second+1];
-
- for (unsigned int j=first; j<last; ++j)
- // end of row reached?
- if (block_sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
- {
- // nothing more to do
- break;
- }
- else
- {
- const unsigned int global_col
- = index_mapping.local_to_global (block_col,
- block_sparsity.get_column_numbers()[j]);
-
- if (distribute[global_col] == -1)
- // distribute entry at irregular
- // row @p{row} and regular column
- // global_col.
- {
- for (unsigned int q=0; q!=lines[distribute[row]].entries.size(); ++q)
- sparsity.add (lines[distribute[row]].entries[q].first,
- global_col);
- }
- else
- // distribute entry at irregular
- // row @p{row} and irregular column
- // @p{global_col}
- {
- for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
- for (unsigned int q=0; q!=lines[distribute[global_col]].entries.size(); ++q)
- sparsity.add (lines[distribute[row]].entries[p].first,
- lines[distribute[global_col]].entries[q].first);
- };
- };
- };
- };
- };
-
- sparsity.compress();
-};
-
-
template<typename number>
void
void
ConstraintMatrix::condense (BlockSparseMatrix<number,blocks,blocks> &uncondensed) const
{
- const BlockSparsityPattern<blocks,blocks> &
+ const BlockSparsityPattern &
sparsity = uncondensed.get_sparsity_pattern ();
Assert (sorted == true, ExcMatrixNotClosed());
Assert (sparsity.is_compressed() == true, ExcMatrixNotClosed());
Assert (sparsity.n_rows() == sparsity.n_cols(),
ExcMatrixNotSquare());
+ Assert (sparsity.n_block_rows() == sparsity.n_block_cols(),
+ ExcMatrixNotSquare());
Assert (sparsity.get_column_indices() == sparsity.get_row_indices(),
ExcMatrixNotSquare());
+void ConstraintMatrix::condense (BlockSparsityPattern &sparsity) const
+{
+ Assert (sorted == true, ExcMatrixNotClosed());
+ Assert (sparsity.is_compressed() == false, ExcMatrixIsClosed());
+ Assert (sparsity.n_rows() == sparsity.n_cols(),
+ ExcMatrixNotSquare());
+ Assert (sparsity.n_block_rows() == sparsity.n_block_cols(),
+ ExcMatrixNotSquare());
+ Assert (sparsity.get_column_indices() == sparsity.get_row_indices(),
+ ExcMatrixNotSquare());
+
+ const BlockIndices &
+ index_mapping = sparsity.get_column_indices();
+
+ const unsigned int n_blocks = sparsity.n_block_rows();
+
+ // store for each index whether it
+ // must be distributed or not. If entry
+ // is -1, no distribution is necessary.
+ // otherwise, the number states which
+ // line in the constraint matrix handles
+ // this index
+ vector<int> distribute (sparsity.n_rows(), -1);
+
+ for (unsigned int c=0; c<lines.size(); ++c)
+ distribute[lines[c].line] = static_cast<signed int>(c);
+
+ const unsigned int n_rows = sparsity.n_rows();
+ for (unsigned int row=0; row<n_rows; ++row)
+ {
+ // get index of this row
+ // within the blocks
+ const pair<unsigned int,unsigned int>
+ block_index = index_mapping.global_to_local(row);
+ const unsigned int block_row = block_index.first;
+
+ if (distribute[row] == -1)
+ // regular line. loop over
+ // all columns and see
+ // whether this column must
+ // be distributed
+ {
+
+ // to loop over all entries
+ // in this row, we have to
+ // loop over all blocks in
+ // this blockrow and the
+ // corresponding row
+ // therein
+ for (unsigned int block_col=0; block_col<n_blocks; ++block_col)
+ {
+ const SparsityPattern &
+ block_sparsity = sparsity.block(block_row, block_col);
+
+ const unsigned int
+ first = block_sparsity.get_rowstart_indices()[block_index.second],
+ last = block_sparsity.get_rowstart_indices()[block_index.second+1];
+ for (unsigned int j=first; j<last; ++j)
+ // end of row reached?
+ if (block_sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
+ {
+ // nothing more
+ // to do
+ break;
+ }
+ else
+ {
+ const unsigned int global_col
+ = index_mapping.local_to_global(block_col,
+ block_sparsity.get_column_numbers()[j]);
+
+ if (distribute[global_col] != -1)
+ // distribute entry at regular
+ // row @p{row} and irregular column
+ // global_col
+ {
+ for (unsigned int q=0;
+ q!=lines[distribute[global_col]]
+ .entries.size(); ++q)
+ sparsity.add (row,
+ lines[distribute[global_col]].entries[q].first);
+ };
+ };
+ };
+ }
+ else
+ {
+ // row must be
+ // distributed. split the
+ // whole row into the
+ // chunks defined by the
+ // blocks
+ for (unsigned int block_col=0; block_col<n_blocks; ++block_col)
+ {
+ const SparsityPattern &
+ block_sparsity = sparsity.block(block_row,block_col);
+
+ const unsigned int
+ first = block_sparsity.get_rowstart_indices()[block_index.second],
+ last = block_sparsity.get_rowstart_indices()[block_index.second+1];
+
+ for (unsigned int j=first; j<last; ++j)
+ // end of row reached?
+ if (block_sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
+ {
+ // nothing more to do
+ break;
+ }
+ else
+ {
+ const unsigned int global_col
+ = index_mapping.local_to_global (block_col,
+ block_sparsity.get_column_numbers()[j]);
+
+ if (distribute[global_col] == -1)
+ // distribute entry at irregular
+ // row @p{row} and regular column
+ // global_col.
+ {
+ for (unsigned int q=0; q!=lines[distribute[row]].entries.size(); ++q)
+ sparsity.add (lines[distribute[row]].entries[q].first,
+ global_col);
+ }
+ else
+ // distribute entry at irregular
+ // row @p{row} and irregular column
+ // @p{global_col}
+ {
+ for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
+ for (unsigned int q=0; q!=lines[distribute[global_col]].entries.size(); ++q)
+ sparsity.add (lines[distribute[row]].entries[p].first,
+ lines[distribute[global_col]].entries[q].first);
+ };
+ };
+ };
+ };
+ };
+
+ sparsity.compress();
+};
+
+
+
unsigned int ConstraintMatrix::n_constraints () const
{
return lines.size();
#undef MatrixType
-#define MatrixType BlockSparsityPattern<2,2>
-matrix_functions_2;
-#undef MatrixType
-
-#define MatrixType BlockSparsityPattern<3,3>
-matrix_functions_2;
-#undef MatrixType
-
template void
DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
- BlockSparsityPattern<2,2> &sparsity);
-template void
-DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
- BlockSparsityPattern<3,3> &sparsity);
+ BlockSparsityPattern &sparsity);
template void
DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
template void
DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
const vector<vector<bool> > &mask,
- BlockSparsityPattern<2,2> &sparsity);
-
-template void
-DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
- const vector<vector<bool> > &mask,
- BlockSparsityPattern<3,3> &sparsity);
+ BlockSparsityPattern &sparsity);
template
void
map<unsigned int,double>::const_iterator dof = boundary_values.begin(),
endd = boundary_values.end();
const unsigned int n_dofs = matrix.m();
- const BlockSparsityPattern<blocks,blocks> &
+ const BlockSparsityPattern &
sparsity_pattern = matrix.get_sparsity_pattern();
// if a diagonal entry is zero
* size of the vector
*/
BlockIndices (const vector<unsigned int> &n);
+
+ /**
+ * Reinitialize the number of
+ * blocks and assign each block
+ * the same number of elements.
+ */
+ void reinit (const unsigned int n_blocks,
+ const unsigned int n_elements_per_block);
/**
* Reinitialize the number of
+inline
+void
+BlockIndices::reinit (const unsigned int n_blocks,
+ const unsigned int n_elements_per_block)
+{
+ const vector<unsigned int> v(n_blocks, n_elements_per_block);
+ reinit (v);
+};
+
+
+
inline
void
BlockIndices::reinit (const vector<unsigned int> &n)
* long as @p{reinit} is not called
* with a new sparsity structure.
*/
- BlockSparseMatrix (const BlockSparsityPattern<rows,columns> &sparsity);
+ BlockSparseMatrix (const BlockSparsityPattern &sparsity);
* The elements of the matrix are
* set to zero by this function.
*/
- virtual void reinit (const BlockSparsityPattern<rows,columns> &sparsity);
+ virtual void reinit (const BlockSparsityPattern &sparsity);
/**
*/
virtual void clear ();
- /**
- * Return the number of blocks in a
- * column.
- */
- unsigned int n_block_rows () const;
-
- /**
- * Return the number of blocks in a
- * row.
- */
- unsigned int n_block_cols () const;
+ /**
+ * Return the number of blocks in a
+ * column.
+ */
+ unsigned int n_block_rows () const;
+
+ /**
+ * Return the number of blocks in a
+ * row.
+ */
+ unsigned int n_block_cols () const;
/**
* Return whether the object is
* function of objects which
* operate on it.
*/
- const BlockSparsityPattern<rows,columns> &
+ const BlockSparsityPattern &
get_sparsity_pattern () const;
/**
* it using the @p{SmartPointer}
* class.
*/
- SmartPointer<const BlockSparsityPattern<rows,columns> > sparsity_pattern;
+ SmartPointer<const BlockSparsityPattern> sparsity_pattern;
/**
* Array of sub-matrices.
template <typename number, int rows, int columns>
BlockSparseMatrix<number,rows,columns>::
-BlockSparseMatrix (const BlockSparsityPattern<rows,columns> &sparsity)
+BlockSparseMatrix (const BlockSparsityPattern &sparsity)
{
reinit (sparsity);
};
template <typename number, int rows, int columns>
void
BlockSparseMatrix<number,rows,columns>::
-reinit (const BlockSparsityPattern<rows,columns> &sparsity)
+reinit (const BlockSparsityPattern &sparsity)
{
sparsity_pattern = &sparsity;
template <typename number, int rows, int columns>
-const BlockSparsityPattern<rows,columns> &
+const BlockSparsityPattern &
BlockSparseMatrix<number,rows,columns>::get_sparsity_pattern () const
{
return *sparsity_pattern;
#include <base/exceptions.h>
#include <base/subscriptor.h>
+#include <base/smartpointer.h>
#include <lac/sparsity_pattern.h>
#include <lac/block_indices.h>
* calls to its member functions to calls to the respective member
* functions of the member sparsity patterns.
*
- * The largest difference between the @p{SparsityPattern} class and this
- * class is probably the absence of several of the constructors as
- * well as the absenceof the @p{reinit} functions. The reason is that
- * mostly, the matrices have different properties and you will want to
- * initialize the blocks making up the matrix separately. You can
- * access the different blocks using the @p{block(row,col)} function.
+ * The largest difference between the @p{SparsityPattern} class and
+ * this class is that mostly, the matrices have different properties
+ * and you will want to work on the blocks making up the matrix rather
+ * than the whole matrix. You can access the different blocks using
+ * the @p{block(row,col)} function.
*
* Attention: this object is not automatically notified if the size of
* one of its subobjects' size is changed. After you initialize the
*
* @author Wolfgang Bangerth, 2000
*/
-template <int rows, int columns=rows>
class BlockSparsityPattern : public Subscriptor
{
public:
BlockSparsityPattern ();
/**
- * Copy operator. For this the same holds
- * as for the copy constructor: it is
- * declared, defined and fine to be called,
- * but the latter only for empty objects.
+ * Initialize the matrix with the
+ * given number of block rows and
+ * columns. The blocks themselves
+ * are still empty, and you have
+ * to call @p{collect_args} after
+ * you assign them sizes.
+ */
+ BlockSparsityPattern (const unsigned int n_rows,
+ const unsigned int n_columns);
+
+ /**
+ * Copy constructor. This
+ * constructor is only allowed to
+ * be called if the sparsity pattern to be
+ * copied is empty, i.e. there
+ * are no block allocated at
+ * present. This is for the same
+ * reason as for the
+ * @p{SparsityPattern}, see there
+ * for the details.
+ */
+ BlockSparsityPattern (const BlockSparsityPattern &bsp);
+
+ /**
+ * Destructor.
+ */
+ ~BlockSparsityPattern ();
+
+ /**
+ * Resize the matrix. This
+ * deletes all previously
+ * existant blocks and replaces
+ * them by unitialized ones,
+ * i.e. ones for which also the
+ * sizes are not yet set. You
+ * have to do that by calling the
+ * @p{reinit} functions of the
+ * blocks themselves. Do not
+ * forget to call
+ * @p{collect_size} after that on
+ * this object.
+ *
+ * The reason that you have to
+ * set sizes of the blocks
+ * yourself is that the sizes may
+ * be varying, the maximum number
+ * of elements per row may be
+ * varying, etc. It is simpler
+ * not to reproduce the interface
+ * of the @p{SparsityPattern}
+ * class here but rather let the
+ * user call whatever function
+ * she desires.
+ */
+ void reinit (const unsigned int n_rows,
+ const unsigned int n_columns);
+
+ /**
+ * Copy operator. For this the
+ * same holds as for the copy
+ * constructor: it is declared,
+ * defined and fine to be called,
+ * but the latter only for empty
+ * objects.
*/
BlockSparsityPattern & operator = (const BlockSparsityPattern &);
*/
void compress ();
- /**
- * Return the number of blocks in a
- * column.
- */
- unsigned int n_block_rows () const;
-
- /**
- * Return the number of blocks in a
- * row.
- */
- unsigned int n_block_cols () const;
+ /**
+ * Return the number of blocks in a
+ * column.
+ */
+ unsigned int n_block_rows () const;
+
+ /**
+ * Return the number of blocks in a
+ * row.
+ */
+ unsigned int n_block_cols () const;
/**
* Return whether the object is
* and then relays to that block.
*/
void add (const unsigned int i, const unsigned int j);
-
+
/**
* Return number of rows of this
* matrix, which equals the
int, int, int, int,
<< "The blocks [" << arg1 << ',' << arg2 << "] and ["
<< arg3 << ',' << arg4 << "] have differing column numbers.");
-
+ /**
+ * Exception
+ */
+ DeclException2 (ExcIncompatibleSizes,
+ int, int,
+ << "The number of blocks " << arg1 << " and " << arg2
+ << " are different.");
private:
+
+ /**
+ * Number of block rows.
+ */
+ unsigned int rows;
+
+ /**
+ * Number of block columns.
+ */
+ unsigned int columns;
+
/**
* Array of sparsity patterns.
*/
- SparsityPattern sub_objects[rows][columns];
+ vector<vector<SmartPointer<SparsityPattern> > > sub_objects;
/**
* Object storing and managing
-template <int rows, int columns>
inline
SparsityPattern &
-BlockSparsityPattern<rows,columns>::block (const unsigned int row,
- const unsigned int column)
+BlockSparsityPattern::block (const unsigned int row,
+ const unsigned int column)
{
- return sub_objects[row][column];
+ return *sub_objects[row][column];
};
-template <int rows, int columns>
inline
const SparsityPattern &
-BlockSparsityPattern<rows,columns>::block (const unsigned int row,
- const unsigned int column) const
+BlockSparsityPattern::block (const unsigned int row,
+ const unsigned int column) const
{
- return sub_objects[row][column];
+ return *sub_objects[row][column];
};
-template <int rows, int columns>
inline
const BlockIndices &
-BlockSparsityPattern<rows,columns>::get_row_indices () const
+BlockSparsityPattern::get_row_indices () const
{
return row_indices;
};
-template <int rows, int columns>
inline
const BlockIndices &
-BlockSparsityPattern<rows,columns>::get_column_indices () const
+BlockSparsityPattern::get_column_indices () const
{
return column_indices;
};
-template <int rows, int columns>
inline
void
-BlockSparsityPattern<rows,columns>::add (const unsigned int i,
- const unsigned int j)
+BlockSparsityPattern::add (const unsigned int i,
+ const unsigned int j)
{
// if you get an error here, are
// you sure you called
const pair<unsigned int,unsigned int>
row_index = row_indices.global_to_local (i),
col_index = column_indices.global_to_local (j);
- sub_objects[row_index.first][col_index.first].add (row_index.second,
- col_index.second);
+ sub_objects[row_index.first][col_index.first]->add (row_index.second,
+ col_index.second);
};
-template <int rows, int columns>
inline
unsigned int
-BlockSparsityPattern<rows,columns>::n_block_cols () const
+BlockSparsityPattern::n_block_cols () const
{
return columns;
}
-template <int rows, int columns>
inline
unsigned int
-BlockSparsityPattern<rows,columns>::n_block_rows () const
+BlockSparsityPattern::n_block_rows () const
{
return rows;
}
+++ /dev/null
-//---------------------------- block_sparsity_pattern.templates.h ---------------------------
-// $Id$
-// Version: $Name$
-//
-// Copyright (C) 1998, 1999, 2000 by the deal.II authors
-//
-// This file is subject to QPL and may not be distributed
-// without copyright and license information. Please refer
-// to the file deal.II/doc/license.html for the text and
-// further information on this license.
-//
-//---------------------------- block_sparsity_pattern.templates.h ---------------------------
-#ifndef __deal2__block_sparsity_pattern_templates_h
-#define __deal2__block_sparsity_pattern_templates_h
-
-
-#include <lac/block_sparsity_pattern.h>
-
-
-template <int rows, int columns>
-BlockSparsityPattern<rows,columns>::BlockSparsityPattern ()
- : row_indices (rows), column_indices(columns)
-{
- Assert (rows>0, ExcInvalidSize (rows));
- Assert (columns>0, ExcInvalidSize (columns));
-};
-
-
-
-template <int rows, int columns>
-BlockSparsityPattern<rows,columns> &
-BlockSparsityPattern<rows,columns>::operator = (const BlockSparsityPattern<rows,columns> &bsp)
-{
- // copy objects
- for (unsigned int i=0; i<rows; ++i)
- for (unsigned int j=0; j<columns; ++j)
- sub_objects[i][j] = bsp.sub_objects[i][j];
- // update index objects
- collect_sizes ();
-
- return *this;
-};
-
-
-
-
-template <int rows, int columns>
-void
-BlockSparsityPattern<rows,columns>::collect_sizes ()
-{
- vector<unsigned int> row_sizes (rows);
- vector<unsigned int> col_sizes (columns);
-
- // first find out the row sizes
- // from the first block column
- for (unsigned int r=0; r<rows; ++r)
- row_sizes[r] = sub_objects[r][0].n_rows();
- // then check that the following
- // block columns have the same
- // sizes
- for (unsigned int c=1; c<columns; ++c)
- for (unsigned int r=0; r<rows; ++r)
- Assert (row_sizes[r] == sub_objects[r][c].n_rows(),
- ExcIncompatibleRowNumbers (r,0,r,c));
-
- // finally initialize the row
- // indices with this array
- row_indices.reinit (row_sizes);
-
-
- // then do the same with the columns
- for (unsigned int c=0; c<columns; ++c)
- col_sizes[c] = sub_objects[0][c].n_cols();
- for (unsigned int r=1; r<rows; ++r)
- for (unsigned int c=0; c<columns; ++c)
- Assert (col_sizes[c] == sub_objects[r][c].n_cols(),
- ExcIncompatibleRowNumbers (0,c,r,c));
-
- // finally initialize the row
- // indices with this array
- column_indices.reinit (col_sizes);
-};
-
-
-
-template <int rows, int columns>
-void
-BlockSparsityPattern<rows,columns>::compress ()
-{
- for (unsigned int i=0; i<rows; ++i)
- for (unsigned int j=0; j<columns; ++j)
- sub_objects[i][j].compress ();
-};
-
-
-
-template <int rows, int columns>
-bool
-BlockSparsityPattern<rows,columns>::empty () const
-{
- for (unsigned int i=0; i<rows; ++i)
- for (unsigned int j=0; j<columns; ++j)
- if (sub_objects[i][j].empty () == false)
- return false;
- return true;
-};
-
-
-
-template <int rows, int columns>
-unsigned int
-BlockSparsityPattern<rows,columns>::max_entries_per_row () const
-{
- unsigned int max_entries = 0;
- for (unsigned int block_row=0; block_row<rows; ++block_row)
- {
- unsigned int this_row = 0;
- for (unsigned int c=0; c<columns; ++c)
- this_row += sub_objects[block_row][c].max_entries_per_row ();
-
- if (this_row > max_entries)
- max_entries = this_row;
- };
- return max_entries;
-};
-
-
-template <int rows, int columns>
-unsigned int
-BlockSparsityPattern<rows,columns>::n_rows () const
-{
- // only count in first column, since
- // all rows should be equivalent
- unsigned int count = 0;
- for (unsigned int r=0; r<rows; ++r)
- count += sub_objects[r][0].n_rows();
- return count;
-};
-
-
-
-template <int rows, int columns>
-unsigned int
-BlockSparsityPattern<rows,columns>::n_cols () const
-{
- // only count in first row, since
- // all rows should be equivalent
- unsigned int count = 0;
- for (unsigned int c=0; c<columns; ++c)
- count += sub_objects[0][c].n_cols();
- return count;
-};
-
-
-
-
-
-template <int rows, int columns>
-unsigned int
-BlockSparsityPattern<rows,columns>::n_nonzero_elements () const
-{
- unsigned int count = 0;
- for (unsigned int i=0; i<rows; ++i)
- for (unsigned int j=0; j<columns; ++j)
- count += sub_objects[i][j].n_nonzero_elements ();
- return count;
-};
-
-
-
-template <int rows, int columns>
-bool
-BlockSparsityPattern<rows,columns>::is_compressed () const
-{
- for (unsigned int i=0; i<rows; ++i)
- for (unsigned int j=0; j<columns; ++j)
- if (sub_objects[i][j].is_compressed () == false)
- return false;
- return true;
-};
-
-
-
-#endif // ifdef block_sparsity_pattern_templates_h
*/
~BlockVector ();
- /**
- * Reinitialize the BlockVector to
- * contain @p{num_blocks} blocks of
- * size @p{block_size} each.
- *
- * If @p{fast==false}, the vector
- * is filled with zeros.
- */
- void reinit (const unsigned int num_blocks,
- const unsigned int block_size,
- const bool fast = false);
+ /**
+ * Reinitialize the BlockVector to
+ * contain @p{num_blocks} blocks of
+ * size @p{block_size} each.
+ *
+ * If @p{fast==false}, the vector
+ * is filled with zeros.
+ */
+ void reinit (const unsigned int num_blocks,
+ const unsigned int block_size,
+ const bool fast = false);
/**
* Reinitialize the BlockVector
//
//---------------------------- block_block_sparsity_pattern.cc ---------------------------
-#include <lac/block_sparsity_pattern.templates.h>
+#include <lac/block_sparsity_pattern.h>
-// explicit instantiations
-template class BlockSparsityPattern<2,2>;
-template class BlockSparsityPattern<2,3>;
-template class BlockSparsityPattern<3,2>;
-template class BlockSparsityPattern<3,3>;
+
+BlockSparsityPattern::BlockSparsityPattern ()
+ :
+ rows (0),
+ columns (0)
+{};
+
+
+
+BlockSparsityPattern::BlockSparsityPattern (const unsigned int r,
+ const unsigned int c)
+ :
+ rows (0),
+ columns (0)
+{
+ reinit (r,c);
+};
+
+
+
+BlockSparsityPattern::~BlockSparsityPattern ()
+{
+ // clear all memory
+ reinit (0,0);
+};
+
+
+
+void
+BlockSparsityPattern::reinit (const unsigned int r,
+ const unsigned int c)
+{
+ // delete previous content
+ for (unsigned int i=0; i<rows; ++i)
+ for (unsigned int j=0; j<columns; ++j)
+ {
+ SparsityPattern * sp = sub_objects[i][j];
+ sub_objects[i][j] = 0;
+ delete sp;
+ };
+ sub_objects.clear ();
+
+ // set new sizes
+ rows = r;
+ columns = c;
+ sub_objects.resize (rows, vector<SmartPointer<SparsityPattern> > (columns));
+
+ // allocate new objects
+ for (unsigned int i=0; i<rows; ++i)
+ for (unsigned int j=0; j<columns; ++j)
+ sub_objects[i][j] = new SparsityPattern;
+};
+
+
+
+BlockSparsityPattern &
+BlockSparsityPattern::operator = (const BlockSparsityPattern &bsp)
+{
+ Assert (rows == bsp.rows, ExcIncompatibleSizes(rows, bsp.rows));
+ Assert (columns == bsp.columns, ExcIncompatibleSizes(columns, bsp.columns));
+ // copy objects
+ for (unsigned int i=0; i<rows; ++i)
+ for (unsigned int j=0; j<columns; ++j)
+ *sub_objects[i][j] = *bsp.sub_objects[i][j];
+ // update index objects
+ collect_sizes ();
+
+ return *this;
+};
+
+
+
+
+void
+BlockSparsityPattern::collect_sizes ()
+{
+ vector<unsigned int> row_sizes (rows);
+ vector<unsigned int> col_sizes (columns);
+
+ // first find out the row sizes
+ // from the first block column
+ for (unsigned int r=0; r<rows; ++r)
+ row_sizes[r] = sub_objects[r][0]->n_rows();
+ // then check that the following
+ // block columns have the same
+ // sizes
+ for (unsigned int c=1; c<columns; ++c)
+ for (unsigned int r=0; r<rows; ++r)
+ Assert (row_sizes[r] == sub_objects[r][c]->n_rows(),
+ ExcIncompatibleRowNumbers (r,0,r,c));
+
+ // finally initialize the row
+ // indices with this array
+ row_indices.reinit (row_sizes);
+
+
+ // then do the same with the columns
+ for (unsigned int c=0; c<columns; ++c)
+ col_sizes[c] = sub_objects[0][c]->n_cols();
+ for (unsigned int r=1; r<rows; ++r)
+ for (unsigned int c=0; c<columns; ++c)
+ Assert (col_sizes[c] == sub_objects[r][c]->n_cols(),
+ ExcIncompatibleRowNumbers (0,c,r,c));
+
+ // finally initialize the row
+ // indices with this array
+ column_indices.reinit (col_sizes);
+};
+
+
+
+void
+BlockSparsityPattern::compress ()
+{
+ for (unsigned int i=0; i<rows; ++i)
+ for (unsigned int j=0; j<columns; ++j)
+ sub_objects[i][j]->compress ();
+};
+
+
+
+bool
+BlockSparsityPattern::empty () const
+{
+ for (unsigned int i=0; i<rows; ++i)
+ for (unsigned int j=0; j<columns; ++j)
+ if (sub_objects[i][j]->empty () == false)
+ return false;
+ return true;
+};
+
+
+
+unsigned int
+BlockSparsityPattern::max_entries_per_row () const
+{
+ unsigned int max_entries = 0;
+ for (unsigned int block_row=0; block_row<rows; ++block_row)
+ {
+ unsigned int this_row = 0;
+ for (unsigned int c=0; c<columns; ++c)
+ this_row += sub_objects[block_row][c]->max_entries_per_row ();
+
+ if (this_row > max_entries)
+ max_entries = this_row;
+ };
+ return max_entries;
+};
+
+
+unsigned int
+BlockSparsityPattern::n_rows () const
+{
+ // only count in first column, since
+ // all rows should be equivalent
+ unsigned int count = 0;
+ for (unsigned int r=0; r<rows; ++r)
+ count += sub_objects[r][0]->n_rows();
+ return count;
+};
+
+
+
+unsigned int
+BlockSparsityPattern::n_cols () const
+{
+ // only count in first row, since
+ // all rows should be equivalent
+ unsigned int count = 0;
+ for (unsigned int c=0; c<columns; ++c)
+ count += sub_objects[0][c]->n_cols();
+ return count;
+};
+
+
+
+
+
+unsigned int
+BlockSparsityPattern::n_nonzero_elements () const
+{
+ unsigned int count = 0;
+ for (unsigned int i=0; i<rows; ++i)
+ for (unsigned int j=0; j<columns; ++j)
+ count += sub_objects[i][j]->n_nonzero_elements ();
+ return count;
+};
+
+
+
+bool
+BlockSparsityPattern::is_compressed () const
+{
+ for (unsigned int i=0; i<rows; ++i)
+ for (unsigned int j=0; j<columns; ++j)
+ if (sub_objects[i][j]->is_compressed () == false)
+ return false;
+ return true;
+};
+
DoFHandler<2> dof_handler;
Sparsity sparsity_pattern;
- Matrix system_matrix;
+ Matrix system_matrix;
Vector system_rhs;
};
template <class Vector, class Matrix, class Sparsity>
LaplaceProblem<Vector,Matrix,Sparsity>::LaplaceProblem () :
dof_handler (triangulation)
+{
+ sparsity_pattern.reinit (system_matrix.n_block_rows(),
+ system_matrix.n_block_cols());
+};
+
+
+
+template <>
+LaplaceProblem<Vector<float>,SparseMatrix<float>,SparsityPattern>::LaplaceProblem () :
+ dof_handler (triangulation)
{};
+
+template <>
+LaplaceProblem<Vector<double>,SparseMatrix<double>,SparsityPattern>::LaplaceProblem () :
+ dof_handler (triangulation)
+{};
+
+
+
template <class Vector, class Matrix, class Sparsity>
void LaplaceProblem<Vector,Matrix,Sparsity>::make_grid_and_dofs ()
{
template <>
-void LaplaceProblem<BlockVector<double>,BlockSparseMatrix<double,2,2>,BlockSparsityPattern<2,2> >::reinit_sparsity ()
+void LaplaceProblem<BlockVector<double>,BlockSparseMatrix<double,2,2>,BlockSparsityPattern>::reinit_sparsity ()
{
const unsigned int n_dofs = dof_handler.n_dofs();
const unsigned int block_size[2] = { n_dofs/3, n_dofs - n_dofs/3 };
template <>
-void LaplaceProblem<BlockVector<double>,BlockSparseMatrix<double,2,2>,BlockSparsityPattern<2,2> >::reinit_vectors ()
+void LaplaceProblem<BlockVector<double>,BlockSparseMatrix<double,2,2>,BlockSparsityPattern>::reinit_vectors ()
{
const unsigned int n_dofs = dof_handler.n_dofs();
const unsigned int block_size_[2] = { n_dofs/3, n_dofs - n_dofs/3 };
template <>
-void LaplaceProblem<BlockVector<double>,BlockSparseMatrix<double,3,3>,BlockSparsityPattern<3,3> >::reinit_sparsity ()
+void LaplaceProblem<BlockVector<double>,BlockSparseMatrix<double,3,3>,BlockSparsityPattern>::reinit_sparsity ()
{
const unsigned int n_dofs = dof_handler.n_dofs();
const unsigned int block_size[3] = { n_dofs/5, n_dofs/7, n_dofs - n_dofs/5 - n_dofs/7 };
template <>
-void LaplaceProblem<BlockVector<double>,BlockSparseMatrix<double,3,3>,BlockSparsityPattern<3,3> >::reinit_vectors ()
+void LaplaceProblem<BlockVector<double>,BlockSparseMatrix<double,3,3>,BlockSparsityPattern>::reinit_vectors ()
{
const unsigned int n_dofs = dof_handler.n_dofs();
const unsigned int block_size_[3] = { n_dofs/5, n_dofs/7, n_dofs - n_dofs/5 - n_dofs/7 };
if (true)
{
- LaplaceProblem<BlockVector<double>,BlockSparseMatrix<double,2,2>,BlockSparsityPattern<2,2> >
+ LaplaceProblem<BlockVector<double>,BlockSparseMatrix<double,2,2>,BlockSparsityPattern>
laplace_problem;
laplace_problem.run ();
if (true)
{
- LaplaceProblem<BlockVector<double>,BlockSparseMatrix<double,3,3>,BlockSparsityPattern<3,3> >
+ LaplaceProblem<BlockVector<double>,BlockSparseMatrix<double,3,3>,BlockSparsityPattern>
laplace_problem;
laplace_problem.run ();
projects / dealii-svn.git / commitdiff