*/
void condense (SparsityPattern &sparsity) const;
+ /**
+ * Same function as above, but
+ * condenses square block sparsity
+ * patterns.
+ */
+ template <int blocks>
+ void condense (BlockSparsityPattern<blocks,blocks> &sparsity) const;
+
/**
* Condense a given matrix. The associated
* matrix struct should be condensed and
template<typename number>
void condense (SparseMatrix<number> &matrix) const;
+ /**
+ * Same function as above, but
+ * condenses square block sparse
+ * matrices.
+ */
+ template <typename number, int blocks>
+ void condense (BlockSparseMatrix<number,blocks,blocks> &sparsity) const;
+
/**
* Condense the given vector #uncondensed#
* into #condensed#. It is the user's
* responsibility to guarantee that all
* entries of #condensed# be zero!
+ *
+ * The #VectorType# may be a
+ * #Vector<float>#,
+ * #Vector<double>#,
+ * #BlockVector<...>#, or any
+ * other type having the same
+ * interface.
*/
- template<typename number>
- void condense (const Vector<number> &uncondensed,
- Vector<number> &condensed) const;
+ template <class VectorType>
+ void condense (const VectorType &uncondensed,
+ VectorType &condensed) const;
/**
- * Condense the given vector in-place.
+ * Condense the given vector
+ * in-place. The #VectorType# may
+ * be a #Vector<float>#,
+ * #Vector<double>#,
+ * #BlockVector<...>#, or any
+ * other type having the same
+ * interface.
*/
- template<typename number>
- void condense (Vector<number> &vec) const;
+ template <class VectorType>
+ void condense (VectorType &vec) const;
/**
- * Re-distribute the elements of the vector
- * #condensed# to #uncondensed#. It is the user's
- * responsibility to guarantee that all
- * entries of #uncondensed# be zero!
+ * Re-distribute the elements of
+ * the vector #condensed# to
+ * #uncondensed#. It is the
+ * user's responsibility to
+ * guarantee that all entries of
+ * #uncondensed# be zero!
+ *
+ * This function undoes the
+ * action of #condense# somehow,
+ * but it should be noted that it
+ * is not the inverse of
+ * #condense#.
*
- * This function undoes the action of
- * #condense# somehow, but it should be noted
- * that it is not the inverse of #condense#-
+ * The #VectorType# may be a
+ * #Vector<float>#,
+ * #Vector<double>#,
+ * #BlockVector<...>#, or any
+ * other type having the same
+ * interface.
*/
- template<typename number>
- void distribute (const Vector<number> &condensed,
- Vector<number> &uncondensed) const;
+ template <class VectorType>
+ void distribute (const VectorType &condensed,
+ VectorType &uncondensed) const;
/**
- * Re-distribute the elements of the vector
- * in-place.
+ * Re-distribute the elements of
+ * the vector in-place. The
+ * #VectorType# may be a
+ * #Vector<float>#,
+ * #Vector<double>#,
+ * #BlockVector<...>#, or any
+ * other type having the same
+ * interface.
*/
- template<typename number>
- void distribute (Vector<number> &vec) const;
+ template <class VectorType>
+ void distribute (VectorType &vec) const;
/**
- * Delete hanging nodes in a vector.
- * Sets all hanging node values to zero.
+ * Delete hanging nodes in a
+ * vector. Sets all hanging node
+ * values to zero. The
+ * #VectorType# may be a
+ * #Vector<float>#,
+ * #Vector<double>#,
+ * #BlockVector<...>#, or any
+ * other type having the same
+ * interface.
*/
- template<typename number>
- void set_zero (Vector<number> &vec) const;
+ template <class VectorType>
+ void set_zero (VectorType &vec) const;
/**
#include <dofs/dof_constraints.h>
-#include <lac/sparse_matrix.h>
template<typename number>
-template<typename number>
+template <typename number, int blocks>
void
-ConstraintMatrix::condense (const Vector<number> &uncondensed,
- Vector<number> &condensed) const
+ConstraintMatrix::condense (BlockSparseMatrix<number,blocks,blocks> &uncondensed) const
+{
+ const BlockSparsityPattern<blocks,blocks> &
+ 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.get_column_indices() == sparsity.get_row_indices(),
+ ExcMatrixNotSquare());
+
+ const BlockIndices<blocks> &
+ 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)
+ {
+ if (distribute[row] == -1)
+ {
+ // get index of this row
+ // within the blocks
+ const pair<unsigned int,unsigned int>
+ block_index = index_mapping.global_to_local(row);
+
+ // regular line. loop over
+ // all columns and see
+ // whether this column must
+ // be distributed
+ for (unsigned int block_col=0; block_col<blocks; ++block_col)
+ {
+ const SparsityPattern &
+ block_sparsity = sparsity.block(block_index.first,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)
+ {
+ // this should not happen, since
+ // we only operate on compressed
+ // matrices!
+ Assert (false, ExcMatrixNotClosed());
+ 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 #row# and irregular column
+ // global_col;
+ // set old entry to zero
+ {
+ const double old_value =
+ uncondensed.block(block_index.first,block_col).global_entry(j);
+
+ for (unsigned int q=0;
+ q!=lines[distribute[global_col]]
+ .entries.size(); ++q)
+ uncondensed.add (row,
+ lines[distribute[global_col]].entries[q].first,
+ old_value *
+ lines[distribute[global_col]].entries[q].second);
+
+ uncondensed.block(block_index.first,block_col).global_entry(j)
+ = 0.;
+ };
+ };
+ };
+ }
+ else
+ {
+ }
+
+// // row must be distributed
+// for (unsigned int j=sparsity.get_rowstart_indices()[row];
+// j<sparsity.get_rowstart_indices()[row+1]; ++j)
+// // end of row reached?
+// if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
+// {
+// // this should not happen, since
+// // we only operate on compressed
+// // matrices!
+// Assert (false, ExcMatrixNotClosed());
+// break;
+// }
+// else
+// {
+// if (distribute[sparsity.get_column_numbers()[j]] == -1)
+// // distribute entry at irregular
+// // row #row# and regular column
+// // sparsity.get_column_numbers()[j]. set old
+// // entry to zero
+// {
+// for (unsigned int q=0;
+// q!=lines[distribute[row]].entries.size(); ++q)
+// uncondensed.add (lines[distribute[row]].entries[q].first,
+// sparsity.get_column_numbers()[j],
+// uncondensed.global_entry(j) *
+// lines[distribute[row]].entries[q].second);
+
+// uncondensed.global_entry(j) = 0.;
+// }
+// else
+// // distribute entry at irregular
+// // row #row# and irregular column
+// // sparsity.get_column_numbers()[j]
+// // set old entry to one if on main
+// // diagonal, zero otherwise
+// {
+// for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
+// for (unsigned int q=0;
+// q!=lines[distribute[sparsity.get_column_numbers()[j]]]
+// .entries.size(); ++q)
+// uncondensed.add (lines[distribute[row]].entries[p].first,
+// lines[distribute[sparsity.get_column_numbers()[j]]]
+// .entries[q].first,
+// uncondensed.global_entry(j) *
+// lines[distribute[row]].entries[p].second *
+// lines[distribute[sparsity.get_column_numbers()[j]]]
+// .entries[q].second);
+
+// uncondensed.global_entry(j) = (row == sparsity.get_column_numbers()[j] ?
+// 1. : 0. );
+// };
+// };
+ };
+};
+
+
+
+template<class VectorType>
+void
+ConstraintMatrix::condense (const VectorType &uncondensed,
+ VectorType &condensed) const
{
Assert (sorted == true, ExcMatrixNotClosed());
Assert (condensed.size()+n_constraints() == uncondensed.size(),
-template<typename number>
+template <class VectorType>
void
-ConstraintMatrix::condense (Vector<number> &vec) const
+ConstraintMatrix::condense (VectorType &vec) const
{
Assert (sorted == true, ExcMatrixNotClosed());
-template<typename number>
+template <class VectorType>
void
-ConstraintMatrix::set_zero (Vector<number> &vec) const
+ConstraintMatrix::set_zero (VectorType &vec) const
{
Assert (sorted == true, ExcMatrixNotClosed());
-template<typename number>
+template<class VectorType>
void
-ConstraintMatrix::distribute (const Vector<number> &condensed,
- Vector<number> &uncondensed) const
+ConstraintMatrix::distribute (const VectorType &condensed,
+ VectorType &uncondensed) const
{
Assert (sorted == true, ExcMatrixNotClosed());
Assert (condensed.size()+n_constraints() == uncondensed.size(),
-template<typename number>
+template<class VectorType>
void
-ConstraintMatrix::distribute (Vector<number> &vec) const
+ConstraintMatrix::distribute (VectorType &vec) const
{
Assert (sorted == true, ExcMatrixNotClosed());
#include <dofs/dof_constraints.h>
+#include <dofs/dof_constraints.templates.h>
+
#include <lac/sparsity_pattern.h>
#include <lac/vector.h>
+#include <lac/block_vector.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
#include <iostream>
#include <algorithm>
#include <numeric>
};
+
ConstraintMatrix::ConstraintMatrix () :
- lines(), sorted(false) {};
+ lines(),
+ sorted(false)
+{};
-void ConstraintMatrix::add_line (const unsigned int line) {
+void ConstraintMatrix::add_line (const unsigned int line)
+{
Assert (sorted==false, ExcMatrixIsClosed());
// check whether line already exists;
};
+
void ConstraintMatrix::add_entry (const unsigned int line,
const unsigned int column,
- const double value) {
+ const double value)
+{
Assert (sorted==false, ExcMatrixIsClosed());
vector<ConstraintLine>::iterator line_ptr;
};
-void ConstraintMatrix::close () {
+
+void ConstraintMatrix::close ()
+{
Assert (sorted==false, ExcMatrixIsClosed());
// sort the entries in the different lines
};
-void ConstraintMatrix::clear () {
+
+void ConstraintMatrix::clear ()
+{
lines = vector<ConstraintLine>();
sorted = false;
};
+
void ConstraintMatrix::condense (const SparsityPattern &uncondensed,
- SparsityPattern &condensed) const {
+ SparsityPattern &condensed) const
+{
Assert (sorted == true, ExcMatrixNotClosed());
Assert (uncondensed.is_compressed() == true, ExcMatrixNotClosed());
Assert (uncondensed.n_rows() == uncondensed.n_cols(),
};
-void ConstraintMatrix::condense (SparsityPattern &sparsity) const {
+
+void ConstraintMatrix::condense (SparsityPattern &sparsity) const
+{
Assert (sorted == true, ExcMatrixNotClosed());
Assert (sparsity.is_compressed() == false, ExcMatrixIsClosed());
Assert (sparsity.n_rows() == sparsity.n_cols(),
// 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);
+
+ unsigned int n_rows = sparsity.n_rows();
+ for (unsigned int row=0; row<n_rows; ++row)
+ {
+ if (distribute[row] == -1)
+ // regular line. loop over cols
+ for (unsigned int j=sparsity.get_rowstart_indices()[row];
+ j<sparsity.get_rowstart_indices()[row+1]; ++j)
+ {
+ const unsigned int column = sparsity.get_column_numbers()[j];
+
+ // end of row reached?
+ if (column == SparsityPattern::invalid_entry)
+ break;
+ else
+ if (distribute[column] != -1)
+ {
+ // distribute entry at regular
+ // row #row# and irregular column
+ // sparsity.colnums[j]
+ for (unsigned int q=0;
+ q!=lines[distribute[column]].entries.size();
+ ++q)
+ sparsity.add (row,
+ lines[distribute[column]].entries[q].first);
+ };
+ }
+ else
+ // row must be distributed
+ for (unsigned int j=sparsity.get_rowstart_indices()[row];
+ j<sparsity.get_rowstart_indices()[row+1]; ++j)
+ // end of row reached?
+ if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
+ break;
+ else
+ {
+ if (distribute[sparsity.get_column_numbers()[j]] == -1)
+ // distribute entry at irregular
+ // row #row# and regular column
+ // sparsity.colnums[j]
+ for (unsigned int q=0;
+ q!=lines[distribute[row]].entries.size(); ++q)
+ sparsity.add (lines[distribute[row]].entries[q].first,
+ sparsity.get_column_numbers()[j]);
+ else
+ // distribute entry at irregular
+ // row #row# and irregular column
+ // sparsity.get_column_numbers()[j]
+ for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
+ for (unsigned int q=0;
+ q!=lines[distribute[sparsity.get_column_numbers()[j]]]
+ .entries.size(); ++q)
+ sparsity.add (lines[distribute[row]].entries[p].first,
+ lines[distribute[sparsity.get_column_numbers()[j]]]
+ .entries[q].first);
+ };
+ };
+
+ sparsity.compress();
+};
+
+
+
+template <int blocks>
+void ConstraintMatrix::condense (BlockSparsityPattern<blocks,blocks> &sparsity) const
+{
+ Assert (false, ExcInternalError());
+/*
+ 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());
+
+ // 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);
};
sparsity.compress();
+*/
};
+
unsigned int ConstraintMatrix::n_constraints () const
{
return lines.size();
};
+
bool ConstraintMatrix::is_constrained (const unsigned int index) const
{
if (sorted == true)
AssertThrow (out, ExcIO());
};
-#include <dofs/dof_constraints.templates.h>
-#define number double
-template void ConstraintMatrix::condense(const SparseMatrix<number> &uncondensed,
- SparseMatrix<number> &condensed) const;
-template void ConstraintMatrix::condense(SparseMatrix<number> &uncondensed) const;
-template void ConstraintMatrix::condense(const Vector<number> &uncondensed,
- Vector<number> &condensed) const;
-template void ConstraintMatrix::condense(Vector<number> &vec) const;
-template void ConstraintMatrix::set_zero(Vector<number> &vec) const;
-template void ConstraintMatrix::distribute(const Vector<number> &condensed,
- Vector<number> &uncondensed) const;
-template void ConstraintMatrix::distribute(Vector<number> &vec) const;
-
-#undef number
-#define number float
-
-template void ConstraintMatrix::condense(const SparseMatrix<number> &uncondensed,
- SparseMatrix<number> &condensed) const;
-template void ConstraintMatrix::condense(SparseMatrix<number> &uncondensed) const;
-template void ConstraintMatrix::condense(const Vector<number> &uncondensed,
- Vector<number> &condensed) const;
-template void ConstraintMatrix::condense(Vector<number> &vec) const;
-template void ConstraintMatrix::set_zero(Vector<number> &vec) const;
-template void ConstraintMatrix::distribute(const Vector<number> &condensed,
- Vector<number> &uncondensed) const;
-template void ConstraintMatrix::distribute(Vector<number> &vec) const;
-
-#undef number
+
+
+
+// explicit instantiations
+//
+// define a list of functions for vectors and matrices, respectively,
+// where the vector/matrix can be replaced using a preprocessor
+// variable VectorType/MatrixType. note that we cannot do so by using
+// a preprocessor function with one arg, since
+// #vector_functions(BlockVector<2,double>)# is not recognized as one
+// arg, and putting parentheses around the arg yields incorrect
+// syntax...
+
+#define vector_functions \
+ template void ConstraintMatrix::condense(const VectorType &uncondensed,\
+ VectorType &condensed) const;\
+ template void ConstraintMatrix::condense(VectorType &vec) const;\
+ template void ConstraintMatrix::set_zero(VectorType &vec) const;\
+ template void ConstraintMatrix::distribute(const VectorType &condensed,\
+ VectorType &uncondensed) const;\
+ template void ConstraintMatrix::distribute(VectorType &vec) const;
+
+
+
+#define matrix_functions_1 \
+ template void ConstraintMatrix::condense(const MatrixType &uncondensed,\
+ MatrixType &condensed) const;
+#define matrix_functions_2 \
+ template void ConstraintMatrix::condense(MatrixType &uncondensed) const;
+
+
+
+#define VectorType Vector<float>
+vector_functions;
+#undef VectorType
+
+#define VectorType Vector<double>
+vector_functions;
+#undef VectorType
+
+#define VectorType BlockVector<2,double>
+vector_functions;
+#undef VectorType
+
+#define VectorType BlockVector<3,double>
+vector_functions;
+#undef VectorType
+
+
+
+#define MatrixType SparseMatrix<float>
+matrix_functions_1;
+matrix_functions_2;
+#undef MatrixType
+
+#define MatrixType SparseMatrix<double>
+matrix_functions_1;
+matrix_functions_2;
+#undef MatrixType
+
+// block sparse matrices are only implemented for one of the two matrix functions
+#define MatrixType BlockSparseMatrix<double,2,2>
+matrix_functions_2;
+#undef MatrixType
+
+#define MatrixType BlockSparseMatrix<double,3,3>
+matrix_functions_2;
+#undef MatrixType
+
+
+#define MatrixType BlockSparsityPattern<2,2>
+matrix_functions_2;
+#undef MatrixType
+
+#define MatrixType BlockSparsityPattern<3,3>
+matrix_functions_2;
+#undef MatrixType
+