]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Work on a few TODOs in DoFCellAccessor::distribute_local_to_global. Let these impleme...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 21 Sep 2009 06:48:08 +0000 (06:48 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 21 Sep 2009 06:48:08 +0000 (06:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@19483 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_accessor.templates.h
deal.II/lac/include/lac/block_vector_base.h
deal.II/lac/include/lac/constraint_matrix.templates.h
deal.II/lac/include/lac/sparse_matrix.h
deal.II/lac/include/lac/vector.h

index 26618816a38cf0acfb61dc9c4b21d356d05a2b3a..17ea766c9781b45e72c54e112cf6b05601f02e4b 100644 (file)
@@ -2163,6 +2163,8 @@ namespace internal
              values(*cache) = local_values(i);
          }      
 
+
+
        template <int dim, int spacedim, class OutputVector, typename number>
        static
        void
@@ -2244,7 +2246,7 @@ namespace internal
        static
        void
        set_active_fe_index (DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim> > &accessor,
-                            const unsigned int                              i)
+                            const unsigned int                                      i)
          {
            typedef
              dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
@@ -2261,6 +2263,165 @@ namespace internal
            accessor.dof_handler->levels[accessor.level()]
              ->active_fe_indices[accessor.present_index] = i;
          }
+
+
+
+        template <int dim, int spacedim, typename number, class OutputVector>
+       static
+       void
+       distribute_local_to_global (const DoFCellAccessor<dealii::DoFHandler<dim,spacedim> > &accessor,
+                                   const dealii::Vector<number> &local_source,
+                                   OutputVector                 &global_destination)
+          {
+           typedef
+             dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
+             BaseClass;
+
+           Assert (accessor.dof_handler != 0,
+                   typename BaseClass::ExcInvalidObject());
+           Assert (&accessor.get_fe() != 0,
+                   typename BaseClass::ExcInvalidObject());
+           Assert (local_source.size() == accessor.get_fe().dofs_per_cell,
+                   typename BaseClass::ExcVectorDoesNotMatch());
+           Assert (accessor.dof_handler->n_dofs() == global_destination.size(),
+                   typename BaseClass::ExcVectorDoesNotMatch());
+
+                                            // check as in documentation that
+                                            // cell is either active, or dofs
+                                            // are only in vertices
+           Assert (!accessor.has_children()
+                   ||
+                   (accessor.get_fe().dofs_per_cell ==
+                    accessor.get_fe().dofs_per_vertex * GeometryInfo<dim>::vertices_per_cell),
+                   ExcMessage ("Cell must either be active, or all DoFs must be in vertices"));
+
+           const unsigned int n_dofs = local_source.size();
+
+           unsigned int * dofs = &accessor.dof_handler->levels[accessor.level()]
+                                 ->cell_dof_indices_cache[accessor.present_index * n_dofs];
+
+                                  // distribute cell vector
+           global_destination.add(n_dofs, dofs, 
+                                  &(*const_cast<dealii::Vector<number>*>(&local_source))(0));
+         }
+
+
+
+        template <int dim, int spacedim, typename number, class OutputVector>
+       static
+       void
+       distribute_local_to_global (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim> > &accessor,
+                                   const dealii::Vector<number> &local_source,
+                                   OutputVector                 &global_destination)
+          {
+           typedef
+             dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
+             BaseClass;
+
+           Assert (accessor->dof_handler != 0,
+                   typename BaseClass::ExcInvalidObject());
+           Assert (&accessor.get_fe() != 0,
+                   typename BaseClass::ExcInvalidObject());
+           Assert (local_source.size() == accessor.get_fe().dofs_per_cell,
+                   typename BaseClass::ExcVectorDoesNotMatch());
+           Assert (accessor.dof_handler->n_dofs() == global_destination.size(),
+                   typename BaseClass::ExcVectorDoesNotMatch());
+
+           const unsigned int n_dofs = local_source.size();
+
+//TODO[WB/MK]: This function could me made more efficient because it allocates memory, which could be avoided by passing in another argument as a scratch array. This should be fixed eventually
+  
+                                  // get indices of dofs
+           std::vector<unsigned int> dofs (n_dofs);
+           accessor.get_dof_indices (dofs);
+  
+                                  // distribute cell vector
+           global_destination.add (dofs, local_source);
+         }
+
+
+
+        template <int dim, int spacedim, typename number, class OutputMatrix>
+       static
+       void
+       distribute_local_to_global (const DoFCellAccessor<dealii::DoFHandler<dim,spacedim> > &accessor,
+                                   const dealii::FullMatrix<number> &local_source,
+                                   OutputMatrix                     &global_destination)
+          {
+           typedef
+             dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
+             BaseClass;
+
+           Assert (accessor.dof_handler != 0,
+                   typename BaseClass::ExcInvalidObject());
+           Assert (&accessor.get_fe() != 0,
+                   typename BaseClass::ExcInvalidObject());
+           Assert (local_source.m() == accessor.get_fe().dofs_per_cell,
+                   typename BaseClass::ExcMatrixDoesNotMatch());
+           Assert (local_source.n() == accessor.get_fe().dofs_per_cell,
+                   typename BaseClass::ExcMatrixDoesNotMatch());
+           Assert (accessor.dof_handler->n_dofs() == global_destination.m(),
+                   typename BaseClass::ExcMatrixDoesNotMatch());
+           Assert (accessor.dof_handler->n_dofs() == global_destination.n(),
+                   typename BaseClass::ExcMatrixDoesNotMatch());
+
+                                            // check as in documentation that
+                                            // cell is either active, or dofs
+                                            // are only in vertices
+           Assert (!accessor.has_children()
+                   ||
+                   (accessor.get_fe().dofs_per_cell ==
+                    accessor.get_fe().dofs_per_vertex * GeometryInfo<dim>::vertices_per_cell),
+                   ExcMessage ("Cell must either be active, or all DoFs must be in vertices"));
+
+           const unsigned int n_dofs = local_source.m();
+
+           unsigned int * dofs = &accessor.dof_handler->levels[accessor.level()]
+                                 ->cell_dof_indices_cache[accessor.present_index * n_dofs];
+
+                                  // distribute cell matrices
+           for (unsigned int i=0; i<n_dofs; ++i)
+             global_destination.add(dofs[i], n_dofs, dofs, 
+                                    &local_source(i,0));
+         }
+
+
+
+        template <int dim, int spacedim, typename number, class OutputMatrix>
+       static
+       void
+       distribute_local_to_global (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim> > &accessor,
+                                   const dealii::FullMatrix<number> &local_source,
+                                   OutputMatrix                     &global_destination)
+          {
+           typedef
+             dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
+             BaseClass;
+
+           Assert (accessor->dof_handler != 0,
+                   typename BaseClass::ExcInvalidObject());
+           Assert (&accessor.get_fe() != 0,
+                   typename BaseClass::ExcInvalidObject());
+           Assert (local_source.m() == accessor.get_fe().dofs_per_cell,
+                   typename BaseClass::ExcMatrixDoesNotMatch());
+           Assert (local_source.n() == accessor.get_fe().dofs_per_cell,
+                   typename BaseClass::ExcVectorDoesNotMatch());
+           Assert (accessor.dof_handler->n_dofs() == global_destination.m(),
+                   typename BaseClass::ExcMatrixDoesNotMatch());
+           Assert (accessor.dof_handler->n_dofs() == global_destination.n(),
+                   typename BaseClass::ExcMatrixDoesNotMatch());
+
+           const unsigned int n_dofs = local_source.size();
+
+//TODO[WB/MK]: This function could me made more efficient because it allocates memory, which could be avoided by passing in another argument as a scratch array.
+  
+                                  // get indices of dofs
+           std::vector<unsigned int> dofs (n_dofs);
+           accessor.get_dof_indices (dofs);
+  
+                                  // distribute cell vector
+           global_destination.add(dofs,local_source);
+         }
     };
   }
 }
@@ -2450,26 +2611,8 @@ DoFCellAccessor<DH>::
 distribute_local_to_global (const Vector<number> &local_source,
                            OutputVector         &global_destination) const
 {
-  Assert (this->dof_handler != 0,
-         typename BaseClass::ExcInvalidObject());
-  Assert (&this->get_fe() != 0,
-         typename BaseClass::ExcInvalidObject());
-  Assert (local_source.size() == this->get_fe().dofs_per_cell,
-         typename BaseClass::ExcVectorDoesNotMatch());
-  Assert (this->dof_handler->n_dofs() == global_destination.size(),
-         typename BaseClass::ExcVectorDoesNotMatch());
-
-  const unsigned int n_dofs = local_source.size();
-
-//TODO[WB]: This function could me made more efficient. First, it allocates memory, which could be avoided by passing in another argument as a scratch array (or by using the dof_indices cache for non-hp DoFHandlers, see the functions in the Implementation class above). second, the elementwise access is really slow if we use PETSc vectors/matrices. This should be fixed eventually
-  
-                                  // get indices of dofs
-  std::vector<unsigned int> dofs (n_dofs);
-  this->get_dof_indices (dofs);
-  
-                                  // distribute cell vector
-  for (unsigned int j=0; j<n_dofs; ++j)
-    global_destination(dofs[j]) += local_source(j);
+  internal::DoFCellAccessor::Implementation::
+    distribute_local_to_global (*this,local_source,global_destination);
 }
 
 
@@ -2482,31 +2625,8 @@ DoFCellAccessor<DH>::
 distribute_local_to_global (const FullMatrix<number> &local_source,
                            OutputMatrix             &global_destination) const
 {
-  Assert (this->dof_handler != 0,
-         typename BaseClass::ExcInvalidObject());
-  Assert (&this->get_fe() != 0,
-         typename BaseClass::ExcInvalidObject());
-  Assert (local_source.m() == this->get_fe().dofs_per_cell,
-         typename BaseClass::ExcVectorDoesNotMatch());
-  Assert (local_source.m() == local_source.n(),
-         typename BaseClass::ExcMatrixDoesNotMatch());
-  Assert (this->dof_handler->n_dofs() == global_destination.m(),
-         typename BaseClass::ExcMatrixDoesNotMatch());
-  Assert (global_destination.m() == global_destination.n(),
-         typename BaseClass::ExcMatrixDoesNotMatch());
-
-  const unsigned int n_dofs = local_source.m();
-
-//TODO[WB]: This function could me made more efficient. First, it allocates memory, which could be avoided by passing in another argument as a scratch array (or by using the dof_indices cache for non-hp DoFHandlers, see the functions in the Implementation class above). second, the elementwise access is really slow if we use PETSc vectors/matrices. This should be fixed eventually
-
-                                  // get indices of dofs
-  std::vector<unsigned int> dofs (n_dofs);
-  this->get_dof_indices (dofs);
-  
-                                  // distribute cell matrix
-  for (unsigned int i=0; i<n_dofs; ++i)
-    for (unsigned int j=0; j<n_dofs; ++j)
-      global_destination.add(dofs[i], dofs[j], local_source(i,j));
+  internal::DoFCellAccessor::Implementation::
+    distribute_local_to_global (*this,local_source,global_destination);
 }
 
 
index 7887ef7cc1d56fdb308c5f27a9174c7ee8638139..5fe707226a07d7fc7cfb51e5aa886cfe3983b537 100644 (file)
@@ -908,6 +908,44 @@ class BlockVectorBase : public Subscriptor
     BlockVectorBase &
     operator -= (const BlockVectorBase &V);
 
+
+                                      /**
+                                        * A collective add operation:
+                                        * This funnction adds a whole
+                                        * set of values stored in @p
+                                        * values to the vector
+                                        * components specified by @p
+                                        * indices.
+                                        */
+    template <typename Number>
+    void add (const std::vector<unsigned int> &indices,
+             const std::vector<Number>       &values);
+
+                                      /**
+                                       * This is a second collective
+                                       * add operation. As a
+                                       * difference, this function
+                                       * takes a deal.II vector of
+                                       * values.
+                                       */
+    template <typename Number>
+    void add (const std::vector<unsigned int> &indices,
+             const Vector<Number>            &values);
+
+                                     /**
+                                      * Take an address where
+                                      * <tt>n_elements</tt> are stored
+                                      * contiguously and add them into
+                                      * the vector. Handles all cases
+                                      * which are not covered by the
+                                      * other two <tt>add()</tt>
+                                      * functions above.
+                                      */
+    template <typename Number>
+    void add (const unsigned int  n_elements,
+             const unsigned int *indices,
+             const Number       *values);
+
                                     /**
                                      * $U(0-DIM)+=s$.  Addition of <tt>s</tt>
                                      * to all components. Note that
@@ -1865,13 +1903,56 @@ BlockVectorBase<VectorType>::operator -= (const BlockVectorBase<VectorType>& v)
     }
   return *this;
 }
+
+
+template <class VectorType>
+template <typename Number>
+inline
+void
+BlockVectorBase<VectorType> ::add (const std::vector<unsigned int> &indices,
+                                  const std::vector<Number>       &values)
+{
+  Assert (indices.size() == values.size(),
+         ExcDimensionMismatch(indices.size(), values.size()));
+  add (indices.size(), &indices[0], &values[0]);
+}
 
 
 
 template <class VectorType>
-void BlockVectorBase<VectorType>::add (const value_type a)
+template <typename Number>
+inline
+void
+BlockVectorBase<VectorType> ::add (const std::vector<unsigned int> &indices,
+                                  const Vector<Number>            &values)
+{
+  Assert (indices.size() == values.size(),
+         ExcDimensionMismatch(indices.size(), values.size()));
+  const unsigned int n_indices = indices.size();
+  for (unsigned int i=0; i<n_indices; ++i)
+    (*this)(indices[i]) += values(i);
+}
+
+
+
+template <class VectorType>
+template <typename Number>
+inline
+void
+BlockVectorBase<VectorType> ::add (const unsigned int  n_indices,
+                                  const unsigned int *indices,
+                                  const Number       *values)
 {
+  for (unsigned int i=0; i<n_indices; ++i)
+    (*this)(indices[i]) += values[i];
+}
+
+
 
+template <class VectorType>
+void BlockVectorBase<VectorType>::add (const value_type a)
+{
   Assert (numbers::is_finite(a), 
           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
index 3fd91ec336be8b18a2352525866819a6d3775f70..a6c2ca41c643688e976bb163d5238a602b842ab9 100644 (file)
@@ -744,68 +744,71 @@ distribute_local_to_global (const Vector<double>            &local_vector,
     }
 
   const unsigned int n_local_dofs = local_vector.size();
-  for (unsigned int i=0; i<n_local_dofs; ++i)
-    {
+  if (lines.size() == 0)
+    global_vector.add(local_dof_indices, local_vector);
+  else
+    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. then, we
                                   // simply copy over the data.
-      const std::vector<ConstraintLine>::const_iterator position =
-       find_constraint(local_dof_indices[i]);
+       const std::vector<ConstraintLine>::const_iterator position =
+         find_constraint(local_dof_indices[i]);
 
-      if (position==lines.end())
-       {
-         global_vector(local_dof_indices[i]) += local_vector(i);
-         continue;
-       }
+       if (position==lines.end())
+         {
+           global_vector(local_dof_indices[i]) += local_vector(i);
+           continue;
+         }
 
-      if (use_matrix)
-       {
-         const double val = position->inhomogeneity;
-         if (val != 0)
-           for (unsigned int j=0; j<n_local_dofs; ++j)
-             {
-               const std::vector<ConstraintLine>::const_iterator
-                 position_j = find_constraint(local_dof_indices[j]);
+       if (use_matrix)
+         {
+           const double val = position->inhomogeneity;
+           if (val != 0)
+             for (unsigned int j=0; j<n_local_dofs; ++j)
+               {
+                 const std::vector<ConstraintLine>::const_iterator
+                   position_j = find_constraint(local_dof_indices[j]);
 
-               if (position_j == lines.end())
-                 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;
-
-                   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[q].first)
-                         -= val * position_j->entries[q].second * matrix_entry;
-                     }
-                 }
-             }
-       }
-      else
+                 if (position_j == lines.end())
+                   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;
+
+                     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[q].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."));
+         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)
-       {
-         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;
-       }
-    }
+       for (unsigned int j=0; j<position->entries.size(); ++j)
+         {
+           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 5a818e8a18869238c9e717b9b68e193d3b332582..62fadd1b5791248d2c0c9c0a3bd08b7fa5c17b39 100644 (file)
@@ -468,7 +468,7 @@ namespace SparseMatrixIterators
  * @<double@></tt>; others can be generated in application programs (see the
  * section on @ref Instantiations in the manual).
  *
- * @author Essentially everyone who has ever worked on deal.II, 1994-2007
+ * @author Essentially everyone who has ever worked on deal.II, 1994-2009
  */
 template <typename number>
 class SparseMatrix : public virtual Subscriptor
index e9d78239a50583aa081e4834fb895a782baebc1c..ccec8767780993883c4d7c37b726fac349a022e8 100644 (file)
@@ -21,6 +21,7 @@
 #include <boost/lambda/lambda.hpp>
 
 #include <cstdio>
+#include <vector>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -653,6 +654,43 @@ class Vector : public Subscriptor
                                      */
     Vector<Number> & operator -= (const Vector<Number> &V);
 
+                                      /**
+                                        * A collective add operation:
+                                        * This funnction adds a whole
+                                        * set of values stored in @p
+                                        * values to the vector
+                                        * components specified by @p
+                                        * indices.
+                                        */
+    template <typename OtherNumber>
+    void add (const std::vector<unsigned int> &indices,
+             const std::vector<OtherNumber>  &values);
+
+                                      /**
+                                       * This is a second collective
+                                       * add operation. As a
+                                       * difference, this function
+                                       * takes a deal.II vector of
+                                       * values.
+                                       */
+    template <typename OtherNumber>
+    void add (const std::vector<unsigned int> &indices,
+             const Vector<OtherNumber>       &values);
+
+                                     /**
+                                      * Take an address where
+                                      * <tt>n_elements</tt> are stored
+                                      * contiguously and add them into
+                                      * the vector. Handles all cases
+                                      * which are not covered by the
+                                      * other two <tt>add()</tt>
+                                      * functions above.
+                                      */
+    template <typename OtherNumber>
+    void add (const unsigned int  n_elements,
+             const unsigned int *indices,
+             const OtherNumber  *values);
+
                                     /**
                                      * Addition of @p s to all
                                      * components. Note that @p s is a
@@ -1254,7 +1292,55 @@ Vector<Number>::scale (const Number factor)
                       (factor*boost::lambda::_1),
                       internal::Vector::minimum_parallel_grain_size);
 }
+
+
+template <typename Number>
+template <typename Number2>
+inline
+void
+Vector<Number>::add (const std::vector<unsigned int> &indices,
+                    const std::vector<Number2>      &values)
+{
+  Assert (indices.size() == values.size(),
+         ExcDimensionMismatch(indices.size(), values.size()));
+  add (indices.size(), &indices[0], &values[0]);
+}
+
+
+
+template <typename Number>
+template <typename Number2>
+inline
+void
+Vector<Number>::add (const std::vector<unsigned int> &indices,
+                    const Vector<Number2>           &values)
+{
+  Assert (indices.size() == values.size(),
+         ExcDimensionMismatch(indices.size(), values.size()));
+  add (indices.size(), &indices[0], values.val);
+}
+
+
+
+template <typename Number>
+template <typename Number2>
+inline
+void
+Vector<Number>::add (const unsigned int  n_indices,
+                    const unsigned int *indices,
+                    const Number2      *values)
+{
+  for (unsigned int i=0; i<n_indices; ++i)
+    {
+      Assert (indices[i] < vec_size, ExcIndexRange(indices[i],0,vec_size));
+      Assert (numbers::is_finite(values[i]),
+             ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
+      val[indices[i]] += values[i];
+    }
+}
 
 
 template <typename Number>

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.