// $Id$
// Version: $Name$
//
-// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
+// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
/**
- * This class implements dealing with linear homogeneous constraints on
- * degrees of freedom. In particular, it handles constraints of the form
- * $x_{i_1} = \sum_{j=2}^M a_{i_j} x_{i_j}$. In the context of adaptive finite
- * elements, such constraints appear most frequently as "hanging nodes". For
- * example, when using Q1 and Q2 elements (i.e. using
+ * This class implements dealing with linear (possibly inhomogeneous)
+ * constraints on degrees of freedom. In particular, it handles constraints
+ * of the form $x_{i_1} = \sum_{j=2}^M a_{i_j} x_{i_j} + b_i$. In the
+ * context of adaptive finite elements, such constraints appear most
+ * frequently as "hanging nodes" and for implementing Dirichlet boundary
+ * conditions in strong form.
+ *
+ *
+ * <h3>Using the ConstraintMatrix for hanging nodes</h3>
+ *
+ * For example, when using Q1 and Q2 elements (i.e. using
* FE_Q<dim,spacedim>(1) and FE_Q<dim,spacedim>(2)) on the two
* marked cells of the mesh
*
* some detail in the @ref hp_paper "hp paper".
*
*
+ * <h3>Using the ConstraintMatrix for Dirichlet boundary conditions</h3>
+ *
+ * The ConstraintMatrix provides an alternative to implementing Dirichlet
+ * boundary conditions (where the alternative is to use the functions
+ * VectorTools::interpolate_boundary_values and
+ * MatrixTools::apply_boundary_values). The general principle of Dirichlet
+ * conditions are algebraic constraints of the form $x_{i} = b_i$, which
+ * fits into the form as described above.
+ *
+ *
* <h3>Description of constraints</h3>
*
- * Each "line" in objects of this class corresponds to one constrained degree
- * of freedom, with the number of the line being $i_1$, and the entries in
- * this line being pairs $(i_j,a_{i_j})$. Note that the constraints are linear
- * in the $x_i$, and that there is no constant (non-homogeneous) term in the
- * constraint. However, this is exactly the form we need for hanging node and
- * certain other constraints, where we need to constrain one degree of freedom
- * in terms of others. The name of the class stems from the fact that these
- * constraints can be represented in matrix form as $X x = 0$, and this object
- * then describes the matrix $X$. The most frequent way to create/fill objects
- * of this type is using the DoFTools::make_hanging_node_constraints()
- * function. The use of these objects is first explained in @ref step_6
- * "step-6".
+ * Each "line" in objects of this class corresponds to one constrained
+ * degree of freedom, with the number of the line being $i_1$, and the
+ * entries in this line being pairs $(i_j,a_{i_j})$. Note that the
+ * constraints are linear in the $x_i$, and that there might be a constant
+ * (non-homogeneous) term in the constraint. This is exactly the form we
+ * need for hanging node constraints, where we need to constrain one degree
+ * of freedom in terms of others. There are other conditions of this form
+ * possible, for example for implementing mean value conditions as is done
+ * in the @ref step_11 "step-11" tutorial program. The name of the class
+ * stems from the fact that these constraints can be represented in matrix
+ * form as $X x = b$, and this object then describes the matrix $X$. The
+ * most frequent way to create/fill objects of this type is using the
+ * DoFTools::make_hanging_node_constraints() function. The use of these
+ * objects is first explained in @ref step_6 "step-6".
*
* Matrices of the present type are organized in lines (rows), but only those
* lines are stored where constraints are present. New constraints are added
* The condensation functions exist for different argument types. The
* in-place functions (i.e. those following the second way) exist for
* arguments of type SparsityPattern, SparseMatrix and
- * BlockSparseMatrix. Note that there are no versions for arguments of
- * type PETScWrappers::SparseMatrix() or any of the other PETSc or
- * Trilinos matrix wrapper classes. This is due to the fact that it is
- * relatively hard to get a representation of the sparsity structure
- * of PETSc matrices, and to modify them; this holds in particular, if
- * the matrix is actually distributed across a cluster of
- * computers. If you want to use PETSc matrices, you can either copy
- * an already condensed deal.II matrix, or build the PETSc matrix in
- * the already condensed form, see the discussion below.
+ * BlockSparseMatrix. Note that there are no versions for arguments of type
+ * PETScWrappers::SparseMatrix() or any of the other PETSc or Trilinos
+ * matrix wrapper classes. This is due to the fact that it is relatively
+ * hard to get a representation of the sparsity structure of PETSc matrices,
+ * and to modify them; this holds in particular, if the matrix is actually
+ * distributed across a cluster of computers. If you want to use
+ * PETSc/Trilinos matrices, you can either copy an already condensed deal.II
+ * matrix, or build the PETSc/Trilinos matrix in the already condensed form,
+ * see the discussion below.
*
*
* <h5>Condensing vectors</h5>
* object has been condensed, further condensation operations don't change it
* any more.
*
- * In contrast to the matrix condensation functions, the vector
- * condensation functions exist in variants for PETSc and Trilinos
- * vectors. However, using them is typically expensive, and should be
- * avoided. You should use the same techniques as mentioned above to
- * avoid their use.
+ * In contrast to the matrix condensation functions, the vector condensation
+ * functions exist in variants for PETSc and Trilinos vectors. However,
+ * using them is typically expensive, and should be avoided. You should use
+ * the same techniques as mentioned above to avoid their use.
+ *
+ *
+ * <h5>Treatment of inhomogeneous constraints</h5>
+ *
+ * In case some constraint lines have inhomogeneities (which is typically
+ * the case if the constraint comes from implementation of inhomogeneous
+ * boundary conditions), the situation is a bit more complicated. This is
+ * because the elimination of the non-diagonal values in the matrix generate
+ * contributions in the eliminated rows in the vector. This means that
+ * inhomogeneities can only be handled with functions that act
+ * simultaneously on a matrix and a vector. This means that all
+ * inhomogeneities are ignored in case the respective condense function is
+ * called without any matrix (or if the matrix has already been condensed
+ * before).
*
*
* <h3>Avoiding explicit condensation</h3>
* paper". This is the case discussed in the hp tutorial program,
* @ref step_27 "step-27", as well as in @ref step_31 "step-31".
*
- * <li>
- * There may not be a condense() function for the matrix you use
- * (this is, for example, the case for the PETSc and Trilinos wrapper
- * classes, where we have no access to the underlying representation
- * of the matrix, and therefore cannot efficiently implement the
- * condense() operation). This is the case discussed in @ref step_17
- * "step-17" and @ref step_18 "step-18".
+ * <li> There may not be a condense() function for the matrix you use (this
+ * is, for example, the case for the PETSc and Trilinos wrapper classes,
+ * where we have no access to the underlying representation of the matrix,
+ * and therefore cannot efficiently implement the condense()
+ * operation). This is the case discussed in @ref step_17 "step-17", @ref
+ * step_18 "step-18", @ref step_31 "step-31", and @ref step_32 "step-32".
* </ul>
*
* In this case, one possibility is to distribute local entries to the final
*
* <h3>Distributing constraints</h3>
*
- * After solving the condensed system of equations, the solution vector has to
- * be redistributed. This is done by the two distribute() functions, one
+ * After solving the condensed system of equations, the solution vector has
+ * to be redistributed. This is done by the two distribute() functions, one
* working with two vectors, one working in-place. The operation of
* distribution undoes the condensation process in some sense, but it should
- * be noted that it is not the inverse operation. Basically, distribution sets
- * the values of the constrained nodes to the value that is computed from the
- * constraint given the values of the unconstrained nodes. This is usually
- * necessary since the condensed linear systems only describe the equations
- * for unconstrained nodes, and constrained nodes need to get their values in
- * a second step.
+ * be noted that it is not the inverse operation. Basically, distribution
+ * sets the values of the constrained nodes to the value that is computed
+ * from the constraint given the values of the unconstrained nodes plus
+ * possible inhomogeneities. This is usually necessary since the condensed
+ * linear systems only describe the equations for unconstrained nodes, and
+ * constrained nodes need to get their values in a second step.
*
* @ingroup dofs
* @author Wolfgang Bangerth, Martin Kronbichler, 1998, 2004, 2008, 2009
* condenses square compressed
* sparsity patterns.
*
- * Given the data structure used
- * by CompressedSparsityPattern,
- * this function becomes
- * quadratic in the number of
- * degrees of freedom for large
- * problems and can dominate
+ * Given the data structure used by
+ * CompressedSparsityPattern, this
+ * function becomes quadratic in the
+ * number of degrees of freedom for
+ * large problems and can dominate
* setting up linear systems when
- * several hundred thousand or
- * millions of unknowns are
- * involved and for problems with
- * many nonzero elements per row
- * (for example for vector-valued
- * problems or hp finite
+ * several hundred thousand or millions
+ * of unknowns are involved and for
+ * problems with many nonzero elements
+ * per row (for example for
+ * vector-valued problems or hp finite
* elements). In this case, it is
* advisable to use the
- * CompressedSetSparsityPattern
- * class instead, see for example
- * @ref step_27 "step-27" and
- * @ref step_31 "step-31".
+ * CompressedSetSparsityPattern class
+ * instead, see for example @ref
+ * step_27 "step-27", or to use the
+ * CompressedSimpleSparsityPattern
+ * class, see for example @ref step_31
+ * "step-31".
*/
void condense (CompressedSparsityPattern &sparsity) const;
*/
void condense (BlockCompressedSimpleSparsityPattern &sparsity) const;
-
+
/**
- * Condense a given matrix. The associated
- * matrix struct should be condensed and
- * compressed. It is the user's
- * responsibility to guarantee that all
- * entries in the @p condensed matrix be
- * zero!
+ * Condense a given matrix. The
+ * associated matrix struct should be
+ * condensed and compressed. It is the
+ * user's responsibility to guarantee
+ * that all entries in the @p condensed
+ * matrix be zero!
*
* The constraint matrix object must be
* closed to call this function.
void condense (BlockSparseMatrix<number> &matrix) const;
/**
- * Condense the given vector
- * @p uncondensed into @p condensed. It
- * is the user's responsibility to
- * guarantee that all entries of
- * @p condensed be zero.
+ * Condense the given vector @p
+ * uncondensed into @p condensed. It is
+ * the user's responsibility to
+ * guarantee that all entries of @p
+ * condensed be zero. Note that this
+ * function does not take any
+ * inhomogeneity into account, use the
+ * function using both a matrix and
+ * vector for that case.
*
* The @p VectorType may be a
* Vector<float>, Vector<double>,
* PETSc or Trilinos vector
* wrapper class, or any other
* type having the same
- * interface.
+ * interface. Note that this
+ * function does not take any
+ * inhomogeneity into account, use the
+ * function using both a matrix and
+ * vector for that case.
*/
template <class VectorType>
void condense (VectorType &vec) const;
+ /**
+ * Condense a given matrix and a given
+ * vector. The associated matrix struct
+ * should be condensed and
+ * compressed. It is the user's
+ * responsibility to guarantee that all
+ * entries in the @p condensed matrix
+ * and vector be zero! This function is
+ * capable of applying inhomogeneous
+ * constraints.
+ *
+ * The constraint matrix object must be
+ * closed to call this function.
+ */
+ template<typename number, class VectorType>
+ void condense (const SparseMatrix<number> &uncondensed_matrix,
+ const VectorType &uncondensed_vector,
+ SparseMatrix<number> &condensed_matrix,
+ VectorType &condensed_vector) const;
+
+ /**
+ * This function does much the same as
+ * the above one, except that it
+ * condenses matrix and vector
+ * 'in-place'. See the general
+ * documentation of this class for more
+ * detailed information.
+ */
+ template<typename number, class VectorType>
+ void condense (SparseMatrix<number> &matrix,
+ VectorType &vector) const;
+
+ /**
+ * Same function as above, but
+ * condenses square block sparse
+ * matrices and vectors.
+ */
+ template <typename number, class BlockVectorType>
+ void condense (BlockSparseMatrix<number> &matrix,
+ BlockVectorType &vector) const;
+
/**
* @}
*/
*/
/**
- * This function takes a vector of local
- * contributions (@p local_vector)
- * corresponding to the degrees of
- * freedom indices given in @p
- * 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
- * @p local_vector and @p
- * local_dof_indices have the same number
- * of elements, this function is happy
- * with whatever it is given.
+ * This function takes a vector of
+ * local contributions (@p
+ * local_vector) corresponding to the
+ * degrees of freedom indices given in
+ * @p 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 @p
+ * local_vector and @p
+ * local_dof_indices have the same
+ * number of elements, this function is
+ * happy with whatever it is
+ * given. Note that this function will
+ * apply all constraints as if they
+ * were homogeneous. For correctly
+ * setting inhomogeneous constraints,
+ * use the function with both matrix
+ * and vector arguments.
*
- * In contrast to the similar function in
- * the DoFAccessor class, this function
- * also takes care of constraints,
- * i.e. if one of the elements of @p
- * local_dof_indices belongs to a
- * constrained node, then rather than
- * writing the corresponding element of
- * @p local_vector into @p
- * global_vector, the element is
- * distributed to the entries in the
- * global vector to which this particular
- * degree of freedom is constrained.
+ * In contrast to the similar function
+ * in the DoFAccessor class, this
+ * function also takes care of
+ * constraints, i.e. if one of the
+ * elements of @p local_dof_indices
+ * belongs to a constrained node, then
+ * rather than writing the
+ * corresponding element of @p
+ * local_vector into @p 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
* vectors and matrices are fully
* assembled.
*
- * In order to do its work
- * properly, this function has to
- * know which degrees of freedom
- * are fixed, for example
- * boundary values. For this, the
- * third argument is a map
- * between the numbers of the
- * DoFs that are fixed and the
- * values they are fixed to. One
- * can pass an empty map in for
- * this argument, but note that
- * you will then have to fix
- * these nodes later on again,
- * for example by using
+ * In order to do its work properly,
+ * this function has to know which
+ * degrees of freedom are fixed, for
+ * example boundary values. For this,
+ * the third argument is a map between
+ * the numbers of the DoFs that are
+ * fixed and the values they are fixed
+ * to. One can pass an empty map in for
+ * this argument, but note that you
+ * will then have to fix these nodes
+ * later on again, for example by using
* MatrixTools::apply_boundary_values
- * to the resulting
- * matrix. However, since the
- * present function was written
- * for the express purpose of not
- * having to use tools that later
- * modify the matrix, it is
- * advisable to have the list of
- * fixed nodes available when
- * calling the present function.
+ * to the resulting matrix. However,
+ * since the present function was
+ * written for the express purpose of
+ * not having to use tools that later
+ * modify the matrix, it is advisable
+ * to have the list of fixed nodes
+ * available when calling the present
+ * function.
*/
template <typename VectorType>
void
VectorType &global_vector) const;
/**
- * This function takes a matrix of local
- * contributions (@p local_matrix)
- * corresponding to the degrees of
- * freedom indices given in @p
- * 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
- * @p local_matrix and @p
- * local_dof_indices have the same number
- * of elements, this function is happy
- * with whatever it is given.
+ * This function takes a matrix of
+ * local contributions (@p
+ * local_matrix) corresponding to the
+ * degrees of freedom indices given in
+ * @p 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 @p
+ * local_matrix and @p
+ * 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. if one of the elements of @p
- * local_dof_indices belongs to a
- * constrained node, then rather than
- * writing the corresponding element of
- * @p local_matrix into @p
- * global_matrix, the element is
- * distributed to the entries in the
- * global matrix to which this particular
- * degree of freedom is constrained.
+ * In contrast to the similar function
+ * in the DoFAccessor class, this
+ * function also takes care of
+ * constraints, i.e. if one of the
+ * elements of @p local_dof_indices
+ * belongs to a constrained node, then
+ * rather than writing the
+ * corresponding element of @p
+ * local_matrix into @p 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 @p local_dof_indices
- * is constrained, we distribute the
+ * 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 @p 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.
+ * 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
+ * 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>
const std::vector<unsigned int> &local_dof_indices,
MatrixType &global_matrix) const;
+ /**
+ * This function simultaneously writes
+ * elements into matrix and vector,
+ * according to the constraints
+ * specified by the calling
+ * ConstraintMatrix. This function can
+ * correctly handle inhomogeneous
+ * constraints as well.
+ */
+ template <typename MatrixType, typename VectorType>
+ void
+ distribute_local_to_global (const FullMatrix<double> &local_matrix,
+ const Vector<double> &local_vector,
+ const std::vector<unsigned int> &local_dof_indices,
+ MatrixType &global_matrix,
+ VectorType &global_vector) const;
+
/**
* Do a similar operation as the
* distribute_local_to_global() function
* second input argument is not necessary
* here.
*
- * The third argument to this
- * function,
- * keep_constrained_entries
- * determines whether the
- * function shall allocate
- * entries in the sparsity
- * pattern at all for entries
- * that will later be set to zero
- * upon condensation of the
- * matrix. These entries are
- * necessary if the matrix is
- * built unconstrained, and only
- * later condensed. They are not
- * necessary if the matrix is
- * built using the
+ * The third argument to this function,
+ * keep_constrained_entries determines
+ * whether the function shall allocate
+ * entries in the sparsity pattern at
+ * all for entries that will later be
+ * set to zero upon condensation of the
+ * matrix. These entries are necessary
+ * if the matrix is built
+ * unconstrained, and only later
+ * condensed. They are not necessary if
+ * the matrix is built using the
* distribute_local_to_global()
* function of this class which
- * distributes entries right away
- * when copying a local matrix
- * into a global object. The
- * default of this argument is
- * true, meaning to allocate the
- * few entries that may later be
- * set to zero.
+ * distributes entries right away when
+ * copying a local matrix into a global
+ * object. The default of this argument
+ * is true, meaning to allocate the few
+ * entries that may later be set to
+ * zero.
*
* By default, the function adds
- * entries for all pairs of
- * indices given in the first
- * argument to the sparsity
- * pattern (unless
+ * entries for all pairs of indices
+ * given in the first argument to the
+ * sparsity pattern (unless
* keep_constrained_entries is
- * false). However, sometimes one
- * would like to only add a
- * subset of all of these
- * pairs. In that case, the last
- * argument can be used which
- * specifies a boolean mask which
- * of the pairs of indices should
- * be considered. If the mask is
- * false for a pair of indices,
- * then no entry will be added to
- * the sparsity pattern for this
- * pair, irrespective of whether
- * one or both of the indices
- * correspond to constrained
- * degrees of freedom.
+ * false). However, sometimes one would
+ * like to only add a subset of all of
+ * these pairs. In that case, the last
+ * argument can be used which specifies
+ * a boolean mask which of the pairs of
+ * indices should be considered. If the
+ * mask is false for a pair of indices,
+ * then no entry will be added to the
+ * sparsity pattern for this pair,
+ * irrespective of whether one or both
+ * of the indices correspond to
+ * constrained degrees of freedom.
*
* This function is not typically called
* from user code, but is used in the
*/
/**
- * Re-distribute the elements of
- * the vector @p condensed to
- * @p uncondensed. It is the
- * user's responsibility to
- * guarantee that all entries of
- * @p uncondensed be zero!
+ * Re-distribute the elements of the
+ * vector @p condensed to @p
+ * uncondensed. It is the user's
+ * responsibility to guarantee that all
+ * entries of @p uncondensed be zero!
*
- * This function undoes the
- * action of @p condense somehow,
- * but it should be noted that it
- * is not the inverse of
- * @p condense.
+ * This function undoes the action of
+ * @p condense somehow, but it should
+ * be noted that it is not the inverse
+ * of @p condense.
*
* The @p VectorType may be a
* Vector<float>, Vector<double>,
- * BlockVector<tt><...></tt>, a
- * PETSc or Trilinos vector
- * wrapper class, or any other
- * type having the same
+ * BlockVector<tt><...></tt>, a PETSc
+ * or Trilinos vector wrapper class, or
+ * any other type having the same
* interface.
*/
template <class VectorType>
VectorType &uncondensed) const;
/**
- * Re-distribute the elements of
- * the vector in-place. The @p
- * VectorType may be a
- * Vector<float>, Vector<double>,
- * BlockVector<tt><...></tt>, a
- * PETSc or Trilinos vector
- * wrapper class, or any other
- * type having the same
+ * Re-distribute the elements of the
+ * vector in-place. The @p VectorType
+ * may be a Vector<float>,
+ * Vector<double>,
+ * BlockVector<tt><...></tt>, a PETSc
+ * or Trilinos vector wrapper class, or
+ * any other type having the same
* interface.
*/
template <class VectorType>
struct ConstraintLine
{
/**
- * Number of this line. Since only very
- * few lines are stored, we can not
- * assume a specific order and have
- * to store the line number explicitly.
+ * Number of this line. Since only
+ * very few lines are stored, we
+ * can not assume a specific order
+ * and have to store the line
+ * number explicitly.
*/
unsigned int line;
/**
- * Row numbers and values of the entries
- * in this line.
+ * Row numbers and values of the
+ * entries in this line.
*
- * For the reason why we use a vector
- * instead of a map and the consequences
- * thereof, the same applies as what is
- * said for ConstraintMatrix@p ::lines.
+ * For the reason why we use a
+ * vector instead of a map and the
+ * consequences thereof, the same
+ * applies as what is said for
+ * ConstraintMatrix@p ::lines.
*/
std::vector<std::pair<unsigned int,double> > entries;
double inhomogeneity;
/**
- * This operator is a bit
- * weird and unintuitive: it
- * compares the line numbers
- * of two lines. We need this
- * to sort the lines; in fact
- * we could do this using a
- * comparison predicate.
- * However, this way, it is
- * easier, albeit unintuitive
- * since two lines really
- * have no god-given order
+ * This operator is a bit weird and
+ * unintuitive: it compares the
+ * line numbers of two lines. We
+ * need this to sort the lines; in
+ * fact we could do this using a
+ * comparison predicate. However,
+ * this way, it is easier, albeit
+ * unintuitive since two lines
+ * really have no god-given order
* relation.
*/
bool operator < (const ConstraintLine &) const;
/**
- * This operator is likewise
- * weird: it checks whether
- * the line indices of the
- * two operands are equal,
- * irrespective of the fact
- * that the contents of the
- * line may be different.
+ * This operator is likewise weird:
+ * it checks whether the line
+ * indices of the two operands are
+ * equal, irrespective of the fact
+ * that the contents of the line
+ * may be different.
*/
bool operator == (const ConstraintLine &) const;
/**
* Determine an estimate for the
- * memory consumption (in bytes)
- * of this object.
+ * memory consumption (in bytes) of
+ * this object.
*/
unsigned int memory_consumption () const;
};
/**
* Store the lines of the matrix.
- * Entries are usually
- * appended in an arbitrary order and
- * insertion into a vector is done best
- * at the end, so the order is
- * unspecified after all entries are
- * inserted. Sorting of the entries takes
- * place when calling the <tt>close()</tt> function.
+ * Entries are usually appended in an
+ * arbitrary order and insertion into a
+ * vector is done best at the end, so
+ * the order is unspecified after all
+ * entries are inserted. Sorting of the
+ * entries takes place when calling the
+ * <tt>close()</tt> function.
*
- * We could, instead of using a vector, use
- * an associative array, like a map to
- * store the lines. This, however, would
- * mean a much more fractioned heap since it
- * allocates many small objects, ans would
- * additionally make usage of this matrix
- * much slower.
+ * We could, instead of using a vector,
+ * use an associative array, like a map
+ * to store the lines. This, however,
+ * would mean a much more fractioned
+ * heap since it allocates many small
+ * objects, ans would additionally make
+ * usage of this matrix much slower.
*/
std::vector<ConstraintLine> lines;
/**
* A list of flags that indicate
- * whether there is a constraint
- * line for a given degree of
- * freedom index. Note that this
- * class has no notion of how
- * many degrees of freedom there
- * really are, so if we check
- * whether there is a constraint
- * line for a given degree of
- * freedom, then this vector may
- * actually be shorter than the
- * index of the DoF we check for.
+ * whether there is a constraint line
+ * for a given degree of freedom
+ * index. Note that this class has no
+ * notion of how many degrees of
+ * freedom there really are, so if we
+ * check whether there is a constraint
+ * line for a given degree of freedom,
+ * then this vector may actually be
+ * shorter than the index of the DoF we
+ * check for.
*
- * This field exists since when adding a
- * new constraint line we have to figure
- * out whether it already
+ * This field exists since when adding
+ * a new constraint line we have to
+ * figure out whether it already
* exists. Previously, we would simply
* walk the unsorted list of constraint
* lines until we either hit the end or
- * found it. This algorithm is O(N) if N
- * is the number of constraints, which
- * makes it O(N^2) when inserting all
- * constraints. For large problems with
- * many constraints, this could easily
- * take 5-10 per cent of the total run
- * time. With this field, we can at least
- * save this time when checking whether a
- * new constraint line already exists.
+ * found it. This algorithm is O(N) if
+ * N is the number of constraints,
+ * which makes it O(N^2) when inserting
+ * all constraints. For large problems
+ * with many constraints, this could
+ * easily take 5-10 per cent of the
+ * total run time. With this field, we
+ * can at least save this time when
+ * checking whether a new constraint
+ * line already exists.
*
* To make things worse, traversing the
- * list of existing constraints requires
- * reads from many different places in
- * memory. Thus, in large 3d
- * applications, the add_line() function
- * showed up very prominently in the
- * overall compute time, mainly because
- * it generated a lot of cache
+ * list of existing constraints
+ * requires reads from many different
+ * places in memory. Thus, in large 3d
+ * applications, the add_line()
+ * function showed up very prominently
+ * in the overall compute time, mainly
+ * because it generated a lot of cache
* misses. This should also be fixed by
- * using the O(1) algorithm to access the
- * fields of this array.
+ * using the O(1) algorithm to access
+ * the fields of this array.
*
* The field is useful in a number of
* other contexts as well, though.
bool sorted;
/**
- * Return @p true if the weight
- * of an entry (the second
- * element of the pair) equals
- * zero. This function is used to
- * delete entries with zero
+ * Return @p true if the weight of an
+ * entry (the second element of the
+ * pair) equals zero. This function is
+ * used to delete entries with zero
* weight.
*/
static bool check_zero_weight (const std::pair<unsigned int, double> &p);
/**
- * Dummy table that serves as
- * default argument for function
+ * Dummy table that serves as default
+ * argument for function
* <tt>add_entries_local_to_global()</tt>.
*/
static const Table<2,bool> default_empty_table;
const std::vector<ConstraintLine>::const_iterator start=lines.begin();
// the usual case is that the line where
- // a value is entered is the one we
+ // the inhomogeneity is set to the one we
// added last, so we search backward
for (line_ptr=(lines.end()-1); line_ptr!=start; --line_ptr)
if (line_ptr->line == line)
next_constraint->entries[q].second *
c->entries[p].second);
};
-
+
++next_constraint;
};
}
// irregular row @p row and
// irregular column
// @p column set old entry
- // to one if on main
+ // to one on main
// diagonal, zero otherwise
{
for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
+= (vec(constraint_line->line) *
constraint_line->entries[q].second);
vec(constraint_line->line) = 0.;
+
+ // in case the constraint is
+ // inhomogeneous, this function is not
+ // appropriate. Throw an exception.
+ Assert (constraint_line->inhomogeneity == 0.,
+ ExcMessage ("Inhomogeneous constraint cannot be condensed "
+ "without any matrix specified."));
}
}
-template <class VectorType>
+template<typename number, class VectorType>
void
-ConstraintMatrix::set_zero (VectorType &vec) const
+ConstraintMatrix::condense (const SparseMatrix<number> &uncondensed,
+ const VectorType &uncondensed_vector,
+ SparseMatrix<number> &condensed,
+ VectorType &condensed_vector) const
{
+ const SparsityPattern &uncondensed_struct = uncondensed.get_sparsity_pattern ();
+
Assert (sorted == true, ExcMatrixNotClosed());
+ Assert (uncondensed_struct.is_compressed() == true, ExcMatrixNotClosed());
+ Assert (condensed.get_sparsity_pattern().is_compressed() == true, ExcMatrixNotClosed());
+ Assert (uncondensed_struct.n_rows() == uncondensed_struct.n_cols(),
+ ExcNotQuadratic());
+ Assert (condensed.n() == condensed.m(),
+ ExcNotQuadratic());
+ Assert (condensed.n()+n_constraints() == uncondensed.n(),
+ ExcDimensionMismatch(condensed.n()+n_constraints(), uncondensed.n()));
+ Assert (condensed_vector.size()+n_constraints() == uncondensed_vector.size(),
+ ExcDimensionMismatch(condensed_vector.size()+n_constraints(),
+ uncondensed_vector.size()));
+ Assert (condensed_vector.size() == condensed.m(),
+ ExcDimensionMismatch(condensed_vector.size(), condensed.m()));
- std::vector<ConstraintLine>::const_iterator constraint_line = lines.begin();
- for (; constraint_line!=lines.end(); ++constraint_line)
- vec(constraint_line->line) = 0.;
-}
-
+ // store for each line of the matrix
+ // its new line number
+ // after compression. If the shift is
+ // -1, this line will be condensed away
+ std::vector<int> new_line;
+ new_line.reserve (uncondensed_struct.n_rows());
-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(),
- ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
- Assert (sorted == true, ExcMatrixNotClosed());
+ std::vector<ConstraintLine>::const_iterator next_constraint = lines.begin();
+ unsigned int shift = 0;
+ const unsigned int n_rows = uncondensed_struct.n_rows();
- 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);
- }
+ if (next_constraint == lines.end())
+ // if no constraint is to be handled
+ for (unsigned int row=0; row!=n_rows; ++row)
+ new_line.push_back (row);
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];
+ for (unsigned int row=0; row!=n_rows; ++row)
+ if (row == next_constraint->line)
+ {
+ // this line is constrained
+ new_line.push_back (-1);
+ // note that @p lines is ordered
+ ++shift;
+ ++next_constraint;
+ if (next_constraint == lines.end())
+ // nothing more to do; finish rest
+ // of loop
+ {
+ for (unsigned int i=row+1; i<n_rows; ++i)
+ new_line.push_back (i-shift);
+ break;
+ };
+ }
+ else
+ new_line.push_back (row-shift);
- 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, but make
- // sure we don't touch the
- // entries of fixed dofs
- //
- // there is one critical
- // point: sometimes a dof
- // may be both constrained
- // and fixed, for example
- // hanging nodes in 3d at
- // the boundary. in that
- // case, we don't quite
- // know what to do --
- // handle the constraint or
- // the fixed
- // value. however, this
- // isn't so hard if all the
- // nodes that this node is
- // constrained to are also
- // fixed nodes, in which
- // case we could do both
- // but opt to copy the
- // element. however, we
- // have to check that all
- // the nodes to which it is
- // constrained are also
- // fixed
- if ((position == lines.end())
- ||
- (position->line != local_dof_indices[i]))
- global_vector(local_dof_indices[i]) += local_vector(i);
- else
+ next_constraint = lines.begin();
+
+ // note: in this loop we need not check
+ // whether @p next_constraint is a valid
+ // iterator, since @p next_constraint is
+ // only evaluated so often as there are
+ // entries in new_line[*] which tells us
+ // which constraints exist
+ for (unsigned int row=0; row<uncondensed_struct.n_rows(); ++row)
+ if (new_line[row] != -1)
+ {
+ // line not constrained
+ // copy entries if column will not
+ // be condensed away, distribute
+ // otherwise
+ for (unsigned int j=uncondensed_struct.get_rowstart_indices()[row];
+ j<uncondensed_struct.get_rowstart_indices()[row+1]; ++j)
+ if (new_line[uncondensed_struct.get_column_numbers()[j]] != -1)
+ condensed.add (new_line[row], new_line[uncondensed_struct.get_column_numbers()[j]],
+ uncondensed.global_entry(j));
+ else
{
- for (unsigned int j=0; j<position->entries.size(); ++j)
- global_vector(position->entries[j].first)
- += local_vector(i) * position->entries[j].second;
+ // let c point to the
+ // constraint of this column
+ std::vector<ConstraintLine>::const_iterator c = lines.begin();
+ while (c->line != uncondensed_struct.get_column_numbers()[j])
+ ++c;
+
+ for (unsigned int q=0; q!=c->entries.size(); ++q)
+ // distribute to rows with
+ // appropriate weight
+ condensed.add (new_line[row], new_line[c->entries[q].first],
+ uncondensed.global_entry(j) * c->entries[q].second);
+
+ // take care of inhomogeneity:
+ // need to subtract this element from the
+ // vector. this corresponds to an
+ // explicit elimination in the respective
+ // row of the inhomogeneous constraint in
+ // the matrix with Gauss elimination
+ condensed_vector(new_line[row]) -= uncondensed.global_entry(j) /
+ uncondensed.diag_element(row) * c->inhomogeneity;
}
- }
- }
+
+ condensed_vector(new_line[row]) += uncondensed_vector(row);
+ }
+ else
+ // line must be distributed
+ {
+ for (unsigned int j=uncondensed_struct.get_rowstart_indices()[row];
+ j<uncondensed_struct.get_rowstart_indices()[row+1]; ++j)
+ // for each column: distribute
+ if (new_line[uncondensed_struct.get_column_numbers()[j]] != -1)
+ // column is not constrained
+ for (unsigned int q=0; q!=next_constraint->entries.size(); ++q)
+ condensed.add (new_line[next_constraint->entries[q].first],
+ new_line[uncondensed_struct.get_column_numbers()[j]],
+ uncondensed.global_entry(j) *
+ next_constraint->entries[q].second);
+
+ else
+ // not only this line but
+ // also this col is constrained
+ {
+ // let c point to the constraint
+ // of this column
+ std::vector<ConstraintLine>::const_iterator c = lines.begin();
+ while (c->line != uncondensed_struct.get_column_numbers()[j])
+ ++c;
+
+ for (unsigned int p=0; p!=c->entries.size(); ++p)
+ for (unsigned int q=0; q!=next_constraint->entries.size(); ++q)
+ condensed.add (new_line[next_constraint->entries[q].first],
+ new_line[c->entries[p].first],
+ uncondensed.global_entry(j) *
+ next_constraint->entries[q].second *
+ c->entries[p].second);
+ };
+
+ // distribute vector
+ for (unsigned int q=0; q!=next_constraint->entries.size(); ++q)
+ condensed_vector(new_line[next_constraint->entries[q].first])
+ +=
+ uncondensed_vector(row) * next_constraint->entries[q].second;
+
+ ++next_constraint;
+ };
}
-template <typename MatrixType>
+template<typename number, class VectorType>
void
-ConstraintMatrix::
-distribute_local_to_global (const FullMatrix<double> &local_matrix,
- const std::vector<unsigned int> &local_dof_indices,
- MatrixType &global_matrix) const
+ConstraintMatrix::condense (SparseMatrix<number> &uncondensed,
+ VectorType &vec) const
{
- Assert (local_matrix.n() == local_dof_indices.size(),
- ExcDimensionMismatch(local_matrix.n(), local_dof_indices.size()));
- Assert (local_matrix.m() == local_dof_indices.size(),
- ExcDimensionMismatch(local_matrix.m(), local_dof_indices.size()));
- Assert (sorted == true, ExcMatrixNotClosed());
+ const SparsityPattern &sparsity = uncondensed.get_sparsity_pattern ();
- const unsigned int n_local_dofs = local_dof_indices.size();
+ Assert (sorted == true, ExcMatrixNotClosed());
+ Assert (sparsity.is_compressed() == true, ExcMatrixNotClosed());
+ Assert (sparsity.n_rows() == sparsity.n_cols(),
+ ExcNotQuadratic());
+ Assert (vec.size() == sparsity.n_rows(),
+ ExcDimensionMismatch(vec.size(), sparsity.n_rows()));
- // A lock that allows only one thread at
- // time to go on in this function.
- mutex.acquire();
+ double average_diagonal = 0;
+ for (unsigned int i=0; i<uncondensed.m(); ++i)
+ average_diagonal += std::fabs (uncondensed.diag_element(i));
+ average_diagonal /= uncondensed.m();
- // have a special case where there are no
- // constraints at all, since then we can be
- // a lot faster
- if (lines.size() == 0)
- global_matrix.add(local_dof_indices, local_matrix);
- else
+ // store for each index whether it must be
+ // distributed or not. If entry is
+ // invalid_unsigned_int, no distribution is
+ // necessary. otherwise, the number states
+ // which line in the constraint matrix
+ // handles this index
+ std::vector<unsigned int> distribute (sparsity.n_rows(),
+ numbers::invalid_unsigned_int);
+
+ for (unsigned int c=0; c<lines.size(); ++c)
+ distribute[lines[c].line] = c;
+
+ const unsigned int n_rows = sparsity.n_rows();
+ for (unsigned int row=0; row<n_rows; ++row)
{
- // 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. Moreover, we
- // store how many entries there are
- // at most in one constrained row in
- // order to set the scratch array for
- // column data to a sufficient size.
- std::vector<const ConstraintLine *>
- constraint_lines (n_local_dofs,
- static_cast<const ConstraintLine *>(0));
- unsigned int n_max_entries_per_row = 0;
- for (unsigned int i=0; i<n_local_dofs; ++i)
+ if (distribute[row] == numbers::invalid_unsigned_int)
+ // regular line. loop over cols
{
- ConstraintLine index_comparison;
- index_comparison.line = local_dof_indices[i];
+ for (typename SparseMatrix<number>::iterator
+ entry = uncondensed.begin(row);
+ entry != uncondensed.end(row); ++entry)
+ {
+ const unsigned int column = entry->column();
+
+ // end of row reached?
+ // this should not
+ // happen, since we only
+ // operate on compressed
+ // matrices!
+ Assert (column != SparsityPattern::invalid_entry,
+ ExcMatrixNotClosed());
+
+ if (distribute[column] != numbers::invalid_unsigned_int)
+ // distribute entry at
+ // regular row @p row
+ // and irregular column
+ // sparsity.get_column_numbers()[j];
+ // set old entry to
+ // zero
+ {
+ for (unsigned int q=0;
+ q!=lines[distribute[column]].entries.size(); ++q)
+ uncondensed.add (row,
+ lines[distribute[column]].entries[q].first,
+ entry->value() *
+ lines[distribute[column]].entries[q].second);
+
+ // need to subtract this element from the
+ // vector. this corresponds to an
+ // explicit elimination in the respective
+ // row of the inhomogeneous constraint in
+ // the matrix with Gauss elimination
+ vec(column) -= entry->value() *
+ lines[distribute[column]].inhomogeneity;
+
+ // set old value to zero
+ entry->value() = 0.;
+ }
+ }
+ }
+ else
+ // row must be distributed
+ {
+ for (typename SparseMatrix<number>::iterator
+ entry = uncondensed.begin(row);
+ entry != uncondensed.end(row); ++entry)
+ {
+ const unsigned int column = entry->column();
+
+ // end of row reached?
+ // this should not
+ // happen, since we only
+ // operate on compressed
+ // matrices!
+ Assert (column != SparsityPattern::invalid_entry,
+ ExcMatrixNotClosed());
+
+ if (distribute[column] == numbers::invalid_unsigned_int)
+ // distribute entry at
+ // irregular row
+ // @p row and regular
+ // column
+ // column. 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,
+ column,
+ entry->value() *
+ lines[distribute[row]].entries[q].second);
+
+ // set old entry to zero
+ entry->value() = 0.;
+ }
+ else
+ // distribute entry at
+ // irregular row @p row and
+ // irregular column
+ // @p column set old entry
+ // to one 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[column]].entries.size(); ++q)
+ uncondensed.add (lines[distribute[row]].entries[p].first,
+ lines[distribute[column]].entries[q].first,
+ entry->value() *
+ lines[distribute[row]].entries[p].second *
+ lines[distribute[column]].entries[q].second);
+
+ // set old entry to correct
+ // value
+ entry->value() = (row == column ? average_diagonal : 0. );
+ }
+ }
+
+ // take care of vector
+ for (unsigned int q=0; q!=lines[distribute[row]].entries.size(); ++q)
+ vec(lines[distribute[row]].entries[q].first)
+ += (vec(row) * lines[distribute[row]].entries[q].second);
+
+ vec(lines[distribute[row]].line) = 0.;
+ }
+ }
+}
+
+
+
+template <typename number, class BlockVectorType>
+void
+ConstraintMatrix::condense (BlockSparseMatrix<number> &uncondensed,
+ BlockVectorType &vec) const
+{
+ const unsigned int blocks = uncondensed.n_block_rows();
+
+ const BlockSparsityPattern &
+ sparsity = uncondensed.get_sparsity_pattern ();
+
+ Assert (sorted == true, ExcMatrixNotClosed());
+ Assert (sparsity.is_compressed() == true, ExcMatrixNotClosed());
+ Assert (sparsity.n_rows() == sparsity.n_cols(),
+ ExcNotQuadratic());
+ Assert (sparsity.n_block_rows() == sparsity.n_block_cols(),
+ ExcNotQuadratic());
+ Assert (sparsity.n_block_rows() == sparsity.n_block_cols(),
+ ExcNotQuadratic());
+ Assert (sparsity.get_column_indices() == sparsity.get_row_indices(),
+ ExcNotQuadratic());
+ Assert (vec.size() == sparsity.n_rows(),
+ ExcDimensionMismatch(vec.size(), sparsity.n_rows()));
+ Assert (vec.n_blocks() == sparsity.n_block_rows(),
+ ExcDimensionMismatch(vec.n_blocks(), sparsity.n_block_rows()));
+
+ double average_diagonal = 0;
+ for (unsigned int b=0; b<uncondensed.n_block_rows(); ++b)
+ for (unsigned int i=0; i<uncondensed.block(b,b).m(); ++i)
+ average_diagonal += std::fabs (uncondensed.block(b,b).diag_element(i));
+ average_diagonal /= uncondensed.m();
+
+ const BlockIndices &
+ index_mapping = sparsity.get_column_indices();
+
+ // store for each index whether it must be
+ // distributed or not. If entry is
+ // numbers::invalid_unsigned_int,
+ // no distribution is necessary.
+ // otherwise, the number states which line
+ // in the constraint matrix handles this
+ // index
+ std::vector<unsigned int> distribute (sparsity.n_rows(),
+ numbers::invalid_unsigned_int);
+
+ for (unsigned int c=0; c<lines.size(); ++c)
+ distribute[lines[c].line] = 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 std::pair<unsigned int,unsigned int>
+ block_index = index_mapping.global_to_local(row);
+ const unsigned int block_row = block_index.first;
+
+ if (distribute[row] == numbers::invalid_unsigned_int)
+ // 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)
+ {
+ for (typename SparseMatrix<number>::iterator
+ entry = uncondensed.block(block_row, block_col).begin(block_index.second);
+ entry != uncondensed.block(block_row, block_col).end(block_index.second);
+ ++entry)
+ {
+ const unsigned int global_col
+ = index_mapping.local_to_global(block_col,entry->column());
+
+ if (distribute[global_col] != numbers::invalid_unsigned_int)
+ // distribute entry at
+ // regular row @p row
+ // and irregular column
+ // global_col; set old
+ // entry to zero
+ {
+ const double old_value = entry->value ();
+
+ 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);
+
+ // need to subtract this element from the
+ // vector. this corresponds to an
+ // explicit elimination in the respective
+ // row of the inhomogeneous constraint in
+ // the matrix with Gauss elimination
+ vec(global_col) -= entry->value() *
+ lines[distribute[global_col]].inhomogeneity;
+
+ entry->value() = 0.;
+ }
+ }
+ }
+ }
+ 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)
+ {
+ for (typename SparseMatrix<number>::iterator
+ entry = uncondensed.block(block_row, block_col).begin(block_index.second);
+ entry != uncondensed.block(block_row, block_col).end(block_index.second);
+ ++entry)
+ {
+ const unsigned int global_col
+ = index_mapping.local_to_global (block_col, entry->column());
+
+ if (distribute[global_col] ==
+ numbers::invalid_unsigned_int)
+ // distribute
+ // entry at
+ // irregular
+ // row @p row
+ // and regular
+ // column
+ // global_col. set
+ // old entry to
+ // zero
+ {
+ const double old_value = entry->value();
+
+ for (unsigned int q=0;
+ q!=lines[distribute[row]].entries.size(); ++q)
+ uncondensed.add (lines[distribute[row]].entries[q].first,
+ global_col,
+ old_value *
+ lines[distribute[row]].entries[q].second);
+
+ entry->value() = 0.;
+ }
+ else
+ // distribute entry at
+ // irregular row @p row
+ // and irregular column
+ // @p global_col set old
+ // entry to one if on
+ // main diagonal, zero
+ // otherwise
+ {
+ const double old_value = entry->value ();
+
+ 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)
+ uncondensed.add (lines[distribute[row]].entries[p].first,
+ lines[distribute[global_col]].entries[q].first,
+ old_value *
+ lines[distribute[row]].entries[p].second *
+ lines[distribute[global_col]].entries[q].second);
+
+ entry->value() = (row == global_col ? average_diagonal : 0. );
+ }
+ }
+ }
+
+ // take care of vector
+ for (unsigned int q=0; q!=lines[distribute[row]].entries.size(); ++q)
+ vec(lines[distribute[row]].entries[q].first)
+ += (vec(row) * lines[distribute[row]].entries[q].second);
+
+ vec(lines[distribute[row]].line) = 0.;
+ }
+ }
+}
+
+
+
+template <class VectorType>
+void
+ConstraintMatrix::set_zero (VectorType &vec) const
+{
+ Assert (sorted == true, ExcMatrixNotClosed());
+
+ std::vector<ConstraintLine>::const_iterator constraint_line = lines.begin();
+ for (; constraint_line!=lines.end(); ++constraint_line)
+ vec(constraint_line->line) = 0.;
+}
+
+
+
+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(),
+ ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
+ 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, but make
+ // sure we don't touch the
+ // entries of fixed dofs
+ //
+ // there is one critical
+ // point: sometimes a dof
+ // may be both constrained
+ // and fixed, for example
+ // hanging nodes in 3d at
+ // the boundary. in that
+ // case, we don't quite
+ // know what to do --
+ // handle the constraint or
+ // the fixed
+ // value. however, this
+ // isn't so hard if all the
+ // nodes that this node is
+ // constrained to are also
+ // fixed nodes, in which
+ // case we could do both
+ // but opt to copy the
+ // element. however, we
+ // have to check that all
+ // the nodes to which it is
+ // constrained are also
+ // fixed
+ if ((position == lines.end())
+ ||
+ (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(),
+ ExcDimensionMismatch(local_matrix.n(), local_dof_indices.size()));
+ Assert (local_matrix.m() == local_dof_indices.size(),
+ ExcDimensionMismatch(local_matrix.m(), local_dof_indices.size()));
+ Assert (global_matrix.m() == global_matrix.n(), ExcNotQuadratic());
+ Assert (sorted == true, ExcMatrixNotClosed());
+
+ const unsigned int n_local_dofs = local_dof_indices.size();
+
+ // A lock that allows only one thread at
+ // time to go on in this function.
+ Threads::ThreadMutex::ScopedLock lock(mutex);
+
+ // have a special case where there are no
+ // constraints at all, since then we can be
+ // a lot faster
+ if (lines.size() == 0)
+ global_matrix.add(local_dof_indices, local_matrix);
+ 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. Moreover, we
+ // store how many entries there are
+ // at most in one constrained row in
+ // order to set the scratch array for
+ // column data to a sufficient size.
+ std::vector<const ConstraintLine *>
+ constraint_lines (n_local_dofs,
+ static_cast<const ConstraintLine *>(0));
+ unsigned int n_max_entries_per_row = 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(),
false);
}
}
- mutex.release();
+}
+
+
+
+template <typename MatrixType, class VectorType>
+void
+ConstraintMatrix::
+distribute_local_to_global (const FullMatrix<double> &local_matrix,
+ const Vector<double> &local_vector,
+ const std::vector<unsigned int> &local_dof_indices,
+ MatrixType &global_matrix,
+ VectorType &global_vector) const
+{
+ Assert (local_matrix.n() == local_dof_indices.size(),
+ ExcDimensionMismatch(local_matrix.n(), local_dof_indices.size()));
+ Assert (local_matrix.m() == local_dof_indices.size(),
+ ExcDimensionMismatch(local_matrix.m(), local_dof_indices.size()));
+ Assert (global_matrix.m() == global_matrix.n(), ExcNotQuadratic());
+ Assert (local_matrix.m() == local_vector.size(),
+ ExcDimensionMismatch(local_matrix.m(), local_vector.size()));
+ Assert (global_matrix.m() == global_vector.size(),
+ ExcDimensionMismatch(global_matrix.m(), global_vector.size()));
+ Assert (sorted == true, ExcMatrixNotClosed());
+
+ const unsigned int n_local_dofs = local_dof_indices.size();
+
+ // A lock that allows only one thread at
+ // time to go on in this function.
+ Threads::ThreadMutex::ScopedLock lock(mutex);
+
+ // have a special case where there are no
+ // constraints at all, since then we can be
+ // a lot faster
+ if (lines.size() == 0)
+ {
+ global_matrix.add(local_dof_indices, local_matrix);
+ for (unsigned int i=0; i<local_dof_indices.size(); ++i)
+ global_vector(local_dof_indices[i]) += local_vector(i);
+ }
+ 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. Moreover, we
+ // store how many entries there are
+ // at most in one constrained row in
+ // order to set the scratch array for
+ // column data to a sufficient size.
+ std::vector<const ConstraintLine *>
+ constraint_lines (n_local_dofs,
+ static_cast<const ConstraintLine *>(0));
+ unsigned int n_max_entries_per_row = 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;
+ n_max_entries_per_row += position->entries.size();
+ }
+ }
+
+ // We need to add the number of
+ // entries in the local matrix in
+ // order to obtain a sufficient size
+ // for the scratch array.
+ n_max_entries_per_row += n_local_dofs;
+ if (column_indices.size() < n_max_entries_per_row)
+ {
+ column_indices.resize(n_max_entries_per_row);
+ column_values.resize(n_max_entries_per_row);
+ }
+
+ // now distribute entries row by row
+ 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);
+
+ unsigned int col_counter = 0;
+
+ for (unsigned int j=0; j<n_local_dofs; ++j)
+ {
+ // we don't need to proceed when the
+ // matrix element is zero
+ if (local_matrix(i,j) == 0)
+ continue;
+
+ 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, so
+ // write the value into
+ // the scratch array
+ column_indices[col_counter] = local_dof_indices[j];
+ column_values[col_counter] = local_matrix(i,j);
+ col_counter++;
+ }
+ else if ((is_constrained_i == true) &&
+ (is_constrained_j == false))
+ {
+ // ok, row is
+ // constrained, but
+ // column is not. This
+ // creates entries in
+ // several rows to the
+ // same column, which is
+ // not covered by the
+ // scratch array. Write
+ // the values directly
+ // into the matrix
+ 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. This
+ // time, we can put
+ // everything into the
+ // scratch array, since
+ // we are in the correct
+ // row.
+ for (unsigned int q=0; q<position_j->entries.size(); ++q)
+ {
+ column_indices[col_counter] = position_j->entries[q].first;
+ column_values[col_counter] = local_matrix(i,j) *
+ position_j->entries[q].second;
+ col_counter++;
+ }
+
+ // need to subtract this element from the
+ // vector. this corresponds to an
+ // explicit elimination in the respective
+ // row of the inhomogeneous constraint in
+ // the matrix with Gauss elimination
+ global_vector(local_dof_indices[i]) -= local_matrix(j,i) *
+ position_j->inhomogeneity;
+ }
+ else if ((is_constrained_i == true) &&
+ (is_constrained_j == true))
+ {
+ // last case: both row
+ // and column are
+ // constrained. Again,
+ // this creates entries
+ // in other rows than the
+ // current one, so write
+ // the values again in
+ // the matrix directly
+ 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
+ //
+ // note that this also
+ // captures the special
+ // case that a dof is
+ // both constrained and
+ // fixed (this can happen
+ // for hanging nodes in
+ // 3d that also happen to
+ // be on the
+ // boundary). in that
+ // case, following the
+ // above program flow, it
+ // is realized that when
+ // distributing the row
+ // and column no elements
+ // of the matrix are
+ // actually touched if
+ // all the degrees of
+ // freedom to which this
+ // dof is constrained are
+ // also constrained (the
+ // usual case with
+ // hanging nodes in
+ // 3d). however, in the
+ // line below, we do
+ // actually do something
+ // with this dof
+ if (i == j)
+ {
+ column_indices[col_counter] = local_dof_indices[j];
+ if (std::fabs (local_matrix(i,j)) < 1e-8)
+ column_values[col_counter] = 1;
+ else
+ column_values[col_counter] = local_matrix(i,j);
+ col_counter++;
+ }
+ }
+ else
+ Assert (false, ExcInternalError());
+ }
+
+ // Check whether we did remain within the
+ // arrays when adding elements into the
+ // scratch arrays. Moreover, there should
+ // be at least one element in the scratch
+ // array (the element diagonal).
+ Assert (col_counter <= n_max_entries_per_row, ExcInternalError());
+
+ // Finally, write the scratch array into
+ // the sparse matrix.
+ if (col_counter > 0)
+ global_matrix.add(local_dof_indices[i], col_counter,
+ &column_indices[0], &column_values[0],
+ false);
+
+ // And we take care of the vector
+ if (is_constrained_i == true)
+ for (unsigned int q=0; q<position_i->entries.size(); ++q)
+ global_vector(position_i->entries[q].first)
+ += local_vector(i) * position_i->entries[q].second;
+ else
+ global_vector(local_dof_indices[i]) += local_vector(i);
+ }
+ }
}
projects / dealii-svn.git / commitdiff