From 29fbe90ab964bb9bf214d2976a62eb6eacd8c360 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Sat, 21 Feb 2009 16:12:26 +0000 Subject: [PATCH] Complete the implementation of inhomogeneous constraints in the ConstraintMatrix class. Need to write a function in VectorTools in order to enable the direct use with Dirichlet b.c. git-svn-id: https://svn.dealii.org/trunk@18403 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/dofs/dof_constraints.h | 733 ++++++------ .../include/dofs/dof_constraints.templates.h | 1008 +++++++++++++++-- .../deal.II/source/dofs/dof_constraints.cc | 59 +- 3 files changed, 1342 insertions(+), 458 deletions(-) diff --git a/deal.II/deal.II/include/dofs/dof_constraints.h b/deal.II/deal.II/include/dofs/dof_constraints.h index ea3ec36ac8..bb0c96144f 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.h @@ -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. + * + * + *

Using the ConstraintMatrix for hanging nodes

+ * + * For example, when using Q1 and Q2 elements (i.e. using * FE_Q<dim,spacedim>(1) and FE_Q<dim,spacedim>(2)) on the two * marked cells of the mesh * @@ -69,20 +75,32 @@ class BlockIndices; * some detail in the @ref hp_paper "hp paper". * * + *

Using the ConstraintMatrix for Dirichlet boundary conditions

+ * + * 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. + * + * *

Description of constraints

* - * 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. * * *
Condensing vectors
@@ -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. + * + * + *
Treatment of inhomogeneous constraints
+ * + * 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). * * *

Avoiding explicit condensation

@@ -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". * - *
  • - * 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". + *
  • 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". * * * In this case, one possibility is to distribute local entries to the final @@ -258,16 +288,16 @@ class BlockIndices; * *

    Distributing constraints

    * - * 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 &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, Vector, @@ -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 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 + void condense (const SparseMatrix &uncondensed_matrix, + const VectorType &uncondensed_vector, + SparseMatrix &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 + void condense (SparseMatrix &matrix, + VectorType &vector) const; + + /** + * Same function as above, but + * condenses square block sparse + * matrices and vectors. + */ + template + void condense (BlockSparseMatrix &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 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 @@ -986,6 +1073,23 @@ class ConstraintMatrix : public Subscriptor const std::vector &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 + void + distribute_local_to_global (const FullMatrix &local_matrix, + const Vector &local_vector, + const std::vector &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, Vector, - * BlockVector<...>, a - * PETSc or Trilinos vector - * wrapper class, or any other - * type having the same + * BlockVector<...>, a PETSc + * or Trilinos vector wrapper class, or + * any other type having the same * interface. */ template @@ -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, Vector, - * BlockVector<...>, 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, + * Vector, + * BlockVector<...>, a PETSc + * or Trilinos vector wrapper class, or + * any other type having the same * interface. */ template @@ -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 > 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 close() 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 + * close() 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 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 &p); /** - * Dummy table that serves as - * default argument for function + * Dummy table that serves as default + * argument for function * add_entries_local_to_global(). */ static const Table<2,bool> default_empty_table; @@ -1506,7 +1595,7 @@ ConstraintMatrix::set_inhomogeneity (const unsigned int line, const std::vector::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) diff --git a/deal.II/deal.II/include/dofs/dof_constraints.templates.h b/deal.II/deal.II/include/dofs/dof_constraints.templates.h index 69c40bb38c..60c4b5716d 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.templates.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.templates.h @@ -145,7 +145,7 @@ ConstraintMatrix::condense (const SparseMatrix &uncondensed, next_constraint->entries[q].second * c->entries[p].second); }; - + ++next_constraint; }; } @@ -260,7 +260,7 @@ ConstraintMatrix::condense (SparseMatrix &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 +template void -ConstraintMatrix::set_zero (VectorType &vec) const +ConstraintMatrix::condense (const SparseMatrix &uncondensed, + const VectorType &uncondensed_vector, + SparseMatrix &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::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 new_line; + new_line.reserve (uncondensed_struct.n_rows()); -template -void -ConstraintMatrix:: -distribute_local_to_global (const Vector &local_vector, - const std::vector &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::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; iline) + { + // 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::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; rowentries.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::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]; + jentries.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::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 +template void -ConstraintMatrix:: -distribute_local_to_global (const FullMatrix &local_matrix, - const std::vector &local_dof_indices, - MatrixType &global_matrix) const +ConstraintMatrix::condense (SparseMatrix &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 distribute (sparsity.n_rows(), + numbers::invalid_unsigned_int); + + for (unsigned int c=0; c - constraint_lines (n_local_dofs, - static_cast(0)); - unsigned int n_max_entries_per_row = 0; - for (unsigned int i=0; i::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::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 +void +ConstraintMatrix::condense (BlockSparseMatrix &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 distribute (sparsity.n_rows(), + numbers::invalid_unsigned_int); + + for (unsigned int c=0; c + 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::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::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 +void +ConstraintMatrix::set_zero (VectorType &vec) const +{ + Assert (sorted == true, ExcMatrixNotClosed()); + + std::vector::const_iterator constraint_line = lines.begin(); + for (; constraint_line!=lines.end(); ++constraint_line) + vec(constraint_line->line) = 0.; +} + + + +template +void +ConstraintMatrix:: +distribute_local_to_global (const Vector &local_vector, + const std::vector &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::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; jentries.size(); ++j) + global_vector(position->entries[j].first) + += local_vector(i) * position->entries[j].second; + } + } + } +} + + + +template +void +ConstraintMatrix:: +distribute_local_to_global (const FullMatrix &local_matrix, + const std::vector &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 + constraint_lines (n_local_dofs, + static_cast(0)); + unsigned int n_max_entries_per_row = 0; + for (unsigned int i=0; i::const_iterator position = std::lower_bound (lines.begin(), @@ -868,7 +1362,273 @@ distribute_local_to_global (const FullMatrix &local_matrix, false); } } - mutex.release(); +} + + + +template +void +ConstraintMatrix:: +distribute_local_to_global (const FullMatrix &local_matrix, + const Vector &local_vector, + const std::vector &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 + constraint_lines (n_local_dofs, + static_cast(0)); + unsigned int n_max_entries_per_row = 0; + for (unsigned int i=0; i::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; ientries.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; qentries.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; pentries.size(); ++p) + for (unsigned int q=0; qentries.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; qentries.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); + } + } } diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index 846b9df307..46e2b50137 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -1997,6 +1997,14 @@ ConstraintMatrix::memory_consumption () const template void ConstraintMatrix::condense(const VectorType &uncondensed,\ VectorType &condensed) const;\ template void ConstraintMatrix::condense(VectorType &vec) const;\ + template void ConstraintMatrix::condense(const SparseMatrix &uncondensed, \ + const VectorType &uncondensed_vector, \ + SparseMatrix &condensed, \ + VectorType &condensed_vector) const; \ + template void ConstraintMatrix::condense(const SparseMatrix &uncondensed, \ + const VectorType &uncondensed_vector, \ + SparseMatrix &condensed, \ + VectorType &condensed_vector) const; \ template void ConstraintMatrix::set_zero(VectorType &vec) const;\ template void ConstraintMatrix:: \ distribute_local_to_global (const Vector &, \ @@ -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(MatrixType &uncondensed, \ + VectorType &vec) const \ + + +CONDENSE_FUNCTIONS(Vector,float,SparseMatrix); +CONDENSE_FUNCTIONS(Vector,float,SparseMatrix); +CONDENSE_FUNCTIONS(Vector,double,SparseMatrix); +CONDENSE_FUNCTIONS(Vector,double,SparseMatrix); +CONDENSE_FUNCTIONS(BlockVector,float,BlockSparseMatrix); +CONDENSE_FUNCTIONS(BlockVector,float,BlockSparseMatrix); +CONDENSE_FUNCTIONS(BlockVector,double,BlockSparseMatrix); +CONDENSE_FUNCTIONS(BlockVector,double,BlockSparseMatrix); + template void @@ -2056,28 +2078,41 @@ void ConstraintMatrix::condense(BlockSparseMatrix &uncondensed) const; -#define MATRIX_FUNCTIONS(MatrixType) \ +#define MATRIX_FUNCTIONS(MatrixType, VectorType) \ template void ConstraintMatrix:: \ distribute_local_to_global (const FullMatrix &, \ const std::vector &, \ - MatrixType &) const + MatrixType &) const; \ +template void ConstraintMatrix:: \ +distribute_local_to_global (const FullMatrix &, \ + const Vector &, \ + const std::vector &, \ + MatrixType &, \ + VectorType &) const -MATRIX_FUNCTIONS(SparseMatrix); -MATRIX_FUNCTIONS(SparseMatrix); +MATRIX_FUNCTIONS(SparseMatrix, Vector); +MATRIX_FUNCTIONS(SparseMatrix, Vector); -MATRIX_FUNCTIONS(BlockSparseMatrix); -MATRIX_FUNCTIONS(BlockSparseMatrix); +MATRIX_FUNCTIONS(BlockSparseMatrix, BlockVector); +MATRIX_FUNCTIONS(BlockSparseMatrix, BlockVector); #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 + + (const FullMatrix &, + const Vector &, + const std::vector &, + TrilinosWrappers::BlockSparseMatrix &, + TrilinosWrappers::MPI::BlockVector &) const; #endif template void ConstraintMatrix:: -- 2.39.5