]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement a function that can distribute vectors with inhomogeneous constraints,...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 29 Apr 2009 16:31:59 +0000 (16:31 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 29 Apr 2009 16:31:59 +0000 (16:31 +0000)
git-svn-id: https://svn.dealii.org/trunk@18790 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/constraint_matrix.h
deal.II/lac/include/lac/constraint_matrix.templates.h
deal.II/lac/source/constraint_matrix.cc

index b26c2b26163da42d779f7ed515b88d40b0b37bbb..0d0027a4e2468784b7268c43102aa70a4aae323b 100644 (file)
@@ -928,12 +928,7 @@ class ConstraintMatrix : public Subscriptor
                                       * 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.
+                                      * given.
                                       *
                                       * In contrast to the similar function
                                       * in the DoFAccessor class, this
@@ -956,32 +951,47 @@ 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
-                                     * 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.
+                                     * Note that, unless an optional
+                                      * FullMatrix object is provided, this
+                                      * function will apply all constraints
+                                      * as if they were homogeneous and
+                                      * throw an exception in case it
+                                      * encounters inhomogeneous
+                                      * constraints. For correctly setting
+                                      * inhomogeneous constraints, you
+                                      * should provide an additional matrix
+                                      * argument or use one of the functions
+                                      * with both matrix and vector
+                                      * arguments.
+                                     *
+                                     * The optional argument
+                                     * <tt>local_matrix</tt> is intended to
+                                     * be used in case one wants to apply
+                                     * inhomogeneous constraints on the
+                                     * vector only. Such a situation could
+                                     * be where one wants to assemble of a
+                                     * right hand side vector on a problem
+                                     * with inhomogeneous constraints, but
+                                     * the global matrix has been assembled
+                                     * previously. A typical example of
+                                     * this is a time stepping algorithm
+                                     * where the stiffness matrix is
+                                     * assembled once, and the right hand
+                                     * side updated every time step. Note
+                                     * that, however, the entries in the
+                                     * columns of the local matrix have to
+                                     * be exactly the same as those that
+                                     * have been written into the global
+                                     * matrix. Otherwise, this function
+                                     * will not be able to correctly handle
+                                     * inhomogeneities.
                                       */
     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;
+                                VectorType                      &global_vector,
+                               const FullMatrix<double>        &local_matrix = FullMatrix<double>()) const;
 
                                      /**
                                       * This function takes a matrix of
index e5b608499450407c426ef90fe01c0a4a58c88c2c..10e10a3b34d6d4824b834d4846bc0aaf5125a05b 100644 (file)
@@ -727,35 +727,39 @@ void
 ConstraintMatrix::
 distribute_local_to_global (const Vector<double>            &local_vector,
                             const std::vector<unsigned int> &local_dof_indices,
-                            VectorType                      &global_vector) const
+                            VectorType                      &global_vector,
+                           const FullMatrix<double>        &local_matrix) const
 {
   Assert (local_vector.size() == local_dof_indices.size(),
           ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
   Assert (sorted == true, ExcMatrixNotClosed());
+  const bool use_matrix = local_matrix.m() != 0 ? true : false;
+  if (local_matrix.m() != 0)
+    {
+      Assert (local_matrix.m() == local_dof_indices.size(),
+             ExcDimensionMismatch(local_matrix.m(), local_dof_indices.size()));
+      Assert (local_matrix.n() == local_dof_indices.size(),
+             ExcDimensionMismatch(local_matrix.n(), local_dof_indices.size()));
+    }
 
   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
   Threads::ThreadMutex::ScopedLock lock(mutex);
   for (unsigned int i=0; i<n_local_dofs; ++i)
     {
                                   // let's see if we can use the bool
                                   // vector that tells about whether a
-                                  // certain constraint exists.
-      if (constraint_line_exists.size() <= local_dof_indices[i])
-       {
-         global_vector(local_dof_indices[i]) += local_vector(i);
-         continue;
-       }
-      if (constraint_line_exists[local_dof_indices[i]] == false)
+                                  // certain constraint exists. then, we
+                                  // simply copy over the data.
+      if (is_constrained (local_dof_indices[i]) == false)
        {
          global_vector(local_dof_indices[i]) += local_vector(i);
          continue;
        }
-                                           // first figure out whether this
-                                           // dof is constrained
+
+                                  // if the dof is constrained, we
+                                  // have to find it in the list of
+                                  // all constraints.
       ConstraintLine index_comparison;
       index_comparison.line = local_dof_indices[i];
 
@@ -764,40 +768,61 @@ distribute_local_to_global (const Vector<double>            &local_vector,
                                     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
+                                  // check whether we've really found the
+                                  // right constraint.
       Assert (position->line == local_dof_indices[i],
              ExcInternalError());
+
+      if (use_matrix)
+       {
+         const double val = position->inhomogeneity;
+         if (val != 0)
+           for (unsigned int j=0; j<n_local_dofs; ++j)
+             if (is_constrained (local_dof_indices[j]) == false)
+               {
+                 global_vector(local_dof_indices[j]) -= val * local_matrix(j,i);
+               }
+             else
+               {
+                 const double matrix_entry = local_matrix(j,i);
+                 if (matrix_entry == 0)
+                   continue;
+
+                 ConstraintLine index_comparison_j;
+                 index_comparison_j.line = local_dof_indices[j];
+
+                 const std::vector<ConstraintLine>::const_iterator
+                   position_j = std::lower_bound (lines.begin(),
+                                                  lines.end(),
+                                                  index_comparison_j);
+                 for (unsigned int q=0; q<position_j->entries.size(); ++q)
+                   {
+                     Assert (is_constrained(position_j->entries[q].first) == false,
+                             ExcMessage ("Tried to distribute to a fixed dof."));
+                     global_vector(position_j->entries[j].first)
+                       -= val * position_j->entries[q].second * matrix_entry;
+                   }
+               }
+       }
+      else
+                                  // in case the constraint is
+                                  // inhomogeneous and we have no matrix
+                                  // available, this function is not
+                                  // appropriate. Throw an exception.
+       Assert (position->inhomogeneity == 0.,
+               ExcMessage ("Inhomogeneous constraint cannot be condensed "
+                           "without any matrix specified."));
+
+                                  // now distribute the constraint,
+                                  // but make sure we don't touch
+                                  // the entries of fixed dofs
       for (unsigned int j=0; j<position->entries.size(); ++j)
-       global_vector(position->entries[j].first)
-         += local_vector(i) * position->entries[j].second;
+       {
+         Assert (is_constrained(position->entries[j].first) == false,
+                 ExcMessage ("Tried to distribute to a fixed dof."));
+         global_vector(position->entries[j].first)
+           += local_vector(i) * position->entries[j].second;
+       }
     }
 }
 
index 7d165e7dc5a2dd1d11ad68a2bd92837655dcebac..d5d26297f81875f54705882f70184b4bb3cc0f58 100644 (file)
@@ -2004,7 +2004,8 @@ ConstraintMatrix::memory_consumption () const
   template void ConstraintMatrix:: \
     distribute_local_to_global<VectorType > (const Vector<double>            &, \
                                              const std::vector<unsigned int> &, \
-                                             VectorType                      &) const; \
+                                             VectorType                      &, \
+                                            const FullMatrix<double>        &) const; \
   template void ConstraintMatrix::distribute<VectorType >(const VectorType &condensed,\
                                                         VectorType       &uncondensed) const;\
   template void ConstraintMatrix::distribute<VectorType >(VectorType &vec) const

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.