From 8ac3495c4fa71188dcd563bbce967a5581cdecd3 Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 15 Mar 2004 16:40:48 +0000 Subject: [PATCH] Add distribute_local_to_global functions. git-svn-id: https://svn.dealii.org/trunk@8742 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/dofs/dof_constraints.h | 149 +++++++++++-- .../include/dofs/dof_constraints.templates.h | 199 +++++++++++++++++- .../deal.II/source/dofs/dof_constraints.cc | 22 ++ 3 files changed, 356 insertions(+), 14 deletions(-) diff --git a/deal.II/deal.II/include/dofs/dof_constraints.h b/deal.II/deal.II/include/dofs/dof_constraints.h index cd4b7a687e..6efc8fd2e4 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.h @@ -20,6 +20,8 @@ #include #include +template class Vector; +template class FullMatrix; class SparsityPattern; class CompressedSparsityPattern; class BlockSparsityPattern; @@ -80,17 +82,20 @@ class BlockIndices; * those which are subject to constraints. * * @item Use only one sparsity pattern and one matrix: doing it this way, the - * condense functions add nonzero entries to the sparsity pattern of the - * large matrix (with constrained nodes in it) where the condensation process - * of the matrix will create additional nonzero elements. In the condensation - * process itself, lines and rows subject to constraints are distributed to - * the lines and rows of unconstrained nodes. The constrained lines remain in - * place, however, unlike in the first possibility described above. In order - * not to disturb the solution process, these lines and rows are filled with - * zeros and identity on the main diagonal; the appropriate value in the right - * hand sides is set to zero. This way, the constrained node will always get - * the value zero upon solution of the equation system and will not couple to - * other nodes any more. + * condense functions add nonzero entries to the sparsity pattern of the large + * matrix (with constrained nodes in it) where the condensation process of the + * matrix will create additional nonzero elements. In the condensation process + * itself, lines and rows subject to constraints are distributed to the lines + * and rows of unconstrained nodes. The constrained lines remain in place, + * however, unlike in the first possibility described above. In order not to + * disturb the solution process, these lines and rows are filled with zeros + * and an appropriate positive value on the main diagonal (we choose an + * average of the magnitudes of the other diagonal elements, so as to make + * sure that the new diagonal entry has the same order of magnitude as the + * other entries; this preserves the scaling properties of the matrix). The + * appropriate value in the right hand sides is set to zero. This way, the + * constrained node will always get the value zero upon solution of the + * equation system and will not couple to other nodes any more. * * This method has the advantage that only one matrix and sparsity pattern is * needed thus using less memory. Additionally, the condensation process is @@ -138,8 +143,22 @@ class BlockIndices; * functions exist in variants for PETSc vectors. However, using them is * typically expensive, and should be avoided. You should use the same * techniques as mentioned above to avoid their use. + * * + * @sect3{Avoiding explicit condensation} + * + * Sometimes, one wants to avoid condensation at all. This may be the case + * since condensation is an expensive operation, or because no condense() + * function is defined for the matrix you use (this is, for example, the case + * for the PETSc wrapper classes, where we have no access to the underlying + * representation of the matrix, and therefore cannot efficiently implement + * the condense() operation). In this case, one possibility is to distribute + * local entries to the final destinations right at the moment of transferring + * them into the global matrices and vectors. For this, one can use the + * distribute_local_to_global() functions of this class, which make a + * subsequent call to condense() unnecessary. * + * * @sect3{Distributing constraints} * * After solving the condensed system of equations, the solution vector has to @@ -153,7 +172,7 @@ class BlockIndices; * for unconstrained nodes, and constrained nodes need to get their values in * a second step. * - * @author Wolfgang Bangerth, 1998 + * @author Wolfgang Bangerth, 1998, 2004 */ class ConstraintMatrix : public Subscriptor { @@ -598,7 +617,111 @@ class ConstraintMatrix : public Subscriptor template void set_zero (VectorType &vec) const; - + /** + * This function takes a vector of local + * contributions (@arg local_vector) + * corresponding to the degrees of + * freedom indices given in @arg + * local_dof_indices and distributes them + * to the global vector. In most cases, + * these local contributions will be the + * result of an integration over a cell + * or face of a cell. However, as long as + * @arg local_vector and @arg + * local_dof_indices have the same number + * of elements, this function is happy + * with whatever it is given. + * + * In contrast to the similar function in + * the DoFAccessor class, this function + * also takes care of constraints, + * i.e. of one of the elements of @arg + * local_dof_indices belongs to a + * constrained node, then rather than + * writing the corresponding element of + * @arg local_vector into @arg + * global_vector, the element is + * distributed to the entries in the + * global vector to which this particular + * degree of freedom is constrained. + * + * Thus, by using this function to + * distribute local contributions to the + * global object, one saves the call to + * the condense function after the + * vectors and matrices are fully + * assembled. + */ + template + void + distribute_local_to_global (const Vector &local_vector, + const std::vector &local_dof_indices, + VectorType &global_vector) const; + + /** + * This function takes a matrix of local + * contributions (@arg local_matrix) + * corresponding to the degrees of + * freedom indices given in @arg + * local_dof_indices and distributes them + * to the global matrix. In most cases, + * these local contributions will be the + * result of an integration over a cell + * or face of a cell. However, as long as + * @arg local_matrix and @arg + * local_dof_indices have the same number + * of elements, this function is happy + * with whatever it is given. + * + * In contrast to the similar function in + * the DoFAccessor class, this function + * also takes care of constraints, + * i.e. of one of the elements of @arg + * local_dof_indices belongs to a + * constrained node, then rather than + * writing the corresponding element of + * @arg local_matrix into @arg + * global_matrix, the element is + * distributed to the entries in the + * global matrix to which this particular + * degree of freedom is constrained. + * + * With this scheme, we never write into + * rows or columns of constrained degrees + * of freedom. In order to make sure that + * the resulting matrix can still be + * inverted, we need to do something with + * the diagonal elements corresponding to + * constrained nodes. Thus, if a degree + * of freedom in @arg local_dof_indices + * is constrained, we distribute the + * corresponding entries in the matrix, + * but also add the absolute value of the + * diagonal entry of the local matrix to + * the corresponding entry in the global + * matrix. Since the exact value of the + * diagonal element is not important (the + * value of the respective degree of + * freedom will be overwritten by the + * distribute() call later on anyway), + * this guarantees that the diagonal + * entry is always non-zero, positive, + * and of the same order of magnitude as + * the other entries of the matrix. + * + * Thus, by using this function to + * distribute local contributions to the + * global object, one saves the call to + * the condense function after the + * vectors and matrices are fully + * assembled. + */ + template + void + distribute_local_to_global (const FullMatrix &local_matrix, + const std::vector &local_dof_indices, + MatrixType &global_matrix) const; + /** * Print the constraint lines. Mainly for * debugging purposes. 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 1d94667def..fd6938e7b6 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.templates.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.templates.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -16,6 +16,7 @@ #include #include +#include #include #include #include @@ -624,6 +625,202 @@ ConstraintMatrix::set_zero (VectorType &vec) const +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(), + ExcWrongDimension()); + Assert (sorted == true, ExcMatrixNotClosed()); + + const unsigned int n_local_dofs = local_vector.size(); + + // have a special case where there are no + // constraints at all, since then we can be + // a lot faster + if (lines.size() == 0) + { + for (unsigned int i=0; i::const_iterator + position = std::lower_bound (lines.begin(), + lines.end(), + index_comparison); + + // if the line is not constrained, + // then simply copy the + // data. otherwise distribute it + if (position->line != local_dof_indices[i]) + global_vector(local_dof_indices[i]) += local_vector(i); + else + for (unsigned int j=0; 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(), + ExcWrongDimension()); + Assert (local_matrix.m() == local_dof_indices.size(), + ExcWrongDimension()); + Assert (sorted == true, ExcMatrixNotClosed()); + + const unsigned int n_local_dofs = local_dof_indices.size(); + + // have a special case where there are no + // constraints at all, since then we can be + // a lot faster + if (lines.size() == 0) + { + for (unsigned int i=0; i + constraint_lines (n_local_dofs, + static_cast(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; + } + + + // now distribute entries + 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 + for (unsigned int q=0; qentries.size(); ++q) + global_matrix.add (local_dof_indices[i], + position_j->entries[q].first, + local_matrix(i,j) * + position_j->entries[q].second); + } + else if ((is_constrained_i == true) && + (is_constrained_j == true)) + { + // last case: both row and + // column are constrained + for (unsigned int p=0; 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 + if (i == j) + global_matrix.add (local_dof_indices[i], + local_dof_indices[i], + local_matrix(i,i)); + } + else + Assert (false, ExcInternalError()); + } + } + } +} + + + template void ConstraintMatrix::distribute (const VectorType &condensed, diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index a3de04de20..789d3b8b60 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -1253,6 +1254,10 @@ ConstraintMatrix::memory_consumption () const VectorType &condensed) const;\ template void ConstraintMatrix::condense(VectorType &vec) const;\ template void ConstraintMatrix::set_zero(VectorType &vec) const;\ + template void ConstraintMatrix:: \ + distribute_local_to_global (const Vector &, \ + const std::vector &, \ + VectorType &) const; \ template void ConstraintMatrix::distribute(const VectorType &condensed,\ VectorType &uncondensed) const;\ template void ConstraintMatrix::distribute(VectorType &vec) const @@ -1295,3 +1300,20 @@ ConstraintMatrix::condense(BlockSparseMatrix &uncondensed) const template void ConstraintMatrix::condense(BlockSparseMatrix &uncondensed) const; + + +#define MATRIX_FUNCTIONS(MatrixType) \ +template void ConstraintMatrix:: \ +distribute_local_to_global (const FullMatrix &, \ + const std::vector &, \ + MatrixType &) const + +MATRIX_FUNCTIONS(SparseMatrix); +MATRIX_FUNCTIONS(SparseMatrix); + +MATRIX_FUNCTIONS(BlockSparseMatrix); +MATRIX_FUNCTIONS(BlockSparseMatrix); + +#ifdef DEAL_II_USE_PETSC +MATRIX_FUNCTIONS(PETScWrappers::SparseMatrix); +#endif -- 2.39.5