#include <base/exceptions.h>
#include <base/subscriptor.h>
+template <typename> class Vector;
+template <typename> class FullMatrix;
class SparsityPattern;
class CompressedSparsityPattern;
class BlockSparsityPattern;
* those which are subject to constraints.
*
* @item Use only one sparsity pattern and one matrix: doing it this way, the
- * condense functions add nonzero entries to the sparsity pattern of the
- * large matrix (with constrained nodes in it) where the condensation process
- * of the matrix will create additional nonzero elements. In the condensation
- * process itself, lines and rows subject to constraints are distributed to
- * the lines and rows of unconstrained nodes. The constrained lines remain in
- * place, however, unlike in the first possibility described above. In order
- * not to disturb the solution process, these lines and rows are filled with
- * zeros and identity on the main diagonal; the appropriate value in the right
- * hand sides is set to zero. This way, the constrained node will always get
- * the value zero upon solution of the equation system and will not couple to
- * other nodes any more.
+ * condense functions add nonzero entries to the sparsity pattern of the large
+ * matrix (with constrained nodes in it) where the condensation process of the
+ * matrix will create additional nonzero elements. In the condensation process
+ * itself, lines and rows subject to constraints are distributed to the lines
+ * and rows of unconstrained nodes. The constrained lines remain in place,
+ * however, unlike in the first possibility described above. In order not to
+ * disturb the solution process, these lines and rows are filled with zeros
+ * and an appropriate positive value on the main diagonal (we choose an
+ * average of the magnitudes of the other diagonal elements, so as to make
+ * sure that the new diagonal entry has the same order of magnitude as the
+ * other entries; this preserves the scaling properties of the matrix). The
+ * appropriate value in the right hand sides is set to zero. This way, the
+ * constrained node will always get the value zero upon solution of the
+ * equation system and will not couple to other nodes any more.
*
* This method has the advantage that only one matrix and sparsity pattern is
* needed thus using less memory. Additionally, the condensation process is
* functions exist in variants for PETSc vectors. However, using them is
* typically expensive, and should be avoided. You should use the same
* techniques as mentioned above to avoid their use.
+ *
*
+ * @sect3{Avoiding explicit condensation}
+ *
+ * Sometimes, one wants to avoid condensation at all. This may be the case
+ * since condensation is an expensive operation, or because no condense()
+ * function is defined for the matrix you use (this is, for example, the case
+ * for the PETSc wrapper classes, where we have no access to the underlying
+ * representation of the matrix, and therefore cannot efficiently implement
+ * the condense() operation). In this case, one possibility is to distribute
+ * local entries to the final destinations right at the moment of transferring
+ * them into the global matrices and vectors. For this, one can use the
+ * distribute_local_to_global() functions of this class, which make a
+ * subsequent call to condense() unnecessary.
*
+ *
* @sect3{Distributing constraints}
*
* After solving the condensed system of equations, the solution vector has to
* for unconstrained nodes, and constrained nodes need to get their values in
* a second step.
*
- * @author Wolfgang Bangerth, 1998
+ * @author Wolfgang Bangerth, 1998, 2004
*/
class ConstraintMatrix : public Subscriptor
{
template <class VectorType>
void set_zero (VectorType &vec) const;
-
+ /**
+ * This function takes a vector of local
+ * contributions (@arg local_vector)
+ * corresponding to the degrees of
+ * freedom indices given in @arg
+ * local_dof_indices and distributes them
+ * to the global vector. In most cases,
+ * these local contributions will be the
+ * result of an integration over a cell
+ * or face of a cell. However, as long as
+ * @arg local_vector and @arg
+ * local_dof_indices have the same number
+ * of elements, this function is happy
+ * with whatever it is given.
+ *
+ * In contrast to the similar function in
+ * the DoFAccessor class, this function
+ * also takes care of constraints,
+ * i.e. of one of the elements of @arg
+ * local_dof_indices belongs to a
+ * constrained node, then rather than
+ * writing the corresponding element of
+ * @arg local_vector into @arg
+ * global_vector, the element is
+ * distributed to the entries in the
+ * global vector to which this particular
+ * degree of freedom is constrained.
+ *
+ * Thus, by using this function to
+ * distribute local contributions to the
+ * global object, one saves the call to
+ * the condense function after the
+ * vectors and matrices are fully
+ * assembled.
+ */
+ template <typename VectorType>
+ void
+ distribute_local_to_global (const Vector<double> &local_vector,
+ const std::vector<unsigned int> &local_dof_indices,
+ VectorType &global_vector) const;
+
+ /**
+ * This function takes a matrix of local
+ * contributions (@arg local_matrix)
+ * corresponding to the degrees of
+ * freedom indices given in @arg
+ * local_dof_indices and distributes them
+ * to the global matrix. In most cases,
+ * these local contributions will be the
+ * result of an integration over a cell
+ * or face of a cell. However, as long as
+ * @arg local_matrix and @arg
+ * local_dof_indices have the same number
+ * of elements, this function is happy
+ * with whatever it is given.
+ *
+ * In contrast to the similar function in
+ * the DoFAccessor class, this function
+ * also takes care of constraints,
+ * i.e. of one of the elements of @arg
+ * local_dof_indices belongs to a
+ * constrained node, then rather than
+ * writing the corresponding element of
+ * @arg local_matrix into @arg
+ * global_matrix, the element is
+ * distributed to the entries in the
+ * global matrix to which this particular
+ * degree of freedom is constrained.
+ *
+ * With this scheme, we never write into
+ * rows or columns of constrained degrees
+ * of freedom. In order to make sure that
+ * the resulting matrix can still be
+ * inverted, we need to do something with
+ * the diagonal elements corresponding to
+ * constrained nodes. Thus, if a degree
+ * of freedom in @arg local_dof_indices
+ * is constrained, we distribute the
+ * corresponding entries in the matrix,
+ * but also add the absolute value of the
+ * diagonal entry of the local matrix to
+ * the corresponding entry in the global
+ * matrix. Since the exact value of the
+ * diagonal element is not important (the
+ * value of the respective degree of
+ * freedom will be overwritten by the
+ * distribute() call later on anyway),
+ * this guarantees that the diagonal
+ * entry is always non-zero, positive,
+ * and of the same order of magnitude as
+ * the other entries of the matrix.
+ *
+ * Thus, by using this function to
+ * distribute local contributions to the
+ * global object, one saves the call to
+ * the condense function after the
+ * vectors and matrices are fully
+ * assembled.
+ */
+ template <typename MatrixType>
+ void
+ distribute_local_to_global (const FullMatrix<double> &local_matrix,
+ const std::vector<unsigned int> &local_dof_indices,
+ MatrixType &global_matrix) const;
+
/**
* Print the constraint lines. Mainly for
* debugging purposes.
// $Id$
// Version: $Name$
//
-// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 by the deal.II authors
+// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
#include <base/config.h>
#include <dofs/dof_constraints.h>
+#include <lac/full_matrix.h>
#include <lac/sparsity_pattern.h>
#include <lac/sparse_matrix.h>
#include <lac/block_sparsity_pattern.h>
+template <typename VectorType>
+void
+ConstraintMatrix::
+distribute_local_to_global (const Vector<double> &local_vector,
+ const std::vector<unsigned int> &local_dof_indices,
+ VectorType &global_vector) const
+{
+ Assert (local_vector.size() == local_dof_indices.size(),
+ ExcWrongDimension());
+ Assert (sorted == true, ExcMatrixNotClosed());
+
+ const unsigned int n_local_dofs = local_vector.size();
+
+ // have a special case where there are no
+ // constraints at all, since then we can be
+ // a lot faster
+ if (lines.size() == 0)
+ {
+ for (unsigned int i=0; i<n_local_dofs; ++i)
+ global_vector(local_dof_indices[i]) += local_vector(i);
+ }
+ else
+ {
+ for (unsigned int i=0; i<n_local_dofs; ++i)
+ {
+ // first figure out whether this
+ // dof is constrained
+ ConstraintLine index_comparison;
+ index_comparison.line = local_dof_indices[i];
+
+ const std::vector<ConstraintLine>::const_iterator
+ position = std::lower_bound (lines.begin(),
+ lines.end(),
+ index_comparison);
+
+ // if the line is not constrained,
+ // then simply copy the
+ // data. otherwise distribute it
+ if (position->line != local_dof_indices[i])
+ global_vector(local_dof_indices[i]) += local_vector(i);
+ else
+ for (unsigned int j=0; j<position->entries.size(); ++j)
+ global_vector(position->entries[j].first)
+ += local_vector(i) * position->entries[j].second;
+ }
+ }
+}
+
+
+
+template <typename MatrixType>
+void
+ConstraintMatrix::
+distribute_local_to_global (const FullMatrix<double> &local_matrix,
+ const std::vector<unsigned int> &local_dof_indices,
+ MatrixType &global_matrix) const
+{
+ Assert (local_matrix.n() == local_dof_indices.size(),
+ ExcWrongDimension());
+ Assert (local_matrix.m() == local_dof_indices.size(),
+ ExcWrongDimension());
+ Assert (sorted == true, ExcMatrixNotClosed());
+
+ const unsigned int n_local_dofs = local_dof_indices.size();
+
+ // have a special case where there are no
+ // constraints at all, since then we can be
+ // a lot faster
+ if (lines.size() == 0)
+ {
+ for (unsigned int i=0; i<n_local_dofs; ++i)
+ for (unsigned int j=0; j<n_local_dofs; ++j)
+ global_matrix.add(local_dof_indices[i],
+ local_dof_indices[j],
+ local_matrix(i,j));
+ }
+ else
+ {
+ // here we have to do something a
+ // little nastier than in the
+ // respective function for vectors. the
+ // reason is that we have two nested
+ // loops and we don't want to
+ // repeatedly check whether a certain
+ // dof is constrained or not by
+ // searching over all the constrained
+ // dofs. so we have to cache this
+ // knowledge, by storing for each dof
+ // index whether and where the line of
+ // the constraint matrix is located
+ std::vector<const ConstraintLine *>
+ constraint_lines (n_local_dofs,
+ static_cast<const ConstraintLine *>(0));
+ for (unsigned int i=0; i<n_local_dofs; ++i)
+ {
+ ConstraintLine index_comparison;
+ index_comparison.line = local_dof_indices[i];
+
+ const std::vector<ConstraintLine>::const_iterator
+ position = std::lower_bound (lines.begin(),
+ lines.end(),
+ index_comparison);
+
+ // if this dof is constrained, then
+ // set the respective entry in the
+ // array. otherwise leave it at the
+ // invalid position
+ if ((position != lines.end()) &&
+ (position->line == local_dof_indices[i]))
+ constraint_lines[i] = &*position;
+ }
+
+
+ // now distribute entries
+ for (unsigned int i=0; i<n_local_dofs; ++i)
+ {
+ const ConstraintLine *position_i = constraint_lines[i];
+ const bool is_constrained_i = (position_i != 0);
+
+ for (unsigned int j=0; j<n_local_dofs; ++j)
+ {
+ const ConstraintLine *position_j = constraint_lines[j];
+ const bool is_constrained_j = (position_j != 0);
+
+ if ((is_constrained_i == false) &&
+ (is_constrained_j == false))
+ {
+ // neither row nor column
+ // is constrained
+ global_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ local_matrix(i,j));
+ }
+ else if ((is_constrained_i == true) &&
+ (is_constrained_j == false))
+ {
+ // ok, row is constrained,
+ // but column is not
+ for (unsigned int q=0; q<position_i->entries.size(); ++q)
+ global_matrix.add (position_i->entries[q].first,
+ local_dof_indices[j],
+ local_matrix(i,j) *
+ position_i->entries[q].second);
+ }
+ else if ((is_constrained_i == false) &&
+ (is_constrained_j == true))
+ {
+ // simply the other way
+ // round: row ok, column is
+ // constrained
+ for (unsigned int q=0; q<position_j->entries.size(); ++q)
+ global_matrix.add (local_dof_indices[i],
+ position_j->entries[q].first,
+ local_matrix(i,j) *
+ position_j->entries[q].second);
+ }
+ else if ((is_constrained_i == true) &&
+ (is_constrained_j == true))
+ {
+ // last case: both row and
+ // column are constrained
+ for (unsigned int p=0; p<position_i->entries.size(); ++p)
+ for (unsigned int q=0; q<position_j->entries.size(); ++q)
+ global_matrix.add (position_i->entries[p].first,
+ position_j->entries[q].first,
+ local_matrix(i,j) *
+ position_i->entries[p].second *
+ position_j->entries[q].second);
+
+ // to make sure that the
+ // global matrix remains
+ // invertible, we need to
+ // do something with the
+ // diagonal elements. add
+ // the absolute value of
+ // the local matrix, so the
+ // resulting entry will
+ // always be positive and
+ // furthermore be in the
+ // same order of magnitude
+ // as the other elements of
+ // the matrix
+ if (i == j)
+ global_matrix.add (local_dof_indices[i],
+ local_dof_indices[i],
+ local_matrix(i,i));
+ }
+ else
+ Assert (false, ExcInternalError());
+ }
+ }
+ }
+}
+
+
+
template<class VectorType>
void
ConstraintMatrix::distribute (const VectorType &condensed,