]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Change ConstraintMatrix implementation so that it has fast random access to Constrain...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 30 Sep 2009 15:01:07 +0000 (15:01 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 30 Sep 2009 15:01:07 +0000 (15:01 +0000)
git-svn-id: https://svn.dealii.org/trunk@19618 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_accessor.h
deal.II/deal.II/include/dofs/dof_accessor.templates.h
deal.II/doc/news/changes.h
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 f3151108eaffa460e52b717abd9ce915c1bf94db..36cacee7ed29264c3bd598834e5d46e3b6590c68 100644 (file)
@@ -26,6 +26,7 @@ DEAL_II_NAMESPACE_OPEN
 template <typename number> class FullMatrix;
 template <typename number> class SparseMatrix;
 template <typename number> class Vector;
+class ConstraintMatrix;
 
 template <typename Accessor> class TriaRawIterator;
 
@@ -923,6 +924,76 @@ class DoFCellAccessor :  public DoFAccessor<DH::dimension,DH>
     void get_dof_values (const InputVector &values,
                         Vector<number>    &local_values) const;
 
+                                    /**
+                                     * Return the values of the given vector
+                                     * restricted to the dofs of this
+                                     * cell in the standard ordering: dofs
+                                     * on vertex 0, dofs on vertex 1, etc,
+                                     * dofs on line 0, dofs on line 1, etc,
+                                     * dofs on quad 0, etc.
+                                     *
+                                     * The vector has to have the
+                                     * right size before being passed
+                                     * to this function. This
+                                     * function is only callable for
+                                     * active cells.
+                                     *
+                                     * The input vector may be either
+                                     * a <tt>Vector<float></tt>,
+                                     * Vector<double>, or a
+                                     * BlockVector<double>, or a
+                                     * PETSc or Trilinos vector if
+                                     * deal.II is compiled to support
+                                     * these libraries. It is in the
+                                     * responsibility of the caller
+                                     * to assure that the types of
+                                     * the numbers stored in input
+                                     * and output vectors are
+                                     * compatible and with similar
+                                     * accuracy.
+                                     */
+    template <class InputVector, typename ForwardIterator>
+    void get_dof_values (const InputVector &values,
+                        ForwardIterator    local_values_begin,
+                        ForwardIterator    local_values_end) const;
+
+                                    /**
+                                     * Return the values of the given vector
+                                     * restricted to the dofs of this
+                                     * cell in the standard ordering: dofs
+                                     * on vertex 0, dofs on vertex 1, etc,
+                                     * dofs on line 0, dofs on line 1, etc,
+                                     * dofs on quad 0, etc.
+                                     *
+                                     * The vector has to have the
+                                     * right size before being passed
+                                     * to this function. This
+                                     * function is only callable for
+                                     * active cells.
+                                     *
+                                     * The input vector may be either a
+                                     * <tt>Vector<float></tt>,
+                                     * Vector<double>, or a
+                                     * BlockVector<double>, or a PETSc or
+                                     * Trilinos vector if deal.II is
+                                     * compiled to support these
+                                     * libraries. It is in the
+                                     * responsibility of the caller to
+                                     * assure that the types of the numbers
+                                     * stored in input and output vectors
+                                     * are compatible and with similar
+                                     * accuracy. The ConstraintMatrix
+                                     * passed as an argument to this
+                                     * function makes sure that constraints
+                                     * are correctly distributed when the
+                                     * dof values are calculated.
+                                     */
+    template <class InputVector, typename ForwardIterator>
+    void get_dof_values (const ConstraintMatrix &constraints,
+                        const InputVector      &values,
+                        ForwardIterator         local_values_begin,
+                        ForwardIterator         local_values_end) const;
+
                                     /**
                                      * This function is the counterpart to
                                      * get_dof_values(): it takes a vector
@@ -1111,6 +1182,50 @@ class DoFCellAccessor :  public DoFAccessor<DH::dimension,DH>
     distribute_local_to_global (const Vector<number> &local_source,
                                 OutputVector         &global_destination) const;
 
+                                    /**
+                                     * Distribute a local (cell based)
+                                     * vector in iterator format to a
+                                     * global one by mapping the local
+                                     * numbering of the degrees of freedom
+                                     * to the global one and entering the
+                                     * local values into the global vector.
+                                     *
+                                     * The elements are <em>added</em> up
+                                     * to the elements in the global
+                                     * vector, rather than just set, since
+                                     * this is usually what one wants.
+                                     */
+    template <typename ForwardIterator, typename OutputVector>
+    void
+    distribute_local_to_global (ForwardIterator   local_source_begin,
+                               ForwardIterator   local_source_end,
+                               OutputVector     &global_destination) const;
+
+                                    /**
+                                     * Distribute a local (cell based)
+                                     * vector in iterator format to a
+                                     * global one by mapping the local
+                                     * numbering of the degrees of freedom
+                                     * to the global one and entering the
+                                     * local values into the global vector.
+                                     *
+                                     * The elements are <em>added</em> up
+                                     * to the elements in the global
+                                     * vector, rather than just set, since
+                                     * this is usually what one
+                                     * wants. Moreover, the
+                                     * ConstraintMatrix passed to this
+                                     * function makes sure that also
+                                     * constraints are eliminated in this
+                                     * process.
+                                     */
+    template <typename ForwardIterator, typename OutputVector>
+    void
+    distribute_local_to_global (const ConstraintMatrix &constraints,
+                                ForwardIterator         local_source_begin,
+                               ForwardIterator         local_source_end,
+                               OutputVector           &global_destination) const;
+
                                     /**
                                      * This function does much the
                                      * same as the
@@ -1127,6 +1242,19 @@ class DoFCellAccessor :  public DoFAccessor<DH::dimension,DH>
     distribute_local_to_global (const FullMatrix<number> &local_source,
                                 OutputMatrix             &global_destination) const;
     
+                                    /**
+                                     * This function does what the two
+                                     * <tt>distribute_local_to_global</tt>
+                                     * functions with vector and matrix
+                                     * argument do, but all at once.
+                                     */
+    template <typename number, typename OutputMatrix, typename OutputVector>
+    void
+    distribute_local_to_global (const FullMatrix<number> &local_matrix,
+                               const Vector<number>     &local_vector,
+                                OutputMatrix             &global_matrix,
+                               OutputVector             &global_vector) const;
+    
                                     /**
                                      * @}
                                      */
index 312a82e64d84d28a95888ee36d2dec5212ab70de..c95b59a760b7d98e95682beb545c54bdb2b5dbe5 100644 (file)
@@ -15,6 +15,7 @@
 
 
 #include <base/config.h>
+#include <lac/constraint_matrix.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_levels.h>
 #include <dofs/dof_faces.h>
@@ -48,9 +49,10 @@ DoFAccessor (const Triangulation<DH::dimension,DH::space_dimension> *tria,
             const int                 index,
             const DH                 *dof_handler)
                :
-               internal::DoFAccessor::Inheritance<structdim,DH::dimension,DH::space_dimension>::BaseClass (tria,
-                                                                                                           level,
-                                                                                                           index),
+               internal::DoFAccessor::Inheritance<structdim,DH::dimension,
+                                                   DH::space_dimension>::BaseClass (tria,
+                                                                                   level,
+                                                                                   index),
                 dof_handler(const_cast<DH*>(dof_handler))
 {}
 
@@ -191,10 +193,10 @@ namespace internal
        static
        unsigned int
        get_dof_index (const dealii::DoFHandler<1,spacedim>   &dof_handler,
-                      const unsigned int              obj_level,
-                      const unsigned int              obj_index,
-                      const unsigned int              fe_index,
-                      const unsigned int              local_index,
+                      const unsigned int                      obj_level,
+                      const unsigned int                      obj_index,
+                      const unsigned int                      fe_index,
+                      const unsigned int                      local_index,
                       internal::int2type<1>)
          {
            return dof_handler.levels[obj_level]->lines.
@@ -2072,17 +2074,16 @@ namespace internal
                                          * indices that we hold for each
                                          * cell.
                                          */
-       template <int dim, int spacedim, class InputVector, typename number>
+       template <int dim, int spacedim, class InputVector, typename ForwardIterator>
        static
        void
        get_dof_values (const DoFCellAccessor<DoFHandler<dim,spacedim> > &accessor,
                        const InputVector &values,
-                       dealii::Vector<number> &local_values)
+                       ForwardIterator    local_values_begin,
+                       ForwardIterator    local_values_end)
          {
-           typedef
-             dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
-             BaseClass;
-           Assert (local_values.size() == accessor.get_fe().dofs_per_cell,
+           typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > BaseClass;
+           Assert (local_values_end-local_values_begin == accessor.get_fe().dofs_per_cell,
                    typename BaseClass::ExcVectorDoesNotMatch());
            Assert (values.size() == accessor.get_dof_handler().n_dofs(),
                    typename BaseClass::ExcVectorDoesNotMatch());
@@ -2097,9 +2098,10 @@ namespace internal
                    ExcMessage ("Cell must either be active, or all DoFs must be in vertices"));
 
            unsigned int *cache = &accessor.dof_handler->levels[accessor.level()]
-                                 ->cell_dof_indices_cache[accessor.present_index * accessor.get_fe().dofs_per_cell];
-           for (unsigned int i=0; i<accessor.get_fe().dofs_per_cell; ++i, ++cache)
-             local_values(i) = values(*cache);
+                                 ->cell_dof_indices_cache[accessor.present_index * 
+                                                          accessor.get_fe().dofs_per_cell];
+           for ( ; local_values_begin != local_values_end; ++local_values_begin, ++cache)
+             *local_values_begin = values(*cache);
          }
 
                                         /**
@@ -2109,22 +2111,98 @@ namespace internal
                                          * not have a cache for the local
                                          * DoF indices.
                                          */
-       template <int dim, int spacedim, class InputVector, typename number>
+       template <int dim, int spacedim, class InputVector, typename ForwardIterator>
        static
        void
        get_dof_values (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim> > &accessor,
                        const InputVector &values,
-                       dealii::Vector<number>    &local_values)
+                       ForwardIterator    local_values_begin,
+                       ForwardIterator    local_values_end)
          {
+           typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > BaseClass;
+
                                             // no caching for hp::DoFHandler
                                             // implemented
+           Assert (local_values_end-local_values_begin == accessor.get_fe().dofs_per_cell,
+                   typename BaseClass::ExcVectorDoesNotMatch());
            const unsigned int dofs_per_cell = accessor.get_fe().dofs_per_cell;
 
            std::vector<unsigned int> local_dof_indices (dofs_per_cell);
            get_dof_indices (accessor, local_dof_indices);
 
-           for (unsigned int i=0; i<dofs_per_cell; ++i)
-             local_values(i) = values(local_dof_indices[i]);
+           for (unsigned int i=0; i<dofs_per_cell; ++i, ++local_values_begin)
+             *local_values_begin = values(local_dof_indices[i]);
+         }
+
+
+                                        /**
+                                         * A function that collects the
+                                         * values of degrees of freedom. This
+                                         * function works for ::DoFHandler
+                                         * and all template arguments and
+                                         * uses the data from the cache of
+                                         * indices that we hold for each
+                                         * cell.
+                                         */
+       template <int dim, int spacedim, class InputVector, typename ForwardIterator>
+       static
+       void
+       get_dof_values (const DoFCellAccessor<DoFHandler<dim,spacedim> > &accessor,
+                       const ConstraintMatrix &constraints,
+                       const InputVector      &values,
+                       ForwardIterator         local_values_begin,
+                       ForwardIterator         local_values_end)
+         {
+           typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > BaseClass;
+           Assert (local_values_end-local_values_begin == accessor.get_fe().dofs_per_cell,
+                   typename BaseClass::ExcVectorDoesNotMatch());
+           Assert (values.size() == accessor.get_dof_handler().n_dofs(),
+                   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"));
+
+           unsigned int *cache = &accessor.dof_handler->levels[accessor.level()]
+                                 ->cell_dof_indices_cache[accessor.present_index * 
+                                                          accessor.get_fe().dofs_per_cell];
+           constraints.get_dof_values(values, *cache, local_values_begin,
+                                      local_values_end);
+         }
+
+                                        /**
+                                         * Same function as above except
+                                         * that it works for
+                                         * hp::DoFHandler objects that do
+                                         * not have a cache for the local
+                                         * DoF indices.
+                                         */
+       template <int dim, int spacedim, class InputVector, typename ForwardIterator>
+       static
+       void
+       get_dof_values (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim> > &accessor,
+                       const ConstraintMatrix &constraints,
+                       const InputVector      &values,
+                       ForwardIterator         local_values_begin,
+                       ForwardIterator         local_values_end)
+         {
+                                            // no caching for hp::DoFHandler
+                                            // implemented
+           typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > BaseClass;
+           Assert (local_values_end-local_values_begin == accessor.get_fe().dofs_per_cell,
+                   typename BaseClass::ExcVectorDoesNotMatch());
+           const unsigned int dofs_per_cell = accessor.get_fe().dofs_per_cell;
+
+           std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+           get_dof_indices (accessor, local_dof_indices);
+
+           constraints.get_dof_values (values, local_dof_indices.begin(),
+                                       local_values_begin, local_values_end);
          }
 
 
@@ -2140,9 +2218,7 @@ namespace internal
                        const dealii::Vector<number> &local_values,
                        OutputVector &values)
          {
-           typedef
-             dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
-             BaseClass;
+           typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > BaseClass;
            Assert (local_values.size() == accessor.get_fe().dofs_per_cell,
                    typename BaseClass::ExcVectorDoesNotMatch());
            Assert (values.size() == accessor.get_dof_handler().n_dofs(),
@@ -2234,9 +2310,7 @@ namespace internal
                                             // ::DoFHandler only supports a
                                             // single active fe with index
                                             // zero
-           typedef
-             dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
-             BaseClass;
+           typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > BaseClass;
            Assert (i == 0, typename BaseClass::ExcInvalidObject());
          }
 
@@ -2248,9 +2322,7 @@ namespace internal
        set_active_fe_index (DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim> > &accessor,
                             const unsigned int                                      i)
          {
-           typedef
-             dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
-             BaseClass;
+           typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > BaseClass;
            Assert (accessor.dof_handler != 0,
                    typename BaseClass::ExcInvalidObject());
            Assert (static_cast<unsigned int>(accessor.level()) <
@@ -2266,22 +2338,20 @@ namespace internal
 
 
 
-        template <int dim, int spacedim, typename number, class OutputVector>
+        template <int dim, int spacedim, typename ForwardIterator, class OutputVector>
        static
        void
        distribute_local_to_global (const DoFCellAccessor<dealii::DoFHandler<dim,spacedim> > &accessor,
-                                   const dealii::Vector<number> &local_source,
-                                   OutputVector                 &global_destination)
+                                   ForwardIterator local_source_begin,
+                                   ForwardIterator local_source_end,
+                                   OutputVector   &global_destination)
           {
-           typedef
-             dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
-             BaseClass;
-
+           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,
+           Assert (local_source_end-local_source_begin == accessor.get_fe().dofs_per_cell,
                    typename BaseClass::ExcVectorDoesNotMatch());
            Assert (accessor.dof_handler->n_dofs() == global_destination.size(),
                    typename BaseClass::ExcVectorDoesNotMatch());
@@ -2295,38 +2365,109 @@ namespace internal
                     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();
+           const unsigned int n_dofs = local_source_end - local_source_begin;
 
            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,local_source.begin());
+           global_destination.add(n_dofs, dofs, local_source_begin);
          }
 
 
 
-        template <int dim, int spacedim, typename number, class OutputVector>
+        template <int dim, int spacedim, typename ForwardIterator, 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)
+                                   ForwardIterator local_source_begin,
+                                   ForwardIterator local_source_end,
+                                   OutputVector   &global_destination)
           {
-           typedef
-             dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
-             BaseClass;
+           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_end-local_source_begin == 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_end - local_source_begin;
+
+//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 (n_dofs, dofs.begin(), local_source_begin);
+         }
+
+
 
-           Assert (accessor->dof_handler != 0,
+        template <int dim, int spacedim, typename ForwardIterator, class OutputVector>
+       static
+       void
+       distribute_local_to_global (const DoFCellAccessor<dealii::DoFHandler<dim,spacedim> > &accessor,
+                                   const ConstraintMatrix &constraints,
+                                   ForwardIterator         local_source_begin,
+                                   ForwardIterator         local_source_end,
+                                   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,
+           Assert (local_source_end-local_source_begin == 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();
+                                            // 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_end - local_source_begin;
+
+           unsigned int * dofs = &accessor.dof_handler->levels[accessor.level()]
+                                 ->cell_dof_indices_cache[accessor.present_index * n_dofs];
+
+                                  // distribute cell vector
+           constraints.distribute_local_to_global (local_source_begin, local_source_end,
+                                                   dofs, global_destination);
+         }
+
+
+
+        template <int dim, int spacedim, typename ForwardIterator, class OutputVector>
+       static
+       void
+       distribute_local_to_global (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim> > &accessor,
+                                   const ConstraintMatrix &constraints,
+                                   ForwardIterator         local_source_begin,
+                                   ForwardIterator         local_source_end,
+                                   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_end-local_source_begin == 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_end - local_source_begin;
 
 //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
 
@@ -2335,7 +2476,8 @@ namespace internal
            accessor.get_dof_indices (dofs);
 
                                   // distribute cell vector
-           global_destination.add (dofs, local_source);
+           constraints.distribute_local_to_global (local_source_begin, local_source_end,
+                                                   dofs.begin(), global_destination);
          }
 
 
@@ -2347,10 +2489,7 @@ namespace internal
                                    const dealii::FullMatrix<number> &local_source,
                                    OutputMatrix                     &global_destination)
           {
-           typedef
-             dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
-             BaseClass;
-
+           typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > BaseClass;
            Assert (accessor.dof_handler != 0,
                    typename BaseClass::ExcInvalidObject());
            Assert (&accessor.get_fe() != 0,
@@ -2378,7 +2517,7 @@ namespace internal
            unsigned int * dofs = &accessor.dof_handler->levels[accessor.level()]
                                  ->cell_dof_indices_cache[accessor.present_index * n_dofs];
 
-                                  // distribute cell matrices
+                                  // distribute cell matrix
            for (unsigned int i=0; i<n_dofs; ++i)
              global_destination.add(dofs[i], n_dofs, dofs,
                                     &local_source(i,0));
@@ -2393,11 +2532,8 @@ namespace internal
                                    const dealii::FullMatrix<number> &local_source,
                                    OutputMatrix                     &global_destination)
           {
-           typedef
-             dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
-             BaseClass;
-
-           Assert (accessor->dof_handler != 0,
+           typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > BaseClass;
+           Assert (accessor.dof_handler != 0,
                    typename BaseClass::ExcInvalidObject());
            Assert (&accessor.get_fe() != 0,
                    typename BaseClass::ExcInvalidObject());
@@ -2418,9 +2554,103 @@ namespace internal
            std::vector<unsigned int> dofs (n_dofs);
            accessor.get_dof_indices (dofs);
 
-                                  // distribute cell vector
+                                  // distribute cell matrix
            global_destination.add(dofs,local_source);
          }
+
+
+
+        template <int dim, int spacedim, typename number, 
+                 class OutputMatrix, typename OutputVector>
+        static
+       void
+       distribute_local_to_global (const DoFCellAccessor<dealii::DoFHandler<dim,spacedim> > &accessor,
+                                   const dealii::FullMatrix<number> &local_matrix,
+                                   const dealii::Vector<number>     &local_vector,
+                                   OutputMatrix                     &global_matrix,
+                                   OutputVector                     &global_vector)
+          {
+           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_matrix.m() == accessor.get_fe().dofs_per_cell,
+                   typename BaseClass::ExcMatrixDoesNotMatch());
+           Assert (local_matrix.n() == accessor.get_fe().dofs_per_cell,
+                   typename BaseClass::ExcVectorDoesNotMatch());
+           Assert (accessor.dof_handler->n_dofs() == global_matrix.m(),
+                   typename BaseClass::ExcMatrixDoesNotMatch());
+           Assert (accessor.dof_handler->n_dofs() == global_matrix.n(),
+                   typename BaseClass::ExcMatrixDoesNotMatch());
+           Assert (local_vector.size() == accessor.get_fe().dofs_per_cell,
+                   typename BaseClass::ExcVectorDoesNotMatch());
+           Assert (accessor.dof_handler->n_dofs() == global_vector.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 = accessor.get_fe().dofs_per_cell;
+           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_matrix.add(dofs[i], n_dofs, dofs, &local_matrix(i,0));
+               global_vector(dofs[i]) += local_vector(i);
+             }
+         }
+
+
+
+        template <int dim, int spacedim, typename number, 
+                 class OutputMatrix, typename OutputVector>
+       static
+       void
+       distribute_local_to_global (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim> > &accessor,
+                                   const dealii::FullMatrix<number> &local_matrix,
+                                   const dealii::Vector<number>     &local_vector,
+                                   OutputMatrix                     &global_matrix,
+                                   OutputVector                     &global_vector)
+          {
+           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_matrix.m() == accessor.get_fe().dofs_per_cell,
+                   typename BaseClass::ExcMatrixDoesNotMatch());
+           Assert (local_matrix.n() == accessor.get_fe().dofs_per_cell,
+                   typename BaseClass::ExcVectorDoesNotMatch());
+           Assert (accessor.dof_handler->n_dofs() == global_matrix.m(),
+                   typename BaseClass::ExcMatrixDoesNotMatch());
+           Assert (accessor.dof_handler->n_dofs() == global_matrix.n(),
+                   typename BaseClass::ExcMatrixDoesNotMatch());
+           Assert (local_vector.size() == accessor.get_fe().dofs_per_cell,
+                   typename BaseClass::ExcVectorDoesNotMatch());
+           Assert (accessor.dof_handler->n_dofs() == global_vector.size(),
+                   typename BaseClass::ExcVectorDoesNotMatch());
+
+           const unsigned int n_dofs = local_matrix.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 matrix and vector
+           global_matrix.add(dofs,local_matrix);
+           global_vector.add(dofs,local_vector);
+         }
     };
   }
 }
@@ -2554,7 +2784,35 @@ DoFCellAccessor<DH>::get_dof_values (const InputVector &values,
                                     Vector<number>    &local_values) const
 {
   internal::DoFCellAccessor::Implementation
-    ::get_dof_values (*this, values, local_values);
+    ::get_dof_values (*this, values, local_values.begin(), local_values.end());
+}
+
+
+
+template <class DH>
+template <class InputVector, typename ForwardIterator>
+void
+DoFCellAccessor<DH>::get_dof_values (const InputVector &values,
+                                    ForwardIterator    local_values_begin,
+                                    ForwardIterator    local_values_end) const
+{
+  internal::DoFCellAccessor::Implementation
+    ::get_dof_values (*this, values, local_values_begin, local_values_end);
+}
+
+
+
+template <class DH>
+template <class InputVector, typename ForwardIterator>
+void
+DoFCellAccessor<DH>::get_dof_values (const ConstraintMatrix &constraints,
+                                    const InputVector      &values,
+                                    ForwardIterator         local_values_begin,
+                                    ForwardIterator         local_values_end) const
+{
+  internal::DoFCellAccessor::Implementation
+    ::get_dof_values (*this, constraints, values, 
+                     local_values_begin, local_values_end);
 }
 
 
@@ -2611,7 +2869,41 @@ distribute_local_to_global (const Vector<number> &local_source,
                            OutputVector         &global_destination) const
 {
   internal::DoFCellAccessor::Implementation::
-    distribute_local_to_global (*this,local_source,global_destination);
+    distribute_local_to_global (*this, local_source.begin(),
+                               local_source.end(), global_destination);
+}
+
+
+
+template <class DH>
+template <typename ForwardIterator, typename OutputVector>
+inline
+void
+DoFCellAccessor<DH>::
+distribute_local_to_global (ForwardIterator  local_source_begin,
+                           ForwardIterator  local_source_end,
+                           OutputVector    &global_destination) const
+{
+  internal::DoFCellAccessor::Implementation::
+    distribute_local_to_global (*this, local_source_begin,
+                               local_source_end, global_destination);
+}
+
+
+
+template <class DH>
+template <typename ForwardIterator, typename OutputVector>
+inline
+void
+DoFCellAccessor<DH>::
+distribute_local_to_global (const ConstraintMatrix &constraints,
+                           ForwardIterator         local_source_begin,
+                           ForwardIterator         local_source_end,
+                           OutputVector           &global_destination) const
+{
+  internal::DoFCellAccessor::Implementation::
+    distribute_local_to_global (*this, constraints, local_source_begin,
+                               local_source_end, global_destination);
 }
 
 
@@ -2629,6 +2921,23 @@ distribute_local_to_global (const FullMatrix<number> &local_source,
 }
 
 
+
+template <class DH>
+template <typename number, typename OutputMatrix, typename OutputVector>
+inline
+void
+DoFCellAccessor<DH>::
+distribute_local_to_global (const FullMatrix<number> &local_matrix,
+                           const Vector<number>     &local_vector,
+                           OutputMatrix             &global_matrix,
+                           OutputVector             &global_vector) const
+{
+  internal::DoFCellAccessor::Implementation::
+    distribute_local_to_global (*this,local_matrix,local_vector,
+                               global_matrix,global_vector);
+}
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
index d2e07c4bea9dc70593187f3023c1c514f61fecdf..436f7fea2efc4f29260abc9a684c0862040c5c37 100644 (file)
@@ -341,6 +341,18 @@ inconvenience this causes.
 
 <ol>
 
+  <li>
+    <p>
+    Improved: The ConstraintMatrix class now uses a cache for random access to
+    the constraint lines. This considerably increases performance of the 
+    *_local_to_global functions, where such an access pattern is usual. Moreover,
+    the ConstraintMatrix class has now a function get_dof_values that can import 
+    data from a global vector to a cell vector with respecting the constraints.
+    </p>
+  <br>
+  (Martin Kronbichler 2009/09/30)
+  </li>
+
   <li>
     <p>
     Fixed: SparsityTools::reorder_Cuthill_McKee would produce an error if the
index e590cf36d1a64536124649a634dda60f707cc971..393c2fc529eeadef3c3577c4c89849374ef259b6 100644 (file)
@@ -929,6 +929,21 @@ class ConstraintMatrix : public Subscriptor
     void condense (BlockSparseMatrix<number> &matrix,
                   BlockVectorType           &vector) const;
 
+                                    /**
+                                     * Delete hanging nodes in a
+                                     * vector. Sets all hanging node
+                                     * values to zero. The @p
+                                     * VectorType may be a
+                                     * Vector<float>, Vector<double>,
+                                     * BlockVector<tt><...></tt>, a
+                                     * PETSc or Trilinos vector
+                                     * wrapper class, or any other
+                                     * type having the same
+                                     * interface.
+                                     */
+    template <class VectorType>
+    void set_zero (VectorType &vec) const;
+
                                     /**
                                      * @}
                                      */
@@ -976,52 +991,144 @@ class ConstraintMatrix : public Subscriptor
                                       * vectors and matrices are fully
                                       * assembled.
                                      *
-                                     * 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.
+                                     * Note that this function will apply all
+                                      * constraints as if they were
+                                      * homogeneous. For correctly setting
+                                      * inhomogeneous constraints, use the
+                                      * similar function with a matrix
+                                      * argument or the function with both
+                                      * matrix and vector arguments.
                                      *
-                                     * The optional argument
+                                     * Note: This function is not
+                                     * thread-safe, so you will need to make
+                                     * sure that only on process at a time
+                                     * calls this function.
+                                      */
+    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 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.
+                                      *
+                                      * 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
+                                      * global object, one saves the call to
+                                      * the condense function after the
+                                      * vectors and matrices are fully
+                                      * assembled.
+                                     *
+                                     * The fourth 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.
+                                     * 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.
                                      *
                                      * Note: This function is not
-                                     * thread-safe, so you will need to
-                                     * make sure that only on process at a
-                                     * time calls this function.
+                                     * thread-safe, so you will need to make
+                                     * sure that only on process at a time
+                                     * calls this function.
                                       */
     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 FullMatrix<double>        &local_matrix = FullMatrix<double>()) const;
+                               const FullMatrix<double>        &local_matrix) const;
+
+
+                                     /**
+                                      * This function takes a pointer to 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 the
+                                      * entries in @p local_dof_indices
+                                      * indicate reasonable global vector
+                                      * entries, this function is happy with
+                                      * whatever it is given.
+                                      *
+                                      * 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 global object, one saves the
+                                      * call to the condense function after
+                                      * the vectors and matrices are fully
+                                      * assembled. Note that this function
+                                      * completely ignores inhomogeneous
+                                      * constraints.
+                                     *
+                                     * Note: This function is not
+                                     * thread-safe, so you will need to
+                                     * make sure that only on process at a
+                                     * time calls this function.
+                                      */
+    template <typename ForwardIteratorVec, typename ForwardIteratorInd,
+              class VectorType>
+    void
+    distribute_local_to_global (ForwardIteratorVec local_vector_begin,
+                                ForwardIteratorVec local_vector_end,
+                               ForwardIteratorInd local_indices_begin,
+                                VectorType        &global_vector) const;
 
                                      /**
                                       * This function takes a matrix of
@@ -1201,23 +1308,49 @@ class ConstraintMatrix : public Subscriptor
     add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
                                 SparsityType                    &sparsity_pattern,
                                 const bool                       keep_constrained_entries = true,
-                                const Table<2,bool>             &dof_mask = default_empty_table) const;
-
-                                    /**
-                                     * Delete hanging nodes in a
-                                     * vector. Sets all hanging node
-                                     * values to zero. The @p
-                                     * VectorType may be a
-                                     * Vector<float>, Vector<double>,
-                                     * BlockVector<tt><...></tt>, a
-                                     * PETSc or Trilinos vector
-                                     * wrapper class, or any other
-                                     * type having the same
-                                     * interface.
-                                     */
-    template <class VectorType>
-    void set_zero (VectorType &vec) const;
+                                const Table<2,bool>             &dof_mask = Table<2,bool>()) const;
 
+                                     /**
+                                      * This function imports values from a
+                                      * global vector (@p global_vector) by
+                                      * applying the constraints to a vector
+                                      * of local values, expressed in
+                                      * iterator format.  In most cases, the
+                                      * local values will be identified by
+                                      * the local dof values on a
+                                      * cell. However, as long as the
+                                      * entries in @p local_dof_indices
+                                      * indicate reasonable global vector
+                                      * entries, this function is happy with
+                                      * whatever it is given.
+                                      *
+                                      * If one of the elements of @p
+                                      * local_dof_indices belongs to a
+                                      * constrained node, then rather than
+                                      * writing the corresponding element of
+                                      * @p global_vector into @p
+                                      * local_vector, the constraints are
+                                      * resolved as the respective
+                                      * distribute function does, i.e., the
+                                      * local entry is constructed from the
+                                      * global entries to which this
+                                      * particular degree of freedom is
+                                      * constrained.
+                                     *
+                                      * In contrast to the similar function
+                                      * get_dof_values in the DoFAccessor
+                                      * class, this function does not need
+                                      * the constrained values to be
+                                      * correctly set (i.e., distribute to
+                                      * be called).
+                                      */
+    template <typename ForwardIteratorVec, typename ForwardIteratorInd,
+              class VectorType>
+    void
+    get_dof_values (const VectorType  &global_vector,
+                   ForwardIteratorInd local_indices_begin,
+                   ForwardIteratorVec local_vector_begin,
+                   ForwardIteratorVec local_vector_end) const;
 
                                     /**
                                      * @}
@@ -1363,9 +1496,9 @@ class ConstraintMatrix : public Subscriptor
                                         /**
                                          * Value of the inhomogeneity.
                                          */
-        double inhomogeneity;
+       double inhomogeneity;
 
-                                        /**
+                                         /**
                                          * This operator is a bit weird and
                                          * unintuitive: it compares the
                                          * line numbers of two lines. We
@@ -1379,7 +1512,7 @@ class ConstraintMatrix : public Subscriptor
                                          */
        bool operator < (const ConstraintLine &) const;
 
-                                        /**
+                                         /**
                                          * This operator is likewise weird:
                                          * it checks whether the line
                                          * indices of the two operands are
@@ -1412,17 +1545,21 @@ class ConstraintMatrix : public Subscriptor
                                      * to store the lines. This, however,
                                      * would mean a much more fractioned
                                      * heap since it allocates many small
-                                     * objects, ans would additionally make
+                                     * objects, and would additionally make
                                      * usage of this matrix much slower.
                                      */
     std::vector<ConstraintLine> 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
+                                     * A list of pointers that contains the
+                                     * address of the ConstraintLine of a
+                                     * constrained degree of freedom, or NULL
+                                     * if the degree of freedom is not
+                                     * constrained. The NULL return value
+                                     * returns thus 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,
@@ -1430,39 +1567,45 @@ class ConstraintMatrix : public Subscriptor
                                      * 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 save
+                                     * this time since we find any constraint
+                                     * in O(1) time or get to know that it a
+                                     * certain degree of freedom is not
+                                     * constrained.
                                      *
                                      * 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.
+                                     * other contexts as well, e.g. when one
+                                     * needs random access to the constraints
+                                     * as in all the functions that apply
+                                     * constraints on the fly while add cell
+                                     * contributions into vectors and
+                                     * matrices.
                                      */
-    std::vector<bool> constraint_line_exists;
+    std::vector<const ConstraintLine*> lines_cache;
 
                                     /**
                                      * Store whether the arrays are sorted.
@@ -1479,13 +1622,6 @@ class ConstraintMatrix : public Subscriptor
                                      */
     static bool check_zero_weight (const std::pair<unsigned int, double> &p);
 
-                                    /**
-                                     * Dummy table that serves as default
-                                     * argument for function
-                                     * <tt>add_entries_local_to_global()</tt>.
-                                     */
-    static const Table<2,bool> default_empty_table;
-
                                     /**
                                      * This function actually implements
                                      * the local_to_global function for
@@ -1540,26 +1676,6 @@ class ConstraintMatrix : public Subscriptor
                                 const Table<2,bool>             &dof_mask,
                                 internal::bool2type<true>) const;
 
-                                   /**
-                                    * This function returns a pointer to
-                                    * the respective ConstraintLine object
-                                    * if a dof is constrained, and a zero
-                                    * pointer otherwise.
-                                    */
-    std::vector<ConstraintLine>::const_iterator
-    find_constraint (const unsigned int line) const;
-
-                                   /**
-                                    * This function returns a pointer to
-                                    * the respective ConstraintLine object
-                                    * if a dof is constrained, and a zero
-                                    * pointer otherwise. This function is
-                                    * used if the constraint matrix has not
-                                    * been closed yet.
-                                    */
-    std::vector<ConstraintLine>::iterator
-    find_constraint (const unsigned int line);
-
 #ifdef DEAL_II_USE_TRILINOS
 //TODO: Make use of the following member thread safe
                                      /**
@@ -1588,7 +1704,7 @@ ConstraintMatrix::ConstraintMatrix (const ConstraintMatrix &constraint_matrix)
                :
                Subscriptor (),
                lines (constraint_matrix.lines),
-               constraint_line_exists (constraint_matrix.constraint_line_exists),
+               lines_cache (constraint_matrix.lines_cache),
                sorted (constraint_matrix.sorted)
 #ifdef DEAL_II_USE_TRILINOS
                ,vec_distribute ()
@@ -1606,22 +1722,36 @@ ConstraintMatrix::add_line (const unsigned int line)
 
                                   // check whether line already exists; it
                                   // may, in which case we can just quit
-  if ((line < constraint_line_exists.size())
+  if ((line < lines_cache.size())
       &&
-      (constraint_line_exists[line] == true))
-    return;
+      (lines_cache[line] != 0))
+    {
+      Assert (lines_cache[line]->line == line, ExcInternalError());
+      return;
+    }
 
                                   // if necessary enlarge vector of
-                                  // existing entries
-  if (line >= constraint_line_exists.size())
-    constraint_line_exists.resize (line+1, false);
-  constraint_line_exists[line] = true;
+                                  // existing entries for cache
+  if (line >= lines_cache.size())
+    lines_cache.resize (line+1,0);
+                                  // enlarge lines vector if necessary. 
+                                  // need to reset the pointers to the 
+                                  // ConstraintLine entries in that case, 
+                                  // since new memory will be allocated
+  if (lines.size() == lines.capacity())
+    {
+      lines.reserve (2*lines.size()+8);
+      std::vector<ConstraintLine>::const_iterator it = lines.begin();
+      for ( ; it != lines.end(); ++it)
+       lines_cache[it->line] = &*it;
+    }
 
                                   // push a new line to the end of the
                                   // list
   lines.push_back (ConstraintLine());
   lines.back().line = line;
   lines.back().inhomogeneity = 0.;
+  lines_cache[line] = &lines.back();
 }
 
 
@@ -1636,8 +1766,6 @@ ConstraintMatrix::add_entry (const unsigned int line,
   Assert (line != column,
          ExcMessage ("Can't constrain a degree of freedom to itself"));
 
-  std::vector<ConstraintLine>::iterator line_ptr = find_constraint(line);
-
                                   // if in debug mode, check whether an
                                   // entry for this column already
                                   // exists and if its the same as
@@ -1646,6 +1774,9 @@ ConstraintMatrix::add_entry (const unsigned int line,
                                   // in any case: exit the function if an
                                   // entry for this column already exists,
                                   // since we don't want to enter it twice
+  ConstraintLine* line_ptr = const_cast<ConstraintLine*>(lines_cache[line]);
+  Assert (line_ptr != 0, ExcInternalError());
+  Assert (line_ptr->line == line, ExcInternalError());
   for (std::vector<std::pair<unsigned int,double> >::const_iterator
          p=line_ptr->entries.begin();
        p != line_ptr->entries.end(); ++p)
@@ -1666,7 +1797,8 @@ void
 ConstraintMatrix::set_inhomogeneity (const unsigned int line,
                                     const double       value)
 {
-  find_constraint(line)->inhomogeneity = value;
+  ConstraintLine* line_ptr = const_cast<ConstraintLine*>(lines_cache[line]);
+  line_ptr->inhomogeneity = value;
 }
 
 
@@ -1675,9 +1807,9 @@ inline
 bool
 ConstraintMatrix::is_constrained (const unsigned int index) const
 {
-  return ((index < constraint_line_exists.size())
+  return ((index < lines_cache.size())
           &&
-          (constraint_line_exists[index] == true));
+          (lines_cache[index] != 0));
 }
 
 
@@ -1686,68 +1818,89 @@ inline
 bool
 ConstraintMatrix::is_inhomogeneously_constrained (const unsigned int index) const
 {
-  const std::vector<ConstraintLine>::const_iterator position
-    = find_constraint(index);
-  return position!=lines.end() ? position->inhomogeneity != 0 : false;
+  return (is_constrained(index) &&
+         lines_cache[index]->inhomogeneity != 0);
 }
 
 
 
+template <class VectorType>
 inline
-std::vector<ConstraintMatrix::ConstraintLine>::const_iterator
-ConstraintMatrix::find_constraint (const unsigned int line) const
+void
+ConstraintMatrix::
+distribute_local_to_global (const Vector<double>            &local_vector,
+                            const std::vector<unsigned int> &local_dof_indices,
+                            VectorType                      &global_vector) const
 {
-  if (is_constrained(line) == false)
-    return lines.end();
-
-  Assert (sorted==true, ExcMatrixNotClosed());
-  ConstraintLine index_comparison;
-  index_comparison.line = line;
+  Assert (local_vector.size() == local_dof_indices.size(),
+          ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
+  distribute_local_to_global (local_vector.begin(), local_vector.end(),
+                             local_dof_indices.begin(), global_vector);
+}
 
-  const std::vector<ConstraintLine>::const_iterator
-    position = std::lower_bound (lines.begin(),
-                                lines.end(),
-                                index_comparison);
-  Assert (position != lines.end(), ExcInternalError());
 
-                                  // check whether we've really found the
-                                  // right constraint.
-  Assert (position->line == line,
-         ExcInternalError());
 
-  return position;
+template <typename ForwardIteratorVec, typename ForwardIteratorInd,
+          class VectorType>
+inline
+void ConstraintMatrix::
+  distribute_local_to_global (ForwardIteratorVec local_vector_begin,
+                             ForwardIteratorVec local_vector_end,
+                             ForwardIteratorInd local_indices_begin,
+                             VectorType        &global_vector) const
+{
+  Assert (sorted == true, ExcMatrixNotClosed());
+  for ( ; local_vector_begin != local_vector_end;
+       ++local_vector_begin, ++local_indices_begin)
+    {
+      if (is_constrained(*local_indices_begin) == false)
+       global_vector(*local_indices_begin) += *local_vector_begin;
+      else
+       {
+         const ConstraintLine* position =
+           lines_cache[*local_indices_begin];
+         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_begin * position->entries[j].second;
+           }
+       }
+    }
 }
 
 
 
+template <typename ForwardIteratorVec, typename ForwardIteratorInd,
+          class VectorType>
 inline
-std::vector<ConstraintMatrix::ConstraintLine>::iterator
-ConstraintMatrix::find_constraint (const unsigned int line)
+void ConstraintMatrix::get_dof_values (const VectorType  &global_vector,
+                                      ForwardIteratorInd local_indices_begin,
+                                      ForwardIteratorVec local_vector_begin,
+                                      ForwardIteratorVec local_vector_end) const
 {
-  Assert (sorted==false, ExcMatrixIsClosed());
-  std::vector<ConstraintLine>::iterator line_ptr;
-  const std::vector<ConstraintLine>::const_iterator start=lines.begin();
-
-                                  // the usual case is that the line where
-                                  // a value is entered is 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)
-      break;
-
-                                  // if the loop didn't break, then
-                                  // line_ptr must be begin().
-                                  // we have an error if that doesn't
-                                  // point to 'line' then
-  Assert (line_ptr->line==line, ExcLineInexistant(line));
-                                  // second check: the respective
-                                  // flag in the
-                                  // constraint_line_exists field
-                                  // must exist
-  Assert (line < constraint_line_exists.size(), ExcInternalError());
-  Assert (constraint_line_exists[line] == true, ExcInternalError());
-
-  return line_ptr;
+  Assert (sorted == true, ExcMatrixNotClosed());
+  for ( ; local_vector_begin != local_vector_end;
+       ++local_vector_begin, ++local_indices_begin)
+    {
+      if (is_constrained(*local_indices_begin) == false)
+       *local_vector_begin = global_vector(*local_indices_begin);
+      else
+       {
+         const ConstraintLine* position =
+           lines_cache[*local_indices_begin];
+         typename VectorType::value_type value = position->inhomogeneity;
+         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."));
+             value += (global_vector(position->entries[j].first) *
+                       position->entries[j].second);
+           }
+         *local_vector_begin = value;
+       }
+    }
 }
 
 DEAL_II_NAMESPACE_CLOSE
index a6c2ca41c643688e976bb163d5238a602b842ab9..b5d535b48f7ad3dcc83a407de36a4cb3aee05222 100644 (file)
@@ -734,14 +734,10 @@ distribute_local_to_global (const Vector<double>            &local_vector,
   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()));
-    }
+  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();
   if (lines.size() == 0)
@@ -753,50 +749,41 @@ distribute_local_to_global (const Vector<double>            &local_vector,
                                   // 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]);
-
-       if (position==lines.end())
+       if (is_constrained(local_dof_indices[i]) == false)
          {
            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 ConstraintLine * position = 
+         lines_cache.size() <= local_dof_indices[i] ? 0 : 
+         lines_cache[local_dof_indices[i]];
+
+       const double val = position->inhomogeneity;
+       if (val != 0)
+         for (unsigned int j=0; j<n_local_dofs; ++j)
+           {
+             const ConstraintLine * position_j = 
+               lines_cache.size() <= local_dof_indices[j] ? 0 : 
+               lines_cache[local_dof_indices[j]];
+
+             if (position_j == 0)
+               global_vector(local_dof_indices[j]) -= val * local_matrix(j,i);
+             else
                {
-                 const std::vector<ConstraintLine>::const_iterator
-                   position_j = find_constraint(local_dof_indices[j]);
+                 const double matrix_entry = local_matrix(j,i);
+                 if (matrix_entry == 0)
+                   continue;
 
-                 if (position_j == lines.end())
-                   global_vector(local_dof_indices[j]) -= val * local_matrix(j,i);
-                 else
+                 for (unsigned int q=0; q<position_j->entries.size(); ++q)
                    {
-                     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;
-                       }
+                     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."));
+           }
 
                                   // now distribute the constraint,
                                   // but make sure we don't touch
@@ -1259,9 +1246,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                                   // resolved in a second step.
     for (unsigned int i = 0; i<n_local_dofs; ++i)
       {
-       const std::vector<ConstraintLine>::const_iterator position =
-         find_constraint(local_dof_indices[i]);
-       if (position == lines.end())
+       if (is_constrained(local_dof_indices[i]) == false)
          {
            my_indices[added_rows].global_row = local_dof_indices[i];
            my_indices[added_rows].local_row = i;
@@ -1270,7 +1255,9 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
          }
 
        constraint_lines.push_back (std::make_pair<unsigned int,
-                                   const ConstraintLine *>(i,&*position));
+                                   const ConstraintLine *>(i,lines_cache[local_dof_indices[i]]));
+       Assert (lines_cache[local_dof_indices[i]]->line == local_dof_indices[i],
+               ExcInternalError());
       }
     Assert (constraint_lines.size() + added_rows == n_local_dofs,
            ExcInternalError());
@@ -1577,9 +1564,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
     unsigned int added_rows = 0;
     for (unsigned int i = 0; i<n_local_dofs; ++i)
       {
-       const std::vector<ConstraintLine>::const_iterator position
-         = find_constraint (local_dof_indices[i]);
-       if (position == lines.end())
+       if (is_constrained(local_dof_indices[i]) == false)
          {
            my_indices[added_rows].global_row = local_dof_indices[i];
            my_indices[added_rows].local_row = i;
@@ -1588,7 +1573,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
          }
 
        constraint_lines.push_back (std::make_pair<unsigned int,
-                                   const ConstraintLine *>(i,&*position));
+                                   const ConstraintLine *>(i,lines_cache[local_dof_indices[i]]));
       }
     Assert (constraint_lines.size() + added_rows == n_local_dofs,
            ExcInternalError());
@@ -1827,9 +1812,7 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
       std::vector<std::pair<unsigned int, const ConstraintLine *> > constraint_lines;
       for (unsigned int i = 0; i<n_local_dofs; ++i)
        {
-         const std::vector<ConstraintLine>::const_iterator position
-           = find_constraint(local_dof_indices[i]);
-         if (position == lines.end())
+         if (is_constrained(local_dof_indices[i]) == false)
            {
              actual_dof_indices[added_rows] = local_dof_indices[i];
              ++added_rows;
@@ -1837,7 +1820,7 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
            }
 
          constraint_lines.push_back (std::make_pair<unsigned int,
-                                     const ConstraintLine *>(i,&*position));
+                                     const ConstraintLine *>(i,lines_cache[local_dof_indices[i]]));
       }
       Assert (constraint_lines.size() + added_rows == n_local_dofs,
              ExcInternalError());
@@ -1918,9 +1901,7 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
                                   // resolved in a second step.
     for (unsigned int i = 0; i<n_local_dofs; ++i)
       {
-       const std::vector<ConstraintLine>::const_iterator position
-         = find_constraint(local_dof_indices[i]);
-       if (position == lines.end())
+       if (is_constrained(local_dof_indices[i]) == false)
          {
            my_indices[added_rows].global_row = local_dof_indices[i];
            my_indices[added_rows].local_row = i;
@@ -1929,7 +1910,7 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
          }
 
        constraint_lines.push_back (std::make_pair<unsigned int,
-                                   const ConstraintLine *>(i,&*position));
+                                   const ConstraintLine *>(i,lines_cache[local_dof_indices[i]]));
       }
     Assert (constraint_lines.size() + added_rows == n_local_dofs,
            ExcInternalError());
@@ -2144,9 +2125,7 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
       std::vector<std::pair<unsigned int, const ConstraintLine *> > constraint_lines;
       for (unsigned int i = 0; i<n_local_dofs; ++i)
        {
-         const std::vector<ConstraintLine>::const_iterator position
-           = find_constraint(local_dof_indices[i]);
-         if (position == lines.end())
+         if (is_constrained(local_dof_indices[i]) == false)
            {
              actual_dof_indices[added_rows] = local_dof_indices[i];
              ++added_rows;
@@ -2154,7 +2133,7 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
            }
 
          constraint_lines.push_back (std::make_pair<unsigned int,
-                                     const ConstraintLine *>(i,&*position));
+                                     const ConstraintLine *>(i,lines_cache[local_dof_indices[i]]));
       }
       Assert (constraint_lines.size() + added_rows == n_local_dofs,
              ExcInternalError());
@@ -2254,8 +2233,7 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
                                   // resolved in a second step.
     for (unsigned int i = 0; i<n_local_dofs; ++i)
       {
-       if (constraint_line_exists.size() <= local_dof_indices[i] ||
-           constraint_line_exists[local_dof_indices[i]] == false)
+       if (is_constrained(local_dof_indices[i]) == false)
          {
            my_indices[added_rows].global_row = local_dof_indices[i];
            my_indices[added_rows].local_row = i;
@@ -2263,18 +2241,8 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
            continue;
          }
 
-       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);
-       Assert (position->line == local_dof_indices[i],
-               ExcInternalError());
-
        constraint_lines.push_back (std::make_pair<unsigned int,
-                                   const ConstraintLine *>(i,&*position));
+                                   const ConstraintLine *>(i,lines_cache[local_dof_indices[i]]));
       }
     Assert (constraint_lines.size() + added_rows == n_local_dofs,
            ExcInternalError());
index 6d806f72a5897d9308a9b94baeef55c9e739fcc0..e0bc0c15b69a19bea433f07f43d6195b79750775 100644 (file)
@@ -55,11 +55,6 @@ DEAL_II_NAMESPACE_OPEN
 
 
 
-                                       // Static member variable
-const Table<2,bool> ConstraintMatrix::default_empty_table = Table<2,bool>();
-
-
-
 bool
 ConstraintMatrix::check_zero_weight (const std::pair<unsigned int, double> &p)
 {
@@ -98,21 +93,15 @@ ConstraintMatrix::add_entries (const unsigned int                        line,
                                const std::vector<std::pair<unsigned int,double> > &col_val_pairs)
 {
   Assert (sorted==false, ExcMatrixIsClosed());
+  Assert (is_constrained(line), ExcLineInexistant(line));
 
-  std::vector<ConstraintLine>::iterator line_ptr;
-  const std::vector<ConstraintLine>::const_iterator start=lines.begin();
-                                  // the usual case is that the line where
-                                  // a value is entered is 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)
-      break;
+  ConstraintLine * line_ptr = const_cast<ConstraintLine*>(lines_cache[line]);
+  Assert (line_ptr->line == line, ExcInternalError());
 
                                   // if the loop didn't break, then
                                   // line_ptr must be begin().
                                   // we have an error if that doesn't
                                   // point to 'line' then
-  Assert (line_ptr->line==line, ExcLineInexistant(line));
 
                                   // if in debug mode, check whether an
                                   // entry for this column already
@@ -157,6 +146,24 @@ void ConstraintMatrix::close ()
                                   // sort the lines
   std::sort (lines.begin(), lines.end());
 
+                                  // update list of pointers and give the
+                                  // vector a sharp size since we won't
+                                  // modify the size any more after this
+                                  // point.
+  {
+    std::vector<const ConstraintLine*> new_lines (lines_cache.size());
+    for (std::vector<ConstraintLine>::const_iterator line=lines.begin();
+        line!=lines.end(); ++line)
+      new_lines[line->line] = &*line;
+    std::swap (lines_cache, new_lines);
+  }
+
+                                  // in debug mode: check whether we really
+                                  // set the pointers correctly.
+  for (unsigned int i=0; i<lines_cache.size(); ++i)
+    if (lines_cache[i] != 0)
+      Assert (i == lines_cache[i]->line, ExcInternalError());
+
                                   // first, strip zero entries, as we
                                   // have to do that only once
   for (std::vector<ConstraintLine>::iterator line = lines.begin();
@@ -229,19 +236,9 @@ void ConstraintMatrix::close ()
                Assert (dof_index != line->line,
                        ExcMessage ("Cycle in constraints detected!"));
 
-                                                // find the line
-                                                // corresponding to
-                                                // this entry. note
-                                                // that we have
-                                                // already sorted
-                                                // them to make this
-                                                // process faster
-               ConstraintLine test_line;
-               test_line.line = dof_index;
-               const std::vector<ConstraintLine>::const_iterator
-                 constrained_line = std::lower_bound (lines.begin(),
-                                                      lines.end(),
-                                                      test_line);
+               const ConstraintLine * constrained_line = lines_cache[dof_index];
+               Assert (constrained_line->line == dof_index, 
+                       ExcInternalError());
 
                                                 // now we have to
                                                 // replace an entry
@@ -468,15 +465,9 @@ void ConstraintMatrix::close ()
                                         // make sure that
                                         // entry->first is not the
                                         // index of a line itself
-       ConstraintLine test_line;
-       test_line.line = entry->first;
-       const std::vector<ConstraintLine>::const_iterator
-         test_line_position = std::lower_bound (lines.begin(),
-                                                lines.end(),
-                                                test_line);
-       Assert ((test_line_position == lines.end())
-               ||
-               (test_line_position->line != entry->first),
+       const ConstraintLine * it = entry->first < lines_cache.size() ? 
+         lines_cache[entry->first] : 0;
+       Assert (it == 0,
                ExcDoFConstrainedToConstrainedDoF(line->line, entry->first));
       };
 #endif
@@ -543,15 +534,8 @@ void ConstraintMatrix::merge (const ConstraintMatrix &other_constraints)
   const bool object_was_sorted = sorted;
   sorted = false;
 
-                                  // before we even start: merge the
-                                  // two flag arrays
-  if (other_constraints.constraint_line_exists.size() >
-      constraint_line_exists.size())
-    constraint_line_exists.resize (other_constraints.constraint_line_exists.size(),
-                                  false);
-  for (unsigned int i=0; i<other_constraints.constraint_line_exists.size(); ++i)
-    if (other_constraints.constraint_line_exists[i] == true)
-      constraint_line_exists[i] = true;
+  if (other_constraints.lines_cache.size() > lines_cache.size())
+    lines_cache.resize(other_constraints.lines_cache.size());
 
                                   // first action is to fold into the
                                   // present object possible
@@ -705,6 +689,9 @@ void ConstraintMatrix::merge (const ConstraintMatrix &other_constraints)
                                   // afterwards as well. otherwise
                                   // leave everything in the unsorted
                                   // state
+  for (std::vector<ConstraintLine>::const_iterator line=lines.begin();
+       line!=lines.end(); ++line)
+    lines_cache[line->line] = &*line;
   if (object_was_sorted == true)
     close ();
 }
@@ -713,8 +700,7 @@ void ConstraintMatrix::merge (const ConstraintMatrix &other_constraints)
 
 void ConstraintMatrix::shift (const unsigned int offset)
 {
-  constraint_line_exists.insert (constraint_line_exists.begin(), offset,
-                                false);
+  lines_cache.insert (lines_cache.begin(), offset, 0);
 
   for (std::vector<ConstraintLine>::iterator i = lines.begin();
        i != lines.end(); i++)
@@ -737,8 +723,8 @@ void ConstraintMatrix::clear ()
   }
 
   {
-    std::vector<bool> tmp;
-    constraint_line_exists.swap (tmp);
+    std::vector<const ConstraintLine*> tmp;
+    lines_cache.swap (tmp);
   }
 
 #ifdef DEAL_II_USE_TRILINOS
@@ -1871,7 +1857,6 @@ template<>
 void
 ConstraintMatrix::distribute (TrilinosWrappers::MPI::Vector &vec) const
 {
-  Assert (sorted == true, ExcMatrixNotClosed());
   ConstraintLine index_comparison;
   index_comparison.line = vec.local_range().first;
 
@@ -1955,39 +1940,14 @@ bool ConstraintMatrix::is_identity_constrained (const unsigned int index) const
   if (is_constrained(index) == false)
     return false;
 
-  if (sorted == true)
-    {
-      ConstraintLine index_comparison;
-      index_comparison.line = index;
+  const ConstraintLine * p = lines_cache[index];
+  Assert (p->line == index, ExcInternalError());
 
-      const std::vector<ConstraintLine>::const_iterator
-       p = std::lower_bound (lines.begin (),
-                             lines.end (),
-                             index_comparison);
                                       // return if an entry for this
                                       // line was found and if it has
                                       // only one entry equal to 1.0
-                                      //
-                                      // note that lower_bound only
-                                      // returns a valid iterator if
-                                      // 'index' is less than the
-                                      // largest line index in out
-                                      // constraints list
-      return ((p != lines.end()) &&
-             (p->line == index) &&
-             (p->entries.size() == 1) &&
-             (p->entries[0].second == 1.0));
-    }
-  else
-    {
-      for (std::vector<ConstraintLine>::const_iterator i=lines.begin();
-          i!=lines.end(); ++i)
-       if (i->line == index)
-         return ((i->entries.size() == 1) &&
-                 (i->entries[0].second == 1.0));
-
-      return false;
-    }
+  return ((p->entries.size() == 1) &&
+         (p->entries[0].second == 1.0));
 }
 
 
@@ -2087,7 +2047,7 @@ unsigned int
 ConstraintMatrix::memory_consumption () const
 {
   return (MemoryConsumption::memory_consumption (lines) +
-         MemoryConsumption::memory_consumption (constraint_line_exists) +
+         MemoryConsumption::memory_consumption (lines_cache) +
          MemoryConsumption::memory_consumption (sorted));
 }
 

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.