]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Complete the implementation of inhomogeneous constraints in the ConstraintMatrix...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 21 Feb 2009 16:12:26 +0000 (16:12 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 21 Feb 2009 16:12:26 +0000 (16:12 +0000)
git-svn-id: https://svn.dealii.org/trunk@18403 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_constraints.h
deal.II/deal.II/include/dofs/dof_constraints.templates.h
deal.II/deal.II/source/dofs/dof_constraints.cc

index ea3ec36ac8b8dd742ca73a24126f94bc6adf80cd..bb0c96144fc22b6e3c9b887308275546472d0422 100644 (file)
@@ -2,7 +2,7 @@
 //    $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
@@ -43,11 +43,17 @@ class BlockIndices;
 
 
 /**
- * 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&lt;dim,spacedim&gt;(1) and FE_Q&lt;dim,spacedim&gt;(2)) on the two
  * marked cells of the mesh
  *
@@ -69,20 +75,32 @@ class BlockIndices;
  * 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
@@ -186,15 +204,15 @@ class BlockIndices;
  * 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>
@@ -205,11 +223,24 @@ class BlockIndices;
  * 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>
@@ -226,13 +257,12 @@ class BlockIndices;
  * 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
@@ -258,16 +288,16 @@ class BlockIndices;
  * 
  * <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
@@ -701,25 +731,25 @@ class ConstraintMatrix : public Subscriptor
                                      * 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;
 
@@ -780,14 +810,14 @@ class ConstraintMatrix : public Subscriptor
                                      */
     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.
@@ -815,11 +845,15 @@ class ConstraintMatrix : public Subscriptor
     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>,
@@ -842,11 +876,56 @@ class ConstraintMatrix : public Subscriptor
                                      * 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;
+
                                     /**
                                      * @}
                                      */
@@ -857,32 +936,40 @@ class ConstraintMatrix : public Subscriptor
                                      */
     
                                      /**
-                                      * 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
@@ -891,30 +978,26 @@ class ConstraintMatrix : public Subscriptor
                                       * 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
@@ -923,61 +1006,65 @@ class ConstraintMatrix : public Subscriptor
                                 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>
@@ -986,6 +1073,23 @@ class ConstraintMatrix : public Subscriptor
                                 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
@@ -1020,52 +1124,43 @@ class ConstraintMatrix : public Subscriptor
                                      * 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
@@ -1106,25 +1201,22 @@ class ConstraintMatrix : public Subscriptor
                                      */
     
                                     /**
-                                     * 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>
@@ -1132,14 +1224,13 @@ class ConstraintMatrix : public Subscriptor
                     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>
@@ -1221,21 +1312,23 @@ class ConstraintMatrix : public Subscriptor
     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;
 
@@ -1245,101 +1338,98 @@ class ConstraintMatrix : public Subscriptor
         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.
@@ -1353,18 +1443,17 @@ class ConstraintMatrix : public Subscriptor
     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;
@@ -1506,7 +1595,7 @@ ConstraintMatrix::set_inhomogeneity (const unsigned int line,
   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)
index 69c40bb38c2c050f85d2cd985e404df35a5bc7da..60c4b5716d8d32e6899b28b97ba210985a5b929c 100644 (file)
@@ -145,7 +145,7 @@ ConstraintMatrix::condense (const SparseMatrix<number> &uncondensed,
                                 next_constraint->entries[q].second *
                                 c->entries[p].second);
            };
-       
+
        ++next_constraint;
       };
 }
@@ -260,7 +260,7 @@ ConstraintMatrix::condense (SparseMatrix<number> &uncondensed) const
                                                  // 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)
@@ -537,155 +537,649 @@ ConstraintMatrix::condense (VectorType &vec) const
           += (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(),
@@ -868,7 +1362,273 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                              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);
+        }
+    }
 }
 
 
index 846b9df307d230e4d02b4ca3c5f864f337b53de8..46e2b501373a6536c711a45301346c94b9c9acdc 100644 (file)
@@ -1997,6 +1997,14 @@ ConstraintMatrix::memory_consumption () const
   template void ConstraintMatrix::condense<VectorType >(const VectorType &uncondensed,\
                                                       VectorType       &condensed) const;\
   template void ConstraintMatrix::condense<VectorType >(VectorType &vec) const;\
+  template void ConstraintMatrix::condense<float,VectorType >(const SparseMatrix<float> &uncondensed, \
+                                                             const VectorType &uncondensed_vector, \
+                                                             SparseMatrix<float> &condensed, \
+                                                             VectorType       &condensed_vector) const; \
+  template void ConstraintMatrix::condense<double,VectorType >(const SparseMatrix<double> &uncondensed, \
+                                                              const VectorType &uncondensed_vector, \
+                                                              SparseMatrix<double> &condensed, \
+                                                              VectorType       &condensed_vector) const; \
   template void ConstraintMatrix::set_zero<VectorType >(VectorType &vec) const;\
   template void ConstraintMatrix:: \
     distribute_local_to_global<VectorType > (const Vector<double>            &, \
@@ -2026,6 +2034,20 @@ VECTOR_FUNCTIONS(TrilinosWrappers::MPI::Vector);
 VECTOR_FUNCTIONS(TrilinosWrappers::MPI::BlockVector);
 #endif
 
+#define CONDENSE_FUNCTIONS(VectorType, number, MatrixType)             \
+  template void ConstraintMatrix::condense<number,VectorType >(MatrixType &uncondensed, \
+                                                              VectorType &vec) const \
+
+
+CONDENSE_FUNCTIONS(Vector<float>,float,SparseMatrix<float>);
+CONDENSE_FUNCTIONS(Vector<double>,float,SparseMatrix<float>);
+CONDENSE_FUNCTIONS(Vector<double>,double,SparseMatrix<double>);
+CONDENSE_FUNCTIONS(Vector<float>,double,SparseMatrix<double>);
+CONDENSE_FUNCTIONS(BlockVector<double>,float,BlockSparseMatrix<float>);
+CONDENSE_FUNCTIONS(BlockVector<float>,float,BlockSparseMatrix<float>);
+CONDENSE_FUNCTIONS(BlockVector<double>,double,BlockSparseMatrix<double>);
+CONDENSE_FUNCTIONS(BlockVector<float>,double,BlockSparseMatrix<double>);
+
 
 template
 void
@@ -2056,28 +2078,41 @@ void
 ConstraintMatrix::condense<float>(BlockSparseMatrix<float> &uncondensed) const;
 
 
-#define MATRIX_FUNCTIONS(MatrixType) \
+#define MATRIX_FUNCTIONS(MatrixType, VectorType)       \
 template void ConstraintMatrix:: \
 distribute_local_to_global<MatrixType > (const FullMatrix<double>        &, \
                                          const std::vector<unsigned int> &, \
-                                         MatrixType                      &) const
+                                         MatrixType                      &) const; \
+template void ConstraintMatrix:: \
+distribute_local_to_global<MatrixType,VectorType > (const FullMatrix<double>        &, \
+                                                   const Vector<double>            &, \
+                                                   const std::vector<unsigned int> &, \
+                                                   MatrixType                      &, \
+                                                   VectorType                      &) const
 
-MATRIX_FUNCTIONS(SparseMatrix<double>);
-MATRIX_FUNCTIONS(SparseMatrix<float>);
+MATRIX_FUNCTIONS(SparseMatrix<double>, Vector<double>);
+MATRIX_FUNCTIONS(SparseMatrix<float>, Vector<float>);
 
-MATRIX_FUNCTIONS(BlockSparseMatrix<double>);
-MATRIX_FUNCTIONS(BlockSparseMatrix<float>);
+MATRIX_FUNCTIONS(BlockSparseMatrix<double>, BlockVector<double>);
+MATRIX_FUNCTIONS(BlockSparseMatrix<float>, BlockVector<float>);
 
 #ifdef DEAL_II_USE_PETSC
-MATRIX_FUNCTIONS(PETScWrappers::SparseMatrix);
-MATRIX_FUNCTIONS(PETScWrappers::BlockSparseMatrix);
-MATRIX_FUNCTIONS(PETScWrappers::MPI::SparseMatrix);
-MATRIX_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix);
+MATRIX_FUNCTIONS(PETScWrappers::SparseMatrix, PETScWrappers::Vector);
+MATRIX_FUNCTIONS(PETScWrappers::BlockSparseMatrix, PETScWrappers::BlockVector);
+MATRIX_FUNCTIONS(PETScWrappers::MPI::SparseMatrix, PETScWrappers::MPI::Vector);
+MATRIX_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix ,PETScWrappers::MPI::BlockVector);
 #endif
 
 #ifdef DEAL_II_USE_TRILINOS
-MATRIX_FUNCTIONS(TrilinosWrappers::SparseMatrix);
-MATRIX_FUNCTIONS(TrilinosWrappers::BlockSparseMatrix);
+MATRIX_FUNCTIONS(TrilinosWrappers::SparseMatrix, TrilinosWrappers::VectorBase);
+MATRIX_FUNCTIONS(TrilinosWrappers::BlockSparseMatrix, TrilinosWrappers::BlockVector);
+template void ConstraintMatrix::distribute_local_to_global
+<TrilinosWrappers::BlockSparseMatrix,TrilinosWrappers::MPI::BlockVector> 
+  (const FullMatrix<double>        &, 
+   const Vector<double>            &, 
+   const std::vector<unsigned int> &, 
+   TrilinosWrappers::BlockSparseMatrix &, 
+   TrilinosWrappers::MPI::BlockVector  &) const;
 #endif
 
 template void ConstraintMatrix::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.