Implement a cache for the results of get_dof_indices

git-svn-id: https://svn.dealii.org/trunk@12778 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/deal.II/include/dofs/dof_accessor.h b/deal.II/deal.II/include/dofs/dof_accessor.h
index 9302b344e18adf57db9b88d4535e937826de764a..56915ac1b4b8733e852022f7c7b3f24e4be92daf 100644 (file)
--- a/deal.II/deal.II/include/dofs/dof_accessor.h
+++ b/deal.II/deal.II/include/dofs/dof_accessor.h
@@ -329,6 +329,22 @@ class DoFObjectAccessor : public DoFAccessor<DH>,
                                      * The vector has to have the
                                      * right size before being passed
                                      * to this function.
+                                     *
+                                     * This function is most often
+                                     * used on active objects (edges,
+                                     * faces, cells). It can be used
+                                     * on non-active objects as well
+                                     * (i.e. objects that have
+                                     * children), but only if the
+                                     * finite element under
+                                     * consideration has degrees of
+                                     * freedom exclusively on
+                                     * vertices. Otherwise, the
+                                     * function doesn't make much
+                                     * sense, since for example
+                                     * inactive edges do not have
+                                     * degrees of freedom associated
+                                     * with them at all.
                                      */
     void get_dof_indices (std::vector<unsigned int> &dof_indices) const;
 
@@ -649,6 +665,22 @@ class DoFObjectAccessor<1, DH> :
                                      *
                                      * It is assumed that the vector already
                                      * has the right size beforehand.
+                                     *
+                                     * This function is most often
+                                     * used on active objects (edges,
+                                     * faces, cells). It can be used
+                                     * on non-active objects as well
+                                     * (i.e. objects that have
+                                     * children), but only if the
+                                     * finite element under
+                                     * consideration has degrees of
+                                     * freedom exclusively on
+                                     * vertices. Otherwise, the
+                                     * function doesn't make much
+                                     * sense, since for example
+                                     * inactive edges do not have
+                                     * degrees of freedom associated
+                                     * with them at all.
                                      */
     void get_dof_indices (std::vector<unsigned int> &dof_indices) const;
 
@@ -896,6 +928,22 @@ class DoFObjectAccessor<2, DH> :
                                      *
                                      * It is assumed that the vector already
                                      * has the right size beforehand.
+                                     *
+                                     * This function is most often
+                                     * used on active objects (edges,
+                                     * faces, cells). It can be used
+                                     * on non-active objects as well
+                                     * (i.e. objects that have
+                                     * children), but only if the
+                                     * finite element under
+                                     * consideration has degrees of
+                                     * freedom exclusively on
+                                     * vertices. Otherwise, the
+                                     * function doesn't make much
+                                     * sense, since for example
+                                     * inactive edges do not have
+                                     * degrees of freedom associated
+                                     * with them at all.
                                      */
     void get_dof_indices (std::vector<unsigned int> &dof_indices) const;
 
@@ -1151,6 +1199,22 @@ class DoFObjectAccessor<3, DH> :
                                      *
                                      * It is assumed that the vector already
                                      * has the right size beforehand.
+                                     *
+                                     * This function is most often
+                                     * used on active objects (edges,
+                                     * faces, cells). It can be used
+                                     * on non-active objects as well
+                                     * (i.e. objects that have
+                                     * children), but only if the
+                                     * finite element under
+                                     * consideration has degrees of
+                                     * freedom exclusively on
+                                     * vertices. Otherwise, the
+                                     * function doesn't make much
+                                     * sense, since for example
+                                     * inactive edges do not have
+                                     * degrees of freedom associated
+                                     * with them at all.
                                      */
     void get_dof_indices (std::vector<unsigned int> &dof_indices) const;
 
@@ -1583,7 +1647,69 @@ class DoFCellAccessor :  public DoFObjectAccessor<DH::dimension,DH>
     template <class OutputVector, typename number>
     void set_dof_values_by_interpolation (const Vector<number> &local_values,
                                          OutputVector         &values) const;
+
+                                    /**
+                                     * Return the indices of the dofs of this
+                                     * quad 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.
+                                     *
+                                     * It is assumed that the vector already
+                                     * has the right size beforehand.
+                                     *
+                                     * This function reimplements the
+                                     * same function in the base
+                                     * class. The functions in the
+                                     * base classes are available for
+                                     * all geometric objects,
+                                     * i.e. even in 3d they can be
+                                     * used to access the dof indices
+                                     * of edges, for example. On the
+                                     * other hand, the most common
+                                     * case is clearly the use on
+                                     * cells, which is why we cache
+                                     * the array for each cell, but
+                                     * not edge. To retrieve the
+                                     * cached values, rather than
+                                     * collect the necessary
+                                     * information every time, this
+                                     * function overwrites the one in
+                                     * the base class.
+                                     *
+                                     * This function is most often
+                                     * used on active objects (edges,
+                                     * faces, cells). It can be used
+                                     * on non-active objects as well
+                                     * (i.e. objects that have
+                                     * children), but only if the
+                                     * finite element under
+                                     * consideration has degrees of
+                                     * freedom exclusively on
+                                     * vertices. Otherwise, the
+                                     * function doesn't make much
+                                     * sense, since for example
+                                     * inactive edges do not have
+                                     * degrees of freedom associated
+                                     * with them at all.
+                                     */
+    void get_dof_indices (std::vector<unsigned int> &dof_indices) const;
     
+  private:
+                                    /**
+                                     * Update the cache in which we
+                                     * store the dof indices of this
+                                     * cell.
+                                     */
+    void update_cell_dof_indices_cache () const;
+
+                                    /**
+                                     * Make the DoFHandler class a
+                                     * friend so that it can call the
+                                     * update_cell_dof_indices_cache()
+                                     * function
+                                     */
+    template <int> friend class DoFHandler;
 };
 
 
@@ -1677,6 +1803,28 @@ void
 DoFObjectAccessor<3,hp::DoFHandler<3> >::set_active_fe_index (const unsigned int i);
 
 
+template <>
+void
+DoFCellAccessor<DoFHandler<1> >::update_cell_dof_indices_cache () const;
+template <>
+void
+DoFCellAccessor<DoFHandler<2> >::update_cell_dof_indices_cache () const;
+template <>
+void
+DoFCellAccessor<DoFHandler<3> >::update_cell_dof_indices_cache () const;
+
+template <>
+void
+DoFCellAccessor<hp::DoFHandler<1> >::update_cell_dof_indices_cache () const;
+template <>
+void
+DoFCellAccessor<hp::DoFHandler<2> >::update_cell_dof_indices_cache () const;
+template <>
+void
+DoFCellAccessor<hp::DoFHandler<3> >::update_cell_dof_indices_cache () const;
+
+
+
 // include more templates
 #include "dof_accessor.templates.h"
 

diff --git a/deal.II/deal.II/include/dofs/dof_accessor.templates.h b/deal.II/deal.II/include/dofs/dof_accessor.templates.h
index b114d7833d4777e04a728cd5035772db3ca04e16..5c5d63f9095c7694f2db86344ff0b76cd55ee8f4 100644 (file)
--- a/deal.II/deal.II/include/dofs/dof_accessor.templates.h
+++ b/deal.II/deal.II/include/dofs/dof_accessor.templates.h
@@ -134,7 +134,7 @@ unsigned int
 DoFObjectAccessor<1,DH>::dof_index (const unsigned int i) const
 {
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
 
                                    // since the exception classes are
                                   // from a template dependent base
@@ -166,7 +166,7 @@ DoFObjectAccessor<1,DH>::vertex_dof_index (const unsigned int vertex,
                                            const unsigned int i) const
 {
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
 
                                    // since the exception classes are
                                   // from a template dependent base
@@ -198,7 +198,7 @@ void
 DoFObjectAccessor<1,DH>::get_dof_indices (std::vector<unsigned int> &dof_indices) const
 {
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
 
                                    // since the exception classes are
                                   // from a template dependent base
@@ -248,7 +248,7 @@ TriaIterator<DH::dimension,DoFObjectAccessor<1,DH> >
 DoFObjectAccessor<1,DH>::child (const unsigned int i) const
 {
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
 
   TriaIterator<dim,DoFObjectAccessor<1,DH> > q (this->tria,
                                                this->present_level+1,
@@ -451,7 +451,7 @@ DoFObjectAccessor<2,DH>::dof_index (const unsigned int i) const
   Assert (i<this->get_fe().dofs_per_quad,
          ExcIndexRange (i, 0, this->get_fe().dofs_per_quad));
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
 
   return this->dof_handler->levels[this->present_level]
     ->quad_dofs[this->present_index*this->get_fe().dofs_per_quad+i];
@@ -474,7 +474,7 @@ DoFObjectAccessor<2,DH>::vertex_dof_index (const unsigned int vertex,
   Assert (i<this->get_fe().dofs_per_vertex,
          ExcIndexRange (i, 0, this->get_fe().dofs_per_vertex));
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
 
   const unsigned int dof_number = (this->vertex_index(vertex) *
                                   this->get_fe().dofs_per_vertex +
@@ -483,7 +483,6 @@ DoFObjectAccessor<2,DH>::vertex_dof_index (const unsigned int vertex,
 }
 
 
-
 template <class DH>
 inline
 void
@@ -514,7 +513,7 @@ DoFObjectAccessor<2,DH>::get_dof_indices (std::vector<unsigned int> &dof_indices
          typename BaseClass::ExcNotActive());
   
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
   
   const unsigned int dofs_per_vertex = this->get_fe().dofs_per_vertex,
                     dofs_per_line   = this->get_fe().dofs_per_line,
@@ -538,7 +537,7 @@ DoFObjectAccessor<2,DH>::line (const unsigned int i) const
 {
   Assert (i<4, ExcIndexRange (i, 0, 4));
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
 
   return TriaIterator<dim,DoFObjectAccessor<1,DH> >
     (
@@ -557,7 +556,7 @@ TriaIterator<DH::dimension,DoFObjectAccessor<2,DH> >
 DoFObjectAccessor<2,DH>::child (const unsigned int i) const
 {
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
 
   TriaIterator<dim,DoFObjectAccessor<2,DH> > q (this->tria,
                                                this->present_level+1,
@@ -635,7 +634,7 @@ distribute_local_to_global (const Vector<number> &local_source,
          typename BaseClass::ExcVectorDoesNotMatch());
 
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
 
   const unsigned int n_dofs = local_source.size();
   
@@ -686,7 +685,7 @@ distribute_local_to_global (const FullMatrix<number> &local_source,
   Assert (global_destination.m() == global_destination.n(),
          typename BaseClass::ExcMatrixDoesNotMatch());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
 
   const unsigned int n_dofs = local_source.m();
 
@@ -772,7 +771,7 @@ DoFObjectAccessor<3,DH>::dof_index (const unsigned int i) const
   Assert (i<this->get_fe().dofs_per_hex,
          ExcIndexRange (i, 0, this->get_fe().dofs_per_hex));
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
 
   return this->dof_handler->levels[this->present_level]
     ->hex_dofs[this->present_index*this->get_fe().dofs_per_hex+i];
@@ -796,7 +795,7 @@ DoFObjectAccessor<3,DH>::vertex_dof_index (const unsigned int vertex,
   Assert (i<this->get_fe().dofs_per_vertex,
          ExcIndexRange (i, 0, this->get_fe().dofs_per_vertex));
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
 
   const unsigned int dof_number = (this->vertex_index(vertex) *
                                   this->get_fe().dofs_per_vertex +
@@ -822,7 +821,7 @@ DoFObjectAccessor<3,DH>::get_dof_indices (std::vector<unsigned int> &dof_indices
                                 this->get_fe().dofs_per_hex),
          typename BaseClass::ExcVectorDoesNotMatch());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
 
                                   // this function really only makes
                                   // sense on non-active objects if
@@ -837,7 +836,7 @@ DoFObjectAccessor<3,DH>::get_dof_indices (std::vector<unsigned int> &dof_indices
           8*this->get_fe().dofs_per_vertex),
          typename BaseClass::ExcNotActive());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
   
   const unsigned int dofs_per_vertex = this->get_fe().dofs_per_vertex,
                     dofs_per_line   = this->get_fe().dofs_per_line,
@@ -969,7 +968,7 @@ distribute_local_to_global (const Vector<number> &local_source,
   Assert (this->dof_handler->n_dofs() == global_destination.size(),
          typename BaseClass::ExcVectorDoesNotMatch());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
 
   const unsigned int n_dofs = local_source.size();
   
@@ -1021,7 +1020,7 @@ distribute_local_to_global (const FullMatrix<number> &local_source,
   Assert (global_destination.m() == global_destination.n(),
          typename BaseClass::ExcMatrixDoesNotMatch());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
 
   const unsigned int n_dofs = local_source.m();
 
@@ -1090,7 +1089,7 @@ DoFObjectAccessor<1,hp::DoFHandler<1> >::dof_index (const unsigned int i) const
   
   Assert (this->dof_handler != 0, BaseClass::ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
                                   // make sure a FE has been selected
                                   // and enough room was reserved
   Assert (&this->get_fe() != 0, BaseClass::ExcInvalidObject());
@@ -1127,7 +1126,7 @@ DoFObjectAccessor<1,hp::DoFHandler<2> >::dof_index (const unsigned int i) const
   
   Assert (this->dof_handler != 0, BaseClass::ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
                                   // make sure a FE has been selected
                                   // and enough room was reserved
   Assert (&this->get_fe() != 0, BaseClass::ExcInvalidObject());
@@ -1164,7 +1163,7 @@ DoFObjectAccessor<1,hp::DoFHandler<3> >::dof_index (const unsigned int i) const
   
   Assert (this->dof_handler != 0, BaseClass::ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
                                   // make sure a FE has been selected
                                   // and enough room was reserved
   Assert (&this->get_fe() != 0, BaseClass::ExcInvalidObject());
@@ -1195,7 +1194,7 @@ DoFObjectAccessor<1,hp::DoFHandler<1> >::set_dof_index (const unsigned int i,
   Assert (this->dof_handler != 0,
          BaseClass::ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
                                   // make sure a FE has been selected
                                   // and enough room was reserved
   Assert (&this->get_fe() != 0,
@@ -1226,7 +1225,7 @@ DoFObjectAccessor<1, hp::DoFHandler<2> >::set_dof_index (const unsigned int i,
   Assert (this->dof_handler != 0,
          BaseClass::ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
                                   // make sure a FE has been selected
                                   // and enough room was reserved
   Assert (&this->get_fe() != 0,
@@ -1261,7 +1260,7 @@ DoFObjectAccessor<1, hp::DoFHandler<3> >::set_dof_index (const unsigned int i,
   Assert (this->dof_handler != 0,
          BaseClass::ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
                                   // make sure a FE has been selected
                                   // and enough room was reserved
   Assert (&this->get_fe() != 0,
@@ -1319,7 +1318,7 @@ unsigned int
 DoFObjectAccessor<1,hp::DoFHandler<1> >::active_fe_index () const
 {
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
   Assert (static_cast<std::vector<unsigned int>::size_type>(this->present_index) <
          this->dof_handler->levels[this->present_level]->active_fe_indices.size (),
          ExcIndexRange (this->present_index, 0,
@@ -1336,7 +1335,7 @@ unsigned int
 DoFObjectAccessor<1,hp::DoFHandler<2> >::active_fe_index () const
 {
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
   Assert (static_cast<std::vector<unsigned int>::size_type>(this->present_index) <
          this->dof_handler->levels[this->present_level]->active_fe_indices.size (),
          ExcIndexRange (this->present_index, 0,
@@ -1353,7 +1352,7 @@ unsigned int
 DoFObjectAccessor<1,hp::DoFHandler<3> >::active_fe_index () const
 {
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
   Assert (static_cast<std::vector<unsigned int>::size_type>(this->present_index) <
          this->dof_handler->levels[this->present_level]->active_fe_indices.size (),
          ExcIndexRange (this->present_index, 0,
@@ -1372,7 +1371,7 @@ DoFObjectAccessor<1,hp::DoFHandler<1> >::set_active_fe_index (const unsigned int
   Assert (this->dof_handler != 0, ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) <
           this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
   Assert (static_cast<std::vector<unsigned int>::size_type>(this->present_index) <
          this->dof_handler->levels[this->present_level]->active_fe_indices.size (),
          ExcIndexRange (this->present_index, 0,
@@ -1392,7 +1391,7 @@ DoFObjectAccessor<1,hp::DoFHandler<2> >::set_active_fe_index (const unsigned int
 {
   Assert (this->dof_handler != 0, ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
   Assert (static_cast<std::vector<unsigned int>::size_type>(this->present_index) <
          this->dof_handler->levels[this->present_level]->active_fe_indices.size (),
          ExcIndexRange (this->present_index, 0,
@@ -1413,7 +1412,7 @@ DoFObjectAccessor<1,hp::DoFHandler<3> >::set_active_fe_index (const unsigned int
 {
   Assert (this->dof_handler != 0, ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
   Assert (static_cast<std::vector<unsigned int>::size_type>(this->present_index) <
          this->dof_handler->levels[this->present_level]->active_fe_indices.size (),
          ExcIndexRange (this->present_index, 0,
@@ -1437,7 +1436,7 @@ unsigned int DoFObjectAccessor<2,hp::DoFHandler<2> >::dof_index (const unsigned
   Assert (this->dof_handler != 0,
          BaseClass::ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
                                   // make sure a FE has been selected
                                   // and enough room was reserved
   Assert (&this->get_fe() != 0,
@@ -1467,7 +1466,7 @@ unsigned int DoFObjectAccessor<2,hp::DoFHandler<3> >::dof_index (const unsigned
   Assert (this->dof_handler != 0,
          BaseClass::ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
                                   // make sure a FE has been selected
                                   // and enough room was reserved
   Assert (&this->get_fe() != 0,
@@ -1499,7 +1498,7 @@ DoFObjectAccessor<2, hp::DoFHandler<2> >::set_dof_index (const unsigned int i,
   Assert (this->dof_handler != 0,
          BaseClass::ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
                                   // make sure a FE has been selected
                                   // and enough room was reserved
   Assert (&this->get_fe() != 0,
@@ -1531,7 +1530,7 @@ DoFObjectAccessor<2, hp::DoFHandler<3> >::set_dof_index (const unsigned int i,
   Assert (this->dof_handler != 0,
          BaseClass::ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
                                   // make sure a FE has been selected
                                   // and enough room was reserved
   Assert (&this->get_fe() != 0,
@@ -1578,7 +1577,7 @@ unsigned int
 DoFObjectAccessor<2,hp::DoFHandler<2> >::active_fe_index () const
 {
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
   Assert (static_cast<std::vector<unsigned int>::size_type>(this->present_index) <
          this->dof_handler->levels[this->present_level]->active_fe_indices.size (),
          ExcIndexRange (this->present_index, 0,
@@ -1595,7 +1594,7 @@ unsigned int
 DoFObjectAccessor<2,hp::DoFHandler<3> >::active_fe_index () const
 {
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
   Assert (static_cast<std::vector<unsigned int>::size_type>(this->present_index) <
          this->dof_handler->levels[this->present_level]->active_fe_indices.size (),
          ExcIndexRange (this->present_index, 0,
@@ -1612,9 +1611,9 @@ DoFObjectAccessor<2,hp::DoFHandler<2> >::set_active_fe_index (const unsigned int
 {
   Assert (this->dof_handler != 0, ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
   Assert (static_cast<std::vector<unsigned int>::size_type>(this->present_index) <
          this->dof_handler->levels[this->present_level]->active_fe_indices.size (),
          ExcIndexRange (this->present_index, 0,
@@ -1635,7 +1634,7 @@ DoFObjectAccessor<2,hp::DoFHandler<3> >::set_active_fe_index (const unsigned int
 {
   Assert (this->dof_handler != 0, ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
   Assert (static_cast<std::vector<unsigned int>::size_type>(this->present_index) <
          this->dof_handler->levels[this->present_level]->active_fe_indices.size (),
          ExcIndexRange (this->present_index, 0,
@@ -1660,7 +1659,7 @@ DoFObjectAccessor<3,hp::DoFHandler<3> >::dof_index (const unsigned int i) const
   Assert (this->dof_handler != 0,
          BaseClass::ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
                                   // make sure a FE has been selected
                                   // and enough room was reserved
   Assert (&this->get_fe() != 0,
@@ -1691,7 +1690,7 @@ DoFObjectAccessor<3, hp::DoFHandler<3> >::set_dof_index (const unsigned int i,
   Assert (this->dof_handler != 0,
          BaseClass::ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
                                   // make sure a FE has been selected
                                   // and enough room was reserved
   Assert (&this->get_fe() != 0,
@@ -1728,7 +1727,7 @@ unsigned int
 DoFObjectAccessor<3,hp::DoFHandler<3> >::active_fe_index () const
 {
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
   Assert (static_cast<std::vector<unsigned int>::size_type>(this->present_index) <
          this->dof_handler->levels[this->present_level]->active_fe_indices.size (),
          ExcIndexRange (this->present_index, 0,
@@ -1746,7 +1745,7 @@ DoFObjectAccessor<3,hp::DoFHandler<3> >::set_active_fe_index (const unsigned int
 {
   Assert (this->dof_handler != 0, ExcInvalidObject());
   Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
-          ExcMessage ("hp::DoFHandler not initialized"));
+          ExcMessage ("DoFHandler not initialized"));
   Assert (static_cast<std::vector<unsigned int>::size_type>(this->present_index) <
          this->dof_handler->levels[this->present_level]->active_fe_indices.size (),
          ExcIndexRange (this->present_index, 0,
@@ -1840,4 +1839,119 @@ DoFCellAccessor<hp::DoFHandler<3> >::face (const unsigned int i) const
 }
 
 
+
+
+template <>
+inline
+void
+DoFCellAccessor<DoFHandler<1> >::
+get_dof_indices (std::vector<unsigned int> &dof_indices) const
+{
+  Assert (dof_indices.size() == this->get_fe().dofs_per_cell,
+         ExcVectorDoesNotMatch());
+
+                                  // check as in documentation that
+                                  // cell is either active, or dofs
+                                  // are only in vertices
+  Assert (!this->has_children()
+         ||
+         (this->get_fe().dofs_per_cell ==
+          this->get_fe().dofs_per_vertex * GeometryInfo<1>::vertices_per_cell),
+         ExcMessage ("Cell must either be active, or all DoFs must be in vertices"));
+  
+  unsigned int *cache = &this->dof_handler->levels[this->present_level]
+                       ->cell_dof_indices_cache[this->present_index * this->get_fe().dofs_per_cell]; 
+  for (unsigned int i=0; i<this->get_fe().dofs_per_cell; ++i, ++cache)
+    dof_indices[i] = *cache;
+}
+
+
+
+template <>
+inline
+void
+DoFCellAccessor<hp::DoFHandler<1> >::
+get_dof_indices (std::vector<unsigned int> &dof_indices) const
+{
+                                  // no caching for hp::DoFHandler implemented 
+  DoFObjectAccessor<dim,hp::DoFHandler<dim> >::get_dof_indices (dof_indices);
+}
+
+
+
+template <>
+inline
+void
+DoFCellAccessor<DoFHandler<2> >::
+get_dof_indices (std::vector<unsigned int> &dof_indices) const
+{
+  Assert (dof_indices.size() == this->get_fe().dofs_per_cell,
+         ExcVectorDoesNotMatch());
+
+                                  // check as in documentation that
+                                  // cell is either active, or dofs
+                                  // are only in vertices
+  Assert (!this->has_children()
+         ||
+         (this->get_fe().dofs_per_cell ==
+          this->get_fe().dofs_per_vertex * GeometryInfo<2>::vertices_per_cell),
+         ExcMessage ("Cell must either be active, or all DoFs must be in vertices"));
+  
+  unsigned int *cache = &this->dof_handler->levels[this->present_level]
+                       ->cell_dof_indices_cache[this->present_index * this->get_fe().dofs_per_cell]; 
+  for (unsigned int i=0; i<this->get_fe().dofs_per_cell; ++i, ++cache)
+    dof_indices[i] = *cache;
+}
+
+
+
+template <>
+inline
+void
+DoFCellAccessor<hp::DoFHandler<2> >::
+get_dof_indices (std::vector<unsigned int> &dof_indices) const
+{
+                                  // no caching for hp::DoFHandler implemented 
+  DoFObjectAccessor<dim,hp::DoFHandler<dim> >::get_dof_indices (dof_indices);
+}
+
+
+
+template <>
+inline
+void
+DoFCellAccessor<DoFHandler<3> >::
+get_dof_indices (std::vector<unsigned int> &dof_indices) const
+{
+  Assert (dof_indices.size() == this->get_fe().dofs_per_cell,
+         ExcVectorDoesNotMatch());
+
+                                  // check as in documentation that
+                                  // cell is either active, or dofs
+                                  // are only in vertices
+  Assert (!this->has_children()
+         ||
+         (this->get_fe().dofs_per_cell ==
+          this->get_fe().dofs_per_vertex * GeometryInfo<3>::vertices_per_cell),
+         ExcMessage ("Cell must either be active, or all DoFs must be in vertices"));
+  
+  unsigned int *cache = &this->dof_handler->levels[this->present_level]
+                       ->cell_dof_indices_cache[this->present_index * this->get_fe().dofs_per_cell]; 
+  for (unsigned int i=0; i<this->get_fe().dofs_per_cell; ++i, ++cache)
+    dof_indices[i] = *cache;
+}
+
+
+
+template <>
+inline
+void
+DoFCellAccessor<hp::DoFHandler<3> >::
+get_dof_indices (std::vector<unsigned int> &dof_indices) const
+{
+                                  // no caching for hp::DoFHandler implemented 
+  DoFObjectAccessor<dim,hp::DoFHandler<dim> >::get_dof_indices (dof_indices);
+}
+
+
 #endif

diff --git a/deal.II/deal.II/include/dofs/dof_handler.h b/deal.II/deal.II/include/dofs/dof_handler.h
index 2a4af7d8b7f2ccf17cbab73f0086b9fd38407fea..9d5823292fcea0c2c18a5271c0e1cb4c6c2000f6 100644 (file)
--- a/deal.II/deal.II/include/dofs/dof_handler.h
+++ b/deal.II/deal.II/include/dofs/dof_handler.h
@@ -1055,6 +1055,11 @@ class DoFHandler  :  public Subscriptor
                                      * Make accessor objects friends.
                                      */
     template <int dim1, class DH> friend class DoFObjectAccessor;
+
+                                    /*
+                                     * Same with cell accessor
+                                     */
+    template <class DH> friend class DoFCellAccessor;
 };
 
 

diff --git a/deal.II/deal.II/include/dofs/dof_levels.h b/deal.II/deal.II/include/dofs/dof_levels.h
index bbfb1f994e443e3642816d5e7a36588fbe375eeb..a361c9a0fde4d4a5c5affc8d330de97aadec435c 100644 (file)
--- a/deal.II/deal.II/include/dofs/dof_levels.h
+++ b/deal.II/deal.II/include/dofs/dof_levels.h
@@ -46,17 +46,24 @@ namespace internal
  *
  * The storage format of the degrees of freedom indices (short: DoF
  * indices) is somewhat like a mirror of the data structures of the
- * triangulation classes.  There is a hierarchy of
- * DoFLevel<dim> classes for the different dimensions which
- * have objects named @p line_dofs, @p quad_dofs and so on, in which
- * the indices of DoFs located on lines and quads, respectively, are
- * stored. The indices are stored levelwise. The layout in these
- * arrays is as follows: if for a selected finite element (use
- * DoFHandler::distribute_dofs() to select a finite element)
- * the number of DoFs on each line (without those in the vertices) is
- * @p N, then the length of the @p line_dofs array is @p N times
- * the number of lines on this level. The DoF indices for the @p ith
- * line are at the positions <tt>N*i...(N+1)*i-1</tt>.
+ * triangulation classes.  There is a hierarchy of DoFLevel<dim>
+ * classes for the different dimensions which have objects named @p
+ * line_dofs, @p quad_dofs and so on, in which the indices of DoFs
+ * located on lines and quads, respectively, are stored. The indices
+ * are stored levelwise. The layout in these arrays is as follows: if
+ * for a selected finite element (use DoFHandler::distribute_dofs() to
+ * select a finite element) the number of DoFs on each line (without
+ * those in the vertices) is @p N, then the length of the @p line_dofs
+ * array is @p N times the number of lines on this level. The DoF
+ * indices for the @p ith line are at the positions
+ * <tt>N*i...(N+1)*i-1</tt>. It is true that this scheme wastes some
+ * memory since we actually only use the memory so allocated for
+ * active lines, whereas the index entries for line i remain unused if
+ * this line has children. On the other hand, it is easy to see that
+ * at least for dim>1 does not constitute a too big fraction of memory
+ * entries, and we opt for higher memory consumption in return for
+ * less index computations that usually slow down modern processors
+ * due to more complex memory access patterns.
  *
  * Since vertices are not associated with a particular level, the
  * indices associated with vertices are not stored in the DoFLevel
@@ -78,6 +85,43 @@ namespace internal
     };
 
 
+
+/**
+ * Common base class for all the DoFLevel classes. We here store
+ * information that is associated with (logical) cells, rather than
+ * concrete objects such as lines, quads, or hexes.
+ *
+ * At present, all we store are cached values for the DoF indices on
+ * each cell, since this is a frequently requested operation. The
+ * values are set by DoFCellAccessor::update_cell_dof_indices_cache
+ * and are used by DoFCellAccessor::get_dof_indices. Note that for
+ * objects that are not cell, such as lines in 2d and 3d, or quads in
+ * 3d, no caching is done since accessing the degrees of freedom on
+ * these objects is a comparatively rare operation.
+ */
+    template <>
+    class DoFLevel<0>
+    {
+      public:
+                                        /**
+                                         * Cache for the DoF indices
+                                         * on cells. The size of this
+                                         * array equals the number of
+                                         * cells on a given level
+                                         * times
+                                         * selected_fe.dofs_per_cell.
+                                         */
+       std::vector<unsigned int> cell_dof_indices_cache;
+
+                                         /**
+                                          * Determine an estimate for the
+                                          * memory consumption (in bytes)
+                                          * of this object.
+                                          */
+        unsigned int memory_consumption () const;
+    };
+    
+    
 /**
  * Store the indices of the degrees of freedom which are located on
  * the lines. See the general template DoFLevel for more information.
@@ -86,7 +130,7 @@ namespace internal
  * @author Wolfgang Bangerth, 1998
  */
     template <>
-    class DoFLevel<1>
+    class DoFLevel<1> : public DoFLevel<0>
     {
       public:
                                          /**

diff --git a/deal.II/deal.II/source/dofs/dof_accessor.cc b/deal.II/deal.II/source/dofs/dof_accessor.cc
index a3f502d19b85b9219313f65c643a6fccec98f5cb..8c57ae088c0b9d887904523372635e90e6c1ab1f 100644 (file)
--- a/deal.II/deal.II/source/dofs/dof_accessor.cc
+++ b/deal.II/deal.II/source/dofs/dof_accessor.cc
@@ -432,6 +432,68 @@ DoFCellAccessor<DoFHandler<1> >::face (const unsigned int) const
   return TriaIterator<1, DoFObjectAccessor<0,DoFHandler<1> > >();
 }
 
+
+template <>
+void
+DoFCellAccessor<DoFHandler<1> >::update_cell_dof_indices_cache () const
+{
+  Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
+          ExcMessage ("DoFHandler not initialized"));
+  
+  Assert (this->dof_handler != 0, DoFAccessor<DoFHandler<1> >::ExcInvalidObject());
+  Assert (&this->get_fe() != 0, DoFAccessor<DoFHandler<1> >::ExcInvalidObject());
+
+                                  // check as in documentation that
+                                  // cell is either active, or dofs
+                                  // are only in vertices. otherwise
+                                  // simply don't update the cache at
+                                  // all. the get_dof_indices
+                                  // function will then make sure we
+                                  // don't access the invalid data
+  if (this->has_children()
+      &&
+      (this->get_fe().dofs_per_cell !=
+       this->get_fe().dofs_per_vertex * GeometryInfo<1>::vertices_per_cell))
+    return;
+  
+  const unsigned int dofs_per_vertex = this->get_fe().dofs_per_vertex,
+                    dofs_per_line   = this->get_fe().dofs_per_line,
+                    dofs_per_cell   = this->get_fe().dofs_per_cell;
+
+                                  // make sure the cache is at least
+                                  // as big as we need it when
+                                  // writing to the last element of
+                                  // this cell
+  Assert (this->present_index * dofs_per_cell + dofs_per_cell
+         <=
+         this->dof_handler->levels[this->present_level]
+         ->cell_dof_indices_cache.size(),
+         ExcInternalError());
+
+  std::vector<unsigned int>::iterator next
+    = (this->dof_handler->levels[this->present_level]
+       ->cell_dof_indices_cache.begin() + this->present_index * dofs_per_cell);
+  
+  for (unsigned int vertex=0; vertex<2; ++vertex)
+    for (unsigned int d=0; d<dofs_per_vertex; ++d)
+      *next++ = vertex_dof_index(vertex,d);
+  for (unsigned int d=0; d<dofs_per_line; ++d)
+    *next++ = dof_index(d);
+}
+
+
+
+template <>
+void
+DoFCellAccessor<hp::DoFHandler<1> >::update_cell_dof_indices_cache () const
+{
+                                  // not implemented, but should also
+                                  // not be called
+  Assert (false, ExcNotImplemented());
+}
+
+
+
 #endif
 
 
@@ -444,6 +506,72 @@ DoFCellAccessor<DoFHandler<2> >::face (const unsigned int i) const
   return this->line(i);
 }
 
+
+
+template <>
+void
+DoFCellAccessor<DoFHandler<2> >::update_cell_dof_indices_cache () const
+{
+  Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
+          ExcMessage ("DoFHandler not initialized"));
+
+  Assert (this->dof_handler != 0, DoFAccessor<DoFHandler<2> >::ExcInvalidObject());
+  Assert (&this->get_fe() != 0, DoFAccessor<DoFHandler<2> >::ExcInvalidObject());
+
+                                  // check as in documentation that
+                                  // cell is either active, or dofs
+                                  // are only in vertices. otherwise
+                                  // simply don't update the cache at
+                                  // all. the get_dof_indices
+                                  // function will then make sure we
+                                  // don't access the invalid data
+  if (this->has_children()
+      &&
+      (this->get_fe().dofs_per_cell !=
+       this->get_fe().dofs_per_vertex * GeometryInfo<2>::vertices_per_cell))
+    return;
+  
+  const unsigned int dofs_per_vertex = this->get_fe().dofs_per_vertex,
+                    dofs_per_line   = this->get_fe().dofs_per_line,
+                    dofs_per_quad   = this->get_fe().dofs_per_quad,
+                    dofs_per_cell   = this->get_fe().dofs_per_cell;
+
+                                  // make sure the cache is at least
+                                  // as big as we need it when
+                                  // writing to the last element of
+                                  // this cell
+  Assert (this->present_index * dofs_per_cell + dofs_per_cell
+         <=
+         this->dof_handler->levels[this->present_level]
+         ->cell_dof_indices_cache.size(),
+         ExcInternalError());
+
+  std::vector<unsigned int>::iterator next
+    = (this->dof_handler->levels[this->present_level]
+       ->cell_dof_indices_cache.begin() + this->present_index * dofs_per_cell);
+
+  for (unsigned int vertex=0; vertex<4; ++vertex)
+    for (unsigned int d=0; d<dofs_per_vertex; ++d)
+      *next++ = vertex_dof_index(vertex,d);
+  for (unsigned int line=0; line<4; ++line)
+    for (unsigned int d=0; d<dofs_per_line; ++d)
+      *next++ = this->line(line)->dof_index(d);
+  for (unsigned int d=0; d<dofs_per_quad; ++d)
+    *next++ = dof_index(d);
+}
+
+
+
+template <>
+void
+DoFCellAccessor<hp::DoFHandler<2> >::update_cell_dof_indices_cache () const
+{
+                                  // not implemented, but should also
+                                  // not be called
+  Assert (false, ExcNotImplemented());
+}
+
+
 #endif
 
 
@@ -456,6 +584,75 @@ DoFCellAccessor<DoFHandler<3> >::face (const unsigned int i) const
   return this->quad(i);
 }
 
+
+
+template <>
+void
+DoFCellAccessor<DoFHandler<3> >::update_cell_dof_indices_cache () const
+{
+  Assert (static_cast<unsigned int>(this->present_level) < this->dof_handler->levels.size(),
+          ExcMessage ("DoFHandler not initialized"));
+
+  Assert (this->dof_handler != 0, DoFAccessor<DoFHandler<3> >::ExcInvalidObject());
+  Assert (&this->get_fe() != 0, DoFAccessor<DoFHandler<3> >::ExcInvalidObject());
+
+                                  // check as in documentation that
+                                  // cell is either active, or dofs
+                                  // are only in vertices. otherwise
+                                  // simply don't update the cache at
+                                  // all. the get_dof_indices
+                                  // function will then make sure we
+                                  // don't access the invalid data
+  if (this->has_children()
+      &&
+      (this->get_fe().dofs_per_cell !=
+       this->get_fe().dofs_per_vertex * GeometryInfo<3>::vertices_per_cell))
+    return;
+  
+  const unsigned int dofs_per_vertex = this->get_fe().dofs_per_vertex,
+                    dofs_per_line   = this->get_fe().dofs_per_line,
+                    dofs_per_quad   = this->get_fe().dofs_per_quad,
+                    dofs_per_hex    = this->get_fe().dofs_per_hex,
+                    dofs_per_cell   = this->get_fe().dofs_per_cell;
+
+                                  // make sure the cache is at least
+                                  // as big as we need it when
+                                  // writing to the last element of
+                                  // this cell
+  Assert (this->present_index * dofs_per_cell + dofs_per_cell
+         <=
+         this->dof_handler->levels[this->present_level]
+         ->cell_dof_indices_cache.size(),
+         ExcInternalError());
+
+  std::vector<unsigned int>::iterator next
+    = (this->dof_handler->levels[this->present_level]
+       ->cell_dof_indices_cache.begin() + this->present_index * dofs_per_cell);
+
+  for (unsigned int vertex=0; vertex<8; ++vertex)
+    for (unsigned int d=0; d<dofs_per_vertex; ++d)
+      *next++ = vertex_dof_index(vertex,d);
+  for (unsigned int line=0; line<12; ++line)
+    for (unsigned int d=0; d<dofs_per_line; ++d)
+      *next++ = this->line(line)->dof_index(d);
+  for (unsigned int quad=0; quad<6; ++quad)
+    for (unsigned int d=0; d<dofs_per_quad; ++d)
+      *next++ = this->quad(quad)->dof_index(d);
+  for (unsigned int d=0; d<dofs_per_hex; ++d)
+    *next++ = dof_index(d);
+}
+
+
+template <>
+void
+DoFCellAccessor<hp::DoFHandler<3> >::update_cell_dof_indices_cache () const
+{
+                                  // not implemented, but should also
+                                  // not be called
+  Assert (false, ExcNotImplemented());
+}
+
+
 #endif
 
 

diff --git a/deal.II/deal.II/source/dofs/dof_handler.cc b/deal.II/deal.II/source/dofs/dof_handler.cc
index be621d115c59e8aac7e5ea6037b877a14a86c45b..032b5e8afb3fe949262f4eda7fc4720c527bce5a 100644 (file)
--- a/deal.II/deal.II/source/dofs/dof_handler.cc
+++ b/deal.II/deal.II/source/dofs/dof_handler.cc
@@ -1457,13 +1457,20 @@ void DoFHandler<dim>::distribute_dofs (const FiniteElement<dim> &ff,
 
   used_dofs = next_free_dof;
 
+                                  // update the cache used
+                                  // for cell dof indices
+  for (cell_iterator cell = begin(); cell != end(); ++cell)
+    cell->update_cell_dof_indices_cache ();
+  
                                   // finally restore the user flags
   const_cast<Triangulation<dim> &>(*tria).load_user_flags(user_flags);
 }
 
 
+
 template <int dim>
-void DoFHandler<dim>::clear () {
+void DoFHandler<dim>::clear ()
+{
                                   // release lock to old fe
   selected_fe = 0;
 
@@ -1682,6 +1689,11 @@ void DoFHandler<1>::renumber_dofs (const std::vector<unsigned int> &new_numbers)
         i!=levels[level]->line_dofs.end(); ++i)
       if (*i != invalid_dof_index)
        *i = new_numbers[*i];
+
+                                  // finally update the cache used
+                                  // for cell dof indices
+  for (cell_iterator cell = begin(); cell != end(); ++cell)
+    cell->update_cell_dof_indices_cache ();
 }
 
 #endif
@@ -1735,7 +1747,12 @@ void DoFHandler<2>::renumber_dofs (const std::vector<unsigned int> &new_numbers)
           i!=levels[level]->quad_dofs.end(); ++i)
        if (*i != invalid_dof_index)
          *i = new_numbers[*i];
-    };
+    }
+
+                                  // finally update the cache used
+                                  // for cell dof indices
+  for (cell_iterator cell = begin(); cell != end(); ++cell)
+    cell->update_cell_dof_indices_cache ();
 }
 
 #endif
@@ -1793,7 +1810,12 @@ void DoFHandler<3>::renumber_dofs (const std::vector<unsigned int> &new_numbers)
           i!=levels[level]->hex_dofs.end(); ++i)
        if (*i != invalid_dof_index)
          *i = new_numbers[*i];
-    };
+    }
+
+                                  // finally update the cache used
+                                  // for cell dof indices
+  for (cell_iterator cell = begin(); cell != end(); ++cell)
+    cell->update_cell_dof_indices_cache ();
 }
 
 #endif
@@ -1967,7 +1989,8 @@ DoFHandler<3>::max_couplings_between_boundary_dofs () const
 #if deal_II_dimension == 1
 
 template <>
-void DoFHandler<1>::reserve_space () {
+void DoFHandler<1>::reserve_space ()
+{
   Assert (selected_fe != 0, ExcNoFESelected());
   Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
 
@@ -1984,10 +2007,14 @@ void DoFHandler<1>::reserve_space () {
     {
       levels.push_back (new internal::DoFHandler::DoFLevel<1>);
 
-      levels.back()->line_dofs = std::vector<unsigned int>(tria->n_raw_lines(i) *
-                                                          selected_fe->dofs_per_line,
-                                                          invalid_dof_index);
-    };
+      levels.back()->line_dofs.resize (tria->n_raw_lines(i) *
+                                      selected_fe->dofs_per_line,
+                                      invalid_dof_index);
+
+      levels.back()->cell_dof_indices_cache.resize (tria->n_raw_lines(i) *
+                                                   selected_fe->dofs_per_cell,
+                                                   invalid_dof_index);
+    }
 }
 
 
@@ -1997,7 +2024,8 @@ void DoFHandler<1>::reserve_space () {
 #if deal_II_dimension == 2
 
 template <>
-void DoFHandler<2>::reserve_space () {
+void DoFHandler<2>::reserve_space ()
+{
   Assert (selected_fe != 0, ExcNoFESelected());
   Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
   
@@ -2014,13 +2042,17 @@ void DoFHandler<2>::reserve_space () {
     {
       levels.push_back (new internal::DoFHandler::DoFLevel<2>);
 
-      levels.back()->line_dofs = std::vector<unsigned int> (tria->n_raw_lines(i) *
-                                                           selected_fe->dofs_per_line,
-                                                           invalid_dof_index);
-      levels.back()->quad_dofs = std::vector<unsigned int> (tria->n_raw_quads(i) *
-                                                           selected_fe->dofs_per_quad,
-                                                           invalid_dof_index);
-    };
+      levels.back()->line_dofs.resize (tria->n_raw_lines(i) *
+                                      selected_fe->dofs_per_line,
+                                      invalid_dof_index);
+      levels.back()->quad_dofs.resize (tria->n_raw_quads(i) *
+                                      selected_fe->dofs_per_quad,
+                                      invalid_dof_index);
+
+      levels.back()->cell_dof_indices_cache.resize (tria->n_raw_quads(i) *
+                                                   selected_fe->dofs_per_cell,
+                                                   invalid_dof_index);
+    }
 }
 
 #endif
@@ -2029,7 +2061,8 @@ void DoFHandler<2>::reserve_space () {
 #if deal_II_dimension == 3
 
 template <>
-void DoFHandler<3>::reserve_space () {
+void DoFHandler<3>::reserve_space ()
+{
   Assert (selected_fe != 0, ExcNoFESelected());
   Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
   
@@ -2055,7 +2088,11 @@ void DoFHandler<3>::reserve_space () {
       levels.back()->hex_dofs = std::vector<unsigned int> (tria->n_raw_hexs(i) *
                                                           selected_fe->dofs_per_hex,
                                                           invalid_dof_index);
-    };
+
+      levels.back()->cell_dof_indices_cache.resize (tria->n_raw_hexs(i) *
+                                                   selected_fe->dofs_per_cell,
+                                                   invalid_dof_index);
+    }
 }
 
 #endif

diff --git a/deal.II/deal.II/source/dofs/dof_levels.cc b/deal.II/deal.II/source/dofs/dof_levels.cc
index 3cf22e4417743e0c786a98ee949dce373902a7fc..2f97a2661b7cb822ec52ca54ecc5b50f1d6aac79 100644 (file)
--- a/deal.II/deal.II/source/dofs/dof_levels.cc
+++ b/deal.II/deal.II/source/dofs/dof_levels.cc
@@ -21,10 +21,18 @@ namespace internal
   namespace DoFHandler
   {
     
+    unsigned int
+    DoFLevel<0>::memory_consumption () const
+    {
+      return MemoryConsumption::memory_consumption (cell_dof_indices_cache);
+    }
+
+
     unsigned int
     DoFLevel<1>::memory_consumption () const
     {
-      return MemoryConsumption::memory_consumption (line_dofs);
+      return (DoFLevel<0>::memory_consumption () +
+              MemoryConsumption::memory_consumption (line_dofs));
     }
 
 

diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html
index 6cc91c7e7634587ee2496ac1ac7fa41cdfc1ee54..d48dad6febec60dae0a27fc2ff04079c51d92c14 100644 (file)
--- a/deal.II/doc/news/changes.html
+++ b/deal.II/doc/news/changes.html
@@ -552,6 +552,15 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p>
+       New: Since calling <code>cell-&gt;get_dof_indices</code> is a
+       fairly frequent operation (called about 7 times per cell in
+       step-14, but much more often in more complicated programs), the
+       results of this function are now cached for faster access.
+       <br> 
+       (WB 2006/03/29)
+       </p>
+
   <li> <p>
        Fixed: An exception was generated when trying to interpolate
        boundary values for a primitive component of a composite finite