]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add distribute_local_to_global functions.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 15 Mar 2004 16:40:48 +0000 (16:40 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 15 Mar 2004 16:40:48 +0000 (16:40 +0000)
git-svn-id: https://svn.dealii.org/trunk@8742 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 cd4b7a687e60105c5927393290bce9b314812f87..6efc8fd2e4fb3397979fc37dc276142ff02a6b46 100644 (file)
@@ -20,6 +20,8 @@
 #include <base/exceptions.h>
 #include <base/subscriptor.h>
 
+template <typename> class Vector;
+template <typename> 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 <class VectorType>
     void set_zero (VectorType &vec) const;
 
-
+                                     /**
+                                      * This function takes a vector of local
+                                      * contributions (@arg local_vector)
+                                      * corresponding to the degrees of
+                                      * freedom indices given in @arg
+                                      * local_dof_indices and distributes them
+                                      * to the global vector. In most cases,
+                                      * these local contributions will be the
+                                      * result of an integration over a cell
+                                      * or face of a cell. However, as long as
+                                      * @arg local_vector and @arg
+                                      * local_dof_indices have the same number
+                                      * of elements, this function is happy
+                                      * with whatever it is given.
+                                      *
+                                      * In contrast to the similar function in
+                                      * the DoFAccessor class, this function
+                                      * also takes care of constraints,
+                                      * i.e. of one of the elements of @arg
+                                      * local_dof_indices belongs to a
+                                      * constrained node, then rather than
+                                      * writing the corresponding element of
+                                      * @arg local_vector into @arg
+                                      * global_vector, the element is
+                                      * distributed to the entries in the
+                                      * global vector to which this particular
+                                      * degree of freedom is constrained.
+                                      *
+                                      * Thus, by using this function to
+                                      * distribute local contributions to the
+                                      * global object, one saves the call to
+                                      * the condense function after the
+                                      * vectors and matrices are fully
+                                      * assembled.
+                                      */
+    template <typename VectorType>
+    void
+    distribute_local_to_global (const Vector<double>            &local_vector,
+                                const std::vector<unsigned int> &local_dof_indices,
+                                VectorType                      &global_vector) const;
+
+                                     /**
+                                      * This function takes a matrix of local
+                                      * contributions (@arg local_matrix)
+                                      * corresponding to the degrees of
+                                      * freedom indices given in @arg
+                                      * local_dof_indices and distributes them
+                                      * to the global matrix. In most cases,
+                                      * these local contributions will be the
+                                      * result of an integration over a cell
+                                      * or face of a cell. However, as long as
+                                      * @arg local_matrix and @arg
+                                      * local_dof_indices have the same number
+                                      * of elements, this function is happy
+                                      * with whatever it is given.
+                                      *
+                                      * In contrast to the similar function in
+                                      * the DoFAccessor class, this function
+                                      * also takes care of constraints,
+                                      * i.e. of one of the elements of @arg
+                                      * local_dof_indices belongs to a
+                                      * constrained node, then rather than
+                                      * writing the corresponding element of
+                                      * @arg local_matrix into @arg
+                                      * global_matrix, the element is
+                                      * distributed to the entries in the
+                                      * global matrix to which this particular
+                                      * degree of freedom is constrained.
+                                      *
+                                      * With this scheme, we never write into
+                                      * rows or columns of constrained degrees
+                                      * of freedom. In order to make sure that
+                                      * the resulting matrix can still be
+                                      * inverted, we need to do something with
+                                      * the diagonal elements corresponding to
+                                      * constrained nodes. Thus, if a degree
+                                      * of freedom in @arg local_dof_indices
+                                      * is constrained, we distribute the
+                                      * corresponding entries in the matrix,
+                                      * but also add the absolute value of the
+                                      * diagonal entry of the local matrix to
+                                      * the corresponding entry in the global
+                                      * matrix. Since the exact value of the
+                                      * diagonal element is not important (the
+                                      * value of the respective degree of
+                                      * freedom will be overwritten by the
+                                      * distribute() call later on anyway),
+                                      * this guarantees that the diagonal
+                                      * entry is always non-zero, positive,
+                                      * and of the same order of magnitude as
+                                      * the other entries of the matrix.
+                                      *
+                                      * Thus, by using this function to
+                                      * distribute local contributions to the
+                                      * global object, one saves the call to
+                                      * the condense function after the
+                                      * vectors and matrices are fully
+                                      * assembled.
+                                      */
+    template <typename MatrixType>
+    void
+    distribute_local_to_global (const FullMatrix<double>        &local_matrix,
+                                const std::vector<unsigned int> &local_dof_indices,
+                                MatrixType                      &global_matrix) const;
+    
                                     /**
                                      * Print the constraint lines. Mainly for
                                      * debugging purposes.
index 1d94667def95b6030ca6f33ab7b2a890d5b08a4f..fd6938e7b6666d06cee029c0aa2f2b74b33a718f 100644 (file)
@@ -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 <base/config.h>
 #include <dofs/dof_constraints.h>
+#include <lac/full_matrix.h>
 #include <lac/sparsity_pattern.h>
 #include <lac/sparse_matrix.h>
 #include <lac/block_sparsity_pattern.h>
@@ -624,6 +625,202 @@ ConstraintMatrix::set_zero (VectorType &vec) const
 
 
 
+template <typename VectorType>
+void
+ConstraintMatrix::
+distribute_local_to_global (const Vector<double>            &local_vector,
+                            const std::vector<unsigned int> &local_dof_indices,
+                            VectorType                      &global_vector) const
+{
+  Assert (local_vector.size() == local_dof_indices.size(),
+          ExcWrongDimension());
+  Assert (sorted == true, ExcMatrixNotClosed());
+
+  const unsigned int n_local_dofs = local_vector.size();
+  
+                                   // have a special case where there are no
+                                   // constraints at all, since then we can be
+                                   // a lot faster
+  if (lines.size() == 0)
+    {
+      for (unsigned int i=0; i<n_local_dofs; ++i)
+        global_vector(local_dof_indices[i]) += local_vector(i);
+    }
+  else
+    {
+      for (unsigned int i=0; i<n_local_dofs; ++i)
+        {
+                                           // first figure out whether this
+                                           // dof is constrained
+          ConstraintLine index_comparison;
+          index_comparison.line = local_dof_indices[i];
+
+          const std::vector<ConstraintLine>::const_iterator
+            position = std::lower_bound (lines.begin(),
+                                         lines.end(),
+                                         index_comparison);
+
+                                           // if the line is not constrained,
+                                           // then simply copy the
+                                           // data. otherwise distribute it
+          if (position->line != local_dof_indices[i])
+            global_vector(local_dof_indices[i]) += local_vector(i);
+          else
+            for (unsigned int j=0; j<position->entries.size(); ++j)
+              global_vector(position->entries[j].first)
+                += local_vector(i) * position->entries[j].second;
+        }
+    }
+}
+
+
+
+template <typename MatrixType>
+void
+ConstraintMatrix::
+distribute_local_to_global (const FullMatrix<double>        &local_matrix,
+                            const std::vector<unsigned int> &local_dof_indices,
+                            MatrixType                      &global_matrix) const
+{
+  Assert (local_matrix.n() == local_dof_indices.size(),
+          ExcWrongDimension());
+  Assert (local_matrix.m() == local_dof_indices.size(),
+          ExcWrongDimension());
+  Assert (sorted == true, ExcMatrixNotClosed());
+
+  const unsigned int n_local_dofs = local_dof_indices.size();
+  
+                                   // have a special case where there are no
+                                   // constraints at all, since then we can be
+                                   // a lot faster
+  if (lines.size() == 0)
+    {
+      for (unsigned int i=0; i<n_local_dofs; ++i)
+        for (unsigned int j=0; j<n_local_dofs; ++j)
+          global_matrix.add(local_dof_indices[i],
+                            local_dof_indices[j],
+                            local_matrix(i,j));
+    }
+  else
+    {
+                                       // here we have to do something a
+                                       // little nastier than in the
+                                       // respective function for vectors. the
+                                       // reason is that we have two nested
+                                       // loops and we don't want to
+                                       // repeatedly check whether a certain
+                                       // dof is constrained or not by
+                                       // searching over all the constrained
+                                       // dofs. so we have to cache this
+                                       // knowledge, by storing for each dof
+                                       // index whether and where the line of
+                                       // the constraint matrix is located
+      std::vector<const ConstraintLine *>
+        constraint_lines (n_local_dofs,
+                          static_cast<const ConstraintLine *>(0));
+      for (unsigned int i=0; i<n_local_dofs; ++i)
+        {
+          ConstraintLine index_comparison;
+          index_comparison.line = local_dof_indices[i];
+
+          const std::vector<ConstraintLine>::const_iterator
+            position = std::lower_bound (lines.begin(),
+                                         lines.end(),
+                                         index_comparison);
+          
+                                           // if this dof is constrained, then
+                                           // set the respective entry in the
+                                           // array. otherwise leave it at the
+                                           // invalid position
+          if ((position != lines.end()) &&
+              (position->line == local_dof_indices[i]))
+            constraint_lines[i] = &*position;
+        }
+
+
+                                       // now distribute entries
+      for (unsigned int i=0; i<n_local_dofs; ++i)
+        {
+          const ConstraintLine *position_i = constraint_lines[i];
+          const bool is_constrained_i = (position_i != 0);
+          
+          for (unsigned int j=0; j<n_local_dofs; ++j)
+            {
+              const ConstraintLine *position_j = constraint_lines[j];
+              const bool is_constrained_j = (position_j != 0);
+
+              if ((is_constrained_i == false) &&
+                  (is_constrained_j == false))
+                {
+                                                   // neither row nor column
+                                                   // is constrained
+                  global_matrix.add (local_dof_indices[i],
+                                     local_dof_indices[j],
+                                     local_matrix(i,j));
+                }
+              else if ((is_constrained_i == true) &&
+                       (is_constrained_j == false))
+                {
+                                                   // ok, row is constrained,
+                                                   // but column is not
+                  for (unsigned int q=0; q<position_i->entries.size(); ++q)
+                    global_matrix.add (position_i->entries[q].first,
+                                       local_dof_indices[j],
+                                       local_matrix(i,j) *
+                                       position_i->entries[q].second);
+                }
+              else if ((is_constrained_i == false) &&
+                       (is_constrained_j == true))
+                {
+                                                   // simply the other way
+                                                   // round: row ok, column is
+                                                   // constrained
+                  for (unsigned int q=0; q<position_j->entries.size(); ++q)
+                    global_matrix.add (local_dof_indices[i],
+                                       position_j->entries[q].first,
+                                       local_matrix(i,j) *
+                                       position_j->entries[q].second);
+                }
+              else if ((is_constrained_i == true) &&
+                       (is_constrained_j == true))
+                {
+                                                   // last case: both row and
+                                                   // column are constrained
+                  for (unsigned int p=0; p<position_i->entries.size(); ++p)
+                    for (unsigned int q=0; q<position_j->entries.size(); ++q)
+                      global_matrix.add (position_i->entries[p].first,
+                                         position_j->entries[q].first,
+                                         local_matrix(i,j) *
+                                         position_i->entries[p].second *
+                                         position_j->entries[q].second);
+
+                                                   // to make sure that the
+                                                   // global matrix remains
+                                                   // invertible, we need to
+                                                   // do something with the
+                                                   // diagonal elements. add
+                                                   // the absolute value of
+                                                   // the local matrix, so the
+                                                   // resulting entry will
+                                                   // always be positive and
+                                                   // furthermore be in the
+                                                   // same order of magnitude
+                                                   // as the other elements of
+                                                   // the matrix
+                  if (i == j)
+                    global_matrix.add (local_dof_indices[i],
+                                       local_dof_indices[i],
+                                       local_matrix(i,i));
+                }
+              else
+                Assert (false, ExcInternalError());
+            }
+        }
+    }
+}
+
+
+
 template<class VectorType>
 void
 ConstraintMatrix::distribute (const VectorType &condensed,
index a3de04de20d80359c474ff419e466629204b4a46..789d3b8b60007e5ae5d512a55a1edf35374074e2 100644 (file)
@@ -20,6 +20,7 @@
 #include <lac/compressed_sparsity_pattern.h>
 #include <lac/vector.h>
 #include <lac/petsc_vector.h>
+#include <lac/petsc_sparse_matrix.h>
 #include <lac/block_vector.h>
 #include <lac/sparse_matrix.h>
 #include <lac/block_sparse_matrix.h>
@@ -1253,6 +1254,10 @@ ConstraintMatrix::memory_consumption () const
                                                       VectorType       &condensed) const;\
   template void ConstraintMatrix::condense<VectorType >(VectorType &vec) const;\
   template void ConstraintMatrix::set_zero<VectorType >(VectorType &vec) const;\
+  template void ConstraintMatrix:: \
+    distribute_local_to_global<VectorType> (const Vector<double>            &, \
+                                            const std::vector<unsigned int> &, \
+                                            VectorType                      &) const; \
   template void ConstraintMatrix::distribute<VectorType >(const VectorType &condensed,\
                                                         VectorType       &uncondensed) const;\
   template void ConstraintMatrix::distribute<VectorType >(VectorType &vec) const
@@ -1295,3 +1300,20 @@ ConstraintMatrix::condense<double>(BlockSparseMatrix<double> &uncondensed) const
 template
 void
 ConstraintMatrix::condense<float>(BlockSparseMatrix<float> &uncondensed) const;
+
+
+#define MATRIX_FUNCTIONS(MatrixType) \
+template void ConstraintMatrix:: \
+distribute_local_to_global<MatrixType > (const FullMatrix<double>        &, \
+                                         const std::vector<unsigned int> &, \
+                                         MatrixType                      &) const
+
+MATRIX_FUNCTIONS(SparseMatrix<double>);
+MATRIX_FUNCTIONS(SparseMatrix<float>);
+
+MATRIX_FUNCTIONS(BlockSparseMatrix<double>);
+MATRIX_FUNCTIONS(BlockSparseMatrix<float>);
+
+#ifdef DEAL_II_USE_PETSC
+MATRIX_FUNCTIONS(PETScWrappers::SparseMatrix);
+#endif

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.