From 80fd6fb73f81e8021a37003b20b9522bd17a320a Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Mon, 9 Sep 2013 19:10:18 +0000
Subject: [PATCH] Implement the compression scheme for DoF indices stored on
 cells.

git-svn-id: https://svn.dealii.org/trunk@30655 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/include/deal.II/hp/dof_handler.h |   5 +-
 deal.II/include/deal.II/hp/dof_level.h   | 134 +++++++++++----
 deal.II/include/deal.II/hp/dof_objects.h |   2 +-
 deal.II/source/hp/dof_handler.cc         |  60 ++++---
 deal.II/source/hp/dof_level.cc           | 209 ++++++++++++++++++++++-
 5 files changed, 345 insertions(+), 65 deletions(-)

diff --git a/deal.II/include/deal.II/hp/dof_handler.h b/deal.II/include/deal.II/hp/dof_handler.h
index d4ff6b16f7..cf416021fd 100644
--- a/deal.II/include/deal.II/hp/dof_handler.h
+++ b/deal.II/include/deal.II/hp/dof_handler.h
@@ -38,7 +38,7 @@ namespace internal
 {
   namespace hp
   {
-    template <int> class DoFLevel;
+    class DoFLevel;
     template <int> class DoFFaces;
     template <int> class DoFObjects;
 
@@ -790,7 +790,7 @@ namespace hp
      * <tt>levels[]</tt> tree of
      * the Triangulation objects.
      */
-    std::vector<dealii::internal::hp::DoFLevel<dim>*> levels;
+    std::vector<dealii::internal::hp::DoFLevel*> levels;
 
     /**
      * Space to store the DoF
@@ -900,7 +900,6 @@ namespace hp
      * the functions that set and
      * retrieve vertex dof indices.
      */
-    template <int> friend class dealii::internal::hp::DoFLevel;
     template <int> friend class dealii::internal::hp::DoFObjects;
     friend struct dealii::internal::hp::DoFHandler::Implementation;
   };
diff --git a/deal.II/include/deal.II/hp/dof_level.h b/deal.II/include/deal.II/hp/dof_level.h
index 343d81eb09..58da6c6c16 100644
--- a/deal.II/include/deal.II/hp/dof_level.h
+++ b/deal.II/include/deal.II/hp/dof_level.h
@@ -29,6 +29,7 @@ DEAL_II_NAMESPACE_OPEN
 namespace hp
 {
   template <int, int> class DoFHandler;
+  template <int, int> class FECollection;
 }
 
 
@@ -56,11 +57,64 @@ namespace internal
      * where each cell's indices start within the array of indices. This is in
      * contrast to the DoFObjects class where each face or edge may have more than
      * one associated finite element with corresponding degrees of freedom.
+     *
+     * The data stored here is represented by three arrays
+     * - The @p active_fe_indices array stores for each cell which
+     *   finite element is used on this cell. Since some cells are not active
+     *   on the current level, some entries in this array may represent an
+     *   invalid value.
+     * - The @p dof_indices array stores for each active cell on the current
+     *   level the dofs that associated with the <i>interior</i> of the cell,
+     *   i.e., the @p dofs_per_line dofs associated with the line in 1d, and
+     *   @p dofs_per_quad and @p dofs_per_hex in 2d and 3d. These numbers are
+     *   in general smaller than @p dofs_per_cell.
+     * - The @p dof_offsets array stores, for each cell, the starting point
+     *   of the dof indices corresponding to this cell in the @p dof_indices
+     *   array. This is analogous to how we store data in compressed row storage
+     *   for sparse matrices. For cells that are not active on the current level,
+     *   we store an invalid value for the starting index.
+     *
+     * <h3>Compression</h3>
+     *
+     * It is common for the indices stored in @p dof_indices for one cell to be
+     * numbered consecutively. For example, using the standard numbering (without
+     * renumbering DoFs), the quad dofs on the first cell of a mesh when using a
+     * $Q_3$ element will be numbered <tt>12, 13, 14, 15</tt>. This allows for
+     * compression if we only store the first entry and have some way to mark
+     * the DoFs on this object as compressed. Here, compression means that we
+     * know that subsequent DoF indices can be obtained from the previous ones
+     * by just incrementing them by one -- in other words, we use a variant of doing
+     * run-length encoding. The way to do this is that we
+     * use positive FE indices for uncompressed sets of DoFs and if a set of indices
+     * is compressed, then we instead store the FE index in binary complement (which
+     * we can identify by looking at the sign bit when interpreting the number as a
+     * signed one). There are two functions, compress_data() and uncompress_data()
+     * that convert between the two possible representations.
+     *
+     * Note that compression is not always possible. For example, if one renumbered
+     * the example above using DoFRenumbering::downstream with $(1,0)^T$ as
+     * direction, then they would likely be numbered <tt>12, 14, 13, 15</tt>, which
+     * can not be compressed using run-length encoding.
      */
-    template <int dim>
     class DoFLevel
     {
     private:
+      /**
+       * The type in which we store the offsets into the dof_indices array.
+       */
+      typedef unsigned int offset_type;
+
+      /**
+       * The type in which we store the active FE index.
+       */
+      typedef unsigned short int active_fe_index_type;
+
+      /**
+       * A signed type that matches the type in which we store the active FE
+       * index. We use this in computing binary complements.
+       */
+      typedef signed short int signed_active_fe_index_type;
+
       /**
        * Indices specifying the finite element of hp::FECollection to
        * use for the different cells on the current level. The vector
@@ -72,7 +126,7 @@ namespace internal
        * it does not have an associated fe index and we store
        * an invalid fe index marker instead.
        */
-      std::vector<unsigned int> active_fe_indices;
+      std::vector<active_fe_index_type> active_fe_indices;
 
       /**
        * Store the start index for the degrees of freedom of each
@@ -86,7 +140,7 @@ namespace internal
        * The type we store is then obviously the type the @p dof_indices array
        * uses for indexing.
        */
-      std::vector<std::vector<types::global_dof_index>::size_type> dof_offsets;
+      std::vector<offset_type> dof_offsets;
 
       /**
        * Store the global indices of the degrees of freedom.
@@ -182,6 +236,7 @@ namespace internal
       set_active_fe_index (const unsigned int obj_index,
 			   const unsigned int fe_index);
 
+
       /**
        * Determine an estimate for the
        * memory consumption (in bytes)
@@ -189,6 +244,31 @@ namespace internal
        */
       std::size_t memory_consumption () const;
 
+    private:
+      /**
+       * Compress the arrays that store dof indices by using a variant
+       * of run-length encoding. See the general documentation of this
+       * class for more information.
+       *
+       * @param fe_collection The object that can tell us how many
+       * degrees of freedom each of the finite elements has that we
+       * store in this object.
+       */
+      template <int dim, int spacedim>
+      void compress_data (const dealii::hp::FECollection<dim,spacedim> &fe_collection);
+
+      /**
+       * Uncompress the arrays that store dof indices by using a variant
+       * of run-length encoding. See the general documentation of this
+       * class for more information.
+       *
+       * @param fe_collection The object that can tell us how many
+       * degrees of freedom each of the finite elements has that we
+       * store in this object.
+       */
+      template <int dim, int spacedim>
+      void uncompress_data (const dealii::hp::FECollection<dim,spacedim> &fe_collection);
+
       /**
        * Make hp::DoFHandler and its auxiliary class a friend since it
        * is the class that needs to create these data structures.
@@ -200,10 +280,9 @@ namespace internal
 
     // -------------------- template functions --------------------------------
 
-    template <int dim>
     inline
     types::global_dof_index
-    DoFLevel<dim>::
+    DoFLevel::
     get_dof_index (const unsigned int                obj_index,
                    const unsigned int                fe_index,
                    const unsigned int                local_index) const
@@ -221,15 +300,20 @@ namespace internal
 
       Assert (fe_index == active_fe_indices[obj_index],
 	      ExcMessage ("FE index does not match that of the present cell"));
-      return dof_indices[dof_offsets[obj_index]+local_index];
+
+      // see if the dof_indices array has been compressed for this
+      // particular cell
+      if ((signed_active_fe_index_type)active_fe_indices[obj_index]>=0)
+        return dof_indices[dof_offsets[obj_index]+local_index];
+      else
+        return dof_indices[dof_offsets[obj_index]]+local_index;
     }
 
 
 
-    template <int dim>
     inline
     void
-    DoFLevel<dim>::
+    DoFLevel::
     set_dof_index (const unsigned int                obj_index,
                    const unsigned int                fe_index,
                    const unsigned int                local_index,
@@ -245,7 +329,8 @@ namespace internal
               ExcMessage ("You are trying to access degree of freedom "
                           "information for an object on which no such "
                           "information is available"));
-
+      Assert ((signed_active_fe_index_type)active_fe_indices[obj_index]>=0,
+              ExcMessage ("This function can no longer be called after compressing the dof_indices array"));
       Assert (fe_index == active_fe_indices[obj_index],
 	      ExcMessage ("FE index does not match that of the present cell"));
       dof_indices[dof_offsets[obj_index]+local_index] = global_index;
@@ -253,50 +338,35 @@ namespace internal
 
 
 
-    template <int dim>
     inline
     unsigned int
-    DoFLevel<dim>::
+    DoFLevel::
     active_fe_index (const unsigned int obj_index) const
     {
       Assert (obj_index < active_fe_indices.size(),
               ExcIndexRange (obj_index, 0, active_fe_indices.size()));
 
-      return active_fe_indices[obj_index];
+      if (((signed_active_fe_index_type)active_fe_indices[obj_index]) >= 0)
+        return active_fe_indices[obj_index];
+      else
+        return (active_fe_index_type)~(signed_active_fe_index_type)active_fe_indices[obj_index];
     }
 
 
 
-    template <int dim>
     inline
     bool
-    DoFLevel<dim>::
+    DoFLevel::
     fe_index_is_active (const unsigned int                obj_index,
                         const unsigned int                fe_index) const
     {
-      const unsigned int invalid_fe_index = numbers::invalid_unsigned_int;
-      Assert ((fe_index != invalid_fe_index),
-              ExcMessage ("You need to specify a FE index when working "
-                          "with hp DoFHandlers"));
-
-      // make sure we are on an
-      // object for which DoFs have
-      // been allocated at all
-      Assert (dof_offsets[obj_index] != numbers::invalid_dof_index,
-              ExcMessage ("You are trying to access degree of freedom "
-                          "information for an object on which no such "
-                          "information is available"));
-
-      Assert (obj_index < active_fe_indices.size(),
-	      ExcInternalError());
-      return (fe_index == active_fe_indices[obj_index]);
+      return (fe_index == active_fe_index(obj_index));
     }
 
 
-    template <int dim>
     inline
     void
-    DoFLevel<dim>::
+    DoFLevel::
     set_active_fe_index (const unsigned int obj_index,
 			 const unsigned int fe_index)
     {
diff --git a/deal.II/include/deal.II/hp/dof_objects.h b/deal.II/include/deal.II/hp/dof_objects.h
index 59ec451d53..cb21ed7701 100644
--- a/deal.II/include/deal.II/hp/dof_objects.h
+++ b/deal.II/include/deal.II/hp/dof_objects.h
@@ -106,7 +106,7 @@ namespace internal
        * The type we store is then obviously the type the @p dofs array
        * uses for indexing.
        */
-      std::vector<std::vector<types::global_dof_index>::size_type> dof_offsets;
+      std::vector<unsigned int> dof_offsets;
 
       /**
        * Store the global indices of
diff --git a/deal.II/source/hp/dof_handler.cc b/deal.II/source/hp/dof_handler.cc
index 8e2bc5d976..318a1c3d7f 100644
--- a/deal.II/source/hp/dof_handler.cc
+++ b/deal.II/source/hp/dof_handler.cc
@@ -539,7 +539,7 @@ namespace internal
           // active_fe_indices field which
           // we have to backup before
           {
-            std::vector<std::vector<unsigned int> >
+            std::vector<std::vector<DoFLevel::active_fe_index_type> >
             active_fe_backup(dof_handler.levels.size ());
             for (unsigned int level = 0; level<dof_handler.levels.size (); ++level)
               std::swap (dof_handler.levels[level]->active_fe_indices,
@@ -553,7 +553,7 @@ namespace internal
 
             for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
               {
-                dof_handler.levels.push_back (new internal::hp::DoFLevel<dim>);
+                dof_handler.levels.push_back (new internal::hp::DoFLevel);
                 std::swap (active_fe_backup[level],
                            dof_handler.levels[level]->active_fe_indices);
               }
@@ -578,9 +578,9 @@ namespace internal
           for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
             {
               dof_handler.levels[level]->dof_offsets
-                = std::vector<std::vector<types::global_dof_index>::size_type> (
+                = std::vector<DoFLevel::offset_type> (
                     dof_handler.tria->n_raw_lines(level),
-                    (std::vector<types::global_dof_index>::size_type)(-1));
+                    (DoFLevel::offset_type)(-1));
 
               types::global_dof_index next_free_dof = 0;
               for (typename DoFHandler<dim,spacedim>::active_cell_iterator
@@ -619,7 +619,7 @@ namespace internal
               Assert (static_cast<unsigned int>
                       (std::count (dof_handler.levels[level]->dof_offsets.begin(),
                                    dof_handler.levels[level]->dof_offsets.end(),
-                                   (std::vector<types::global_dof_index>::size_type)(-1)))
+                                   (DoFLevel::offset_type)(-1)))
                       ==
                       dof_handler.tria->n_raw_lines(level) - dof_handler.tria->n_active_lines(level),
                       ExcInternalError());
@@ -654,7 +654,7 @@ namespace internal
           // active_fe_indices field which
           // we have to backup before
           {
-            std::vector<std::vector<unsigned int> >
+            std::vector<std::vector<DoFLevel::active_fe_index_type> >
             active_fe_backup(dof_handler.levels.size ());
             for (unsigned int level = 0; level<dof_handler.levels.size (); ++level)
               std::swap (dof_handler.levels[level]->active_fe_indices,
@@ -668,7 +668,7 @@ namespace internal
 
             for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
               {
-                dof_handler.levels.push_back (new internal::hp::DoFLevel<dim>);
+                dof_handler.levels.push_back (new internal::hp::DoFLevel);
                 std::swap (active_fe_backup[level],
                            dof_handler.levels[level]->active_fe_indices);
               }
@@ -695,9 +695,9 @@ namespace internal
           for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
             {
               dof_handler.levels[level]->dof_offsets
-                = std::vector<std::vector<types::global_dof_index>::size_type> (
+                = std::vector<DoFLevel::offset_type> (
                     dof_handler.tria->n_raw_quads(level),
-                    (std::vector<types::global_dof_index>::size_type)(-1));
+                    (DoFLevel::offset_type)(-1));
 
               types::global_dof_index next_free_dof = 0;
               for (typename DoFHandler<dim,spacedim>::active_cell_iterator
@@ -736,7 +736,7 @@ namespace internal
               Assert (static_cast<unsigned int>
                       (std::count (dof_handler.levels[level]->dof_offsets.begin(),
                                    dof_handler.levels[level]->dof_offsets.end(),
-                                   (std::vector<types::global_dof_index>::size_type)(-1)))
+                                   (DoFLevel::offset_type)(-1)))
                       ==
                       dof_handler.tria->n_raw_quads(level) - dof_handler.tria->n_active_quads(level),
                       ExcInternalError());
@@ -750,7 +750,7 @@ namespace internal
           // then allocate as much space as
           // we need and prime the linked
           // list for lines (see the
-          // description in hp::DoFLevels)
+          // description in hp::DoFLevel)
           // with the indices we will
           // need. note that our task is
           // more complicated since two
@@ -847,9 +847,8 @@ namespace internal
             // active ones will have a
             // non-invalid value later on
             dof_handler.faces->lines.dof_offsets
-              = std::vector<std::vector<types::global_dof_index>::size_type>
-                                                     (dof_handler.tria->n_raw_lines(),
-                                                     (std::vector<types::global_dof_index>::size_type)(-1));
+              = std::vector<unsigned int> (dof_handler.tria->n_raw_lines(),
+                                             (unsigned int)(-1));
             dof_handler.faces->lines.dofs
               = std::vector<types::global_dof_index> (n_line_slots,
                                                       DoFHandler<dim,spacedim>::invalid_dof_index);
@@ -1020,7 +1019,7 @@ namespace internal
           // active_fe_indices field which
           // we have to backup before
           {
-            std::vector<std::vector<unsigned int> >
+            std::vector<std::vector<DoFLevel::active_fe_index_type> >
             active_fe_backup(dof_handler.levels.size ());
             for (unsigned int level = 0; level<dof_handler.levels.size (); ++level)
               std::swap (dof_handler.levels[level]->active_fe_indices,
@@ -1034,7 +1033,7 @@ namespace internal
 
             for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
               {
-                dof_handler.levels.push_back (new internal::hp::DoFLevel<dim>);
+                dof_handler.levels.push_back (new internal::hp::DoFLevel);
                 std::swap (active_fe_backup[level],
                            dof_handler.levels[level]->active_fe_indices);
               }
@@ -1061,9 +1060,9 @@ namespace internal
           for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
             {
               dof_handler.levels[level]->dof_offsets
-                = std::vector<std::vector<types::global_dof_index>::size_type>
+                = std::vector<DoFLevel::offset_type>
                                                        (dof_handler.tria->n_raw_hexs(level),
-                                                       (std::vector<types::global_dof_index>::size_type)(-1));
+                                                       (DoFLevel::offset_type)(-1));
 
               types::global_dof_index next_free_dof = 0;
               for (typename DoFHandler<dim,spacedim>::active_cell_iterator
@@ -1102,7 +1101,7 @@ namespace internal
               Assert (static_cast<unsigned int>
                       (std::count (dof_handler.levels[level]->dof_offsets.begin(),
                                    dof_handler.levels[level]->dof_offsets.end(),
-                                   (std::vector<types::global_dof_index>::size_type)(-1)))
+                                   (DoFLevel::offset_type)(-1)))
                       ==
                       dof_handler.tria->n_raw_hexs(level) - dof_handler.tria->n_active_hexs(level),
                       ExcInternalError());
@@ -1116,7 +1115,7 @@ namespace internal
           // then allocate as much space as
           // we need and prime the linked
           // list for quad (see the
-          // description in hp::DoFLevels)
+          // description in hp::DoFLevel)
           // with the indices we will
           // need. note that our task is
           // more complicated since two
@@ -1217,9 +1216,9 @@ namespace internal
             if (true)
               {
                 dof_handler.faces->quads.dof_offsets
-                  = std::vector<std::vector<types::global_dof_index>::size_type>
+                  = std::vector<unsigned int>
                                                          (dof_handler.tria->n_raw_quads(),
-                                                         (std::vector<types::global_dof_index>::size_type)(-1));
+                                                         (unsigned int)(-1));
                 dof_handler.faces->quads.dofs
                   = std::vector<types::global_dof_index> (n_quad_slots,
                                                           DoFHandler<dim,spacedim>::invalid_dof_index);
@@ -2691,6 +2690,9 @@ namespace hp
       = std::vector<IndexSet> (1,
                                number_cache.locally_owned_dofs);
 
+    for (unsigned int level=0; level<levels.size(); ++level)
+      levels[level]->compress_data (*finite_elements);
+
     // finally restore the user flags
     const_cast<Triangulation<dim,spacedim> &>(*tria).load_user_flags(user_flags);
   }
@@ -2727,7 +2729,17 @@ namespace hp
       }
 #endif
 
+    // uncompress the internal storage scheme of dofs on cells
+    // so that we can access dofs in turns
+    for (unsigned int level=0; level<levels.size(); ++level)
+      levels[level]->uncompress_data (*finite_elements);
+
+    // do the renumbering
     renumber_dofs_internal (new_numbers, dealii::internal::int2type<dim>());
+
+    // now re-compress the dof indices
+    for (unsigned int level=0; level<levels.size(); ++level)
+      levels[level]->compress_data (*finite_elements);
   }
 
 
@@ -3088,7 +3100,7 @@ namespace hp
     // Create sufficiently many
     // hp::DoFLevels.
     while (levels.size () < tria->n_levels ())
-      levels.push_back (new dealii::internal::hp::DoFLevel<dim>);
+      levels.push_back (new dealii::internal::hp::DoFLevel);
 
     // then make sure that on each
     // level we have the appropriate
@@ -3164,7 +3176,7 @@ namespace hp
     // Normally only one level is added, but if this Triangulation
     // is created by copy_triangulation, it can be more than one level.
     while (levels.size () < tria->n_levels ())
-      levels.push_back (new dealii::internal::hp::DoFLevel<dim>);
+      levels.push_back (new dealii::internal::hp::DoFLevel);
 
     // Coarsening can lead to the loss
     // of levels. Hence remove them.
diff --git a/deal.II/source/hp/dof_level.cc b/deal.II/source/hp/dof_level.cc
index 02a47d2990..8215acd3da 100644
--- a/deal.II/source/hp/dof_level.cc
+++ b/deal.II/source/hp/dof_level.cc
@@ -16,6 +16,8 @@
 
 #include <deal.II/base/memory_consumption.h>
 #include <deal.II/hp/dof_level.h>
+#include <deal.II/hp/fe_collection.h>
+#include <iostream>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -23,9 +25,196 @@ namespace internal
 {
   namespace hp
   {
-    template <int dim>
+    template <int dim, int spacedim>
+    void
+    DoFLevel::compress_data (const dealii::hp::FECollection<dim,spacedim> &fe_collection)
+    {
+      return;
+
+      if (dof_offsets.size() == 0 || dof_indices.size()==0)
+        return;
+
+      // in a first run through, count how many new slots we need in the
+      // dof_indices array after compression. note that the 'cell'
+      // counter is incremented inside the loop
+      unsigned int new_size = 0;
+      for (unsigned int cell=0; cell<dof_offsets.size(); )
+        // see if this cell is active on the current level
+        if (dof_offsets[cell] != (offset_type)(-1))
+          {
+            // find the next cell active on this level
+            unsigned int next_cell = cell+1;
+            while ((next_cell<dof_offsets.size()) &&
+                    (dof_offsets[next_cell] == (offset_type)(-1)))
+              ++next_cell;
+
+            const unsigned int next_offset = (next_cell < dof_offsets.size() ?
+                                                dof_offsets[next_cell] :
+                                                dof_indices.size());
+
+            Assert (next_offset-dof_offsets[cell] == fe_collection[active_fe_indices[cell]].template n_dofs_per_object<dim>(),
+                    ExcInternalError());
+
+            // see if the range of dofs for this cell can be compressed and if so
+            // how many slots we have to store for them
+            if (next_offset > dof_offsets[cell])
+              {
+                bool compressible = true;
+                for (unsigned int j=dof_offsets[cell]+1; j<next_offset; ++j)
+                  if (dof_indices[j] != dof_indices[j-1]+1)
+                    {
+                      compressible = false;
+                      break;
+                    }
+                if (compressible == true)
+                  new_size += 1;
+                else
+                  new_size += (next_offset-dof_offsets[cell]);
+              }
+
+            // then move on to the next cell
+            cell = next_cell;
+          }
+        else
+          ++cell;
+
+      // now allocate the new array and copy into it whatever we need
+      std::vector<types::global_dof_index> new_dof_indices;
+      new_dof_indices.reserve(new_size);
+      for (unsigned int cell=0; cell<dof_offsets.size(); )
+        // see if this cell is active on the current level
+        if (dof_offsets[cell] != (offset_type)(-1))
+          {
+            // find the next cell active on this level
+            unsigned int next_cell = cell+1;
+            while ((next_cell<dof_offsets.size()) &&
+                    (dof_offsets[next_cell] == (offset_type)(-1)))
+              ++next_cell;
+
+            const unsigned int next_offset = (next_cell < dof_offsets.size() ?
+                                                dof_offsets[next_cell] :
+                                                dof_indices.size());
+
+            Assert (next_offset-dof_offsets[cell] == fe_collection[active_fe_indices[cell]].template n_dofs_per_object<dim>(),
+                    ExcInternalError());
+
+            // see if the range of dofs for this cell can be compressed and if so
+            // how many slots we have to store for them
+            if (next_offset > dof_offsets[cell])
+              {
+                bool compressible = true;
+                for (unsigned int j=dof_offsets[cell]+1; j<next_offset; ++j)
+                  if (dof_indices[j] != dof_indices[j-1]+1)
+                    {
+                      compressible = false;
+                      break;
+                    }
+
+                // if this cell is compressible, then copy the first index and mark this
+                // in the dof_offsets array
+                if (compressible == true)
+                  {
+                    new_dof_indices.push_back (dof_indices[dof_offsets[cell]]);
+
+                    // make sure that the current active_fe_index indicates
+                    // that this entry hasn't been compressed yet
+                    Assert ((signed_active_fe_index_type)active_fe_indices[cell] >= 0, ExcInternalError());
+
+                    // then mark the compression
+                    active_fe_indices[cell] = (active_fe_index_type)~(signed_active_fe_index_type)active_fe_indices[cell];
+                  }
+                else
+                  for (unsigned int i=dof_offsets[cell]; i<next_offset; ++i)
+                    new_dof_indices.push_back (dof_indices[i]);
+              }
+
+            // then move on to the next cell
+            cell = next_cell;
+          }
+        else
+          ++cell;
+
+      // finally swap old and new content
+      Assert (new_dof_indices.size() == new_size, ExcInternalError());
+      dof_indices.swap (new_dof_indices);
+    }
+
+
+
+    template <int dim, int spacedim>
+    void
+    DoFLevel::uncompress_data(const dealii::hp::FECollection<dim,spacedim> &fe_collection)
+    {
+      return;
+
+      if (dof_offsets.size() == 0 || dof_indices.size()==0)
+        return;
+
+      // in a first run through, count how many new slots we need in the
+      // dof_indices array after uncompression.
+      unsigned int new_size = 0;
+      for (unsigned int cell=0; cell<dof_offsets.size(); ++cell)
+        if (dof_offsets[cell] != (offset_type)(-1))
+          {
+            // we know now that the slot for this cell is used. extract the
+            // active_fe_index for it and see how many entries we need
+            new_size += fe_collection[active_fe_index(cell)].template n_dofs_per_object<dim>();
+          }
+
+      // now allocate the new array and copy into it whatever we need
+      std::vector<types::global_dof_index> new_dof_indices;
+      new_dof_indices.reserve(new_size);
+      std::vector<offset_type> new_dof_offsets (dof_offsets.size(), (offset_type)(-1));
+      for (unsigned int cell=0; cell<dof_offsets.size(); )
+        // see if this cell is active on the current level
+        if (dof_offsets[cell] != (offset_type)(-1))
+          {
+            // find the next cell active on this level
+            unsigned int next_cell = cell+1;
+            while ((next_cell<dof_offsets.size()) &&
+                    (dof_offsets[next_cell] == (offset_type)(-1)))
+              ++next_cell;
+
+            const unsigned int next_offset = (next_cell < dof_offsets.size() ?
+                                                dof_offsets[next_cell] :
+                                                dof_indices.size());
+
+            // set offset for this cell
+            new_dof_offsets[cell] = new_dof_indices.size();
+
+            // see if we need to uncompress this set of dofs
+            if ((signed_active_fe_index_type)active_fe_indices[cell]>=0)
+              {
+                // apparently not. simply copy them
+                Assert (next_offset-dof_offsets[cell] == fe_collection[active_fe_indices[cell]].template n_dofs_per_object<dim>(),
+                        ExcInternalError());
+                for (unsigned int i=dof_offsets[cell]; i<next_offset; ++i)
+                  new_dof_indices.push_back (dof_indices[i]);
+              }
+            else
+              {
+                // apparently so. uncompress
+                Assert (next_offset-dof_offsets[cell] == 1,
+                        ExcInternalError());
+                for (unsigned int i=0; i<fe_collection[active_fe_indices[cell]].template n_dofs_per_object<dim>(); ++i)
+                  new_dof_indices.push_back (dof_indices[dof_offsets[cell]]+i);
+              }
+
+            // then move on to the next cell
+            cell = next_cell;
+          }
+        else
+          ++cell;
+
+      // verify correct size, then swap arrays
+      Assert (new_dof_indices.size() == new_size, ExcInternalError());
+      dof_indices.swap (new_dof_indices);
+      dof_offsets.swap (new_dof_offsets);
+    }
+
+
     std::size_t
-    DoFLevel<dim>::memory_consumption () const
+    DoFLevel::memory_consumption () const
     {
       return (MemoryConsumption::memory_consumption (active_fe_indices) +
 	      MemoryConsumption::memory_consumption (dof_indices) +
@@ -34,9 +223,19 @@ namespace internal
 
 
     // explicit instantiations
-    template std::size_t DoFLevel<1>::memory_consumption () const;
-    template std::size_t DoFLevel<2>::memory_consumption () const;
-    template std::size_t DoFLevel<3>::memory_consumption () const;
+    template void DoFLevel::compress_data(const dealii::hp::FECollection<1,1> &);
+    template void DoFLevel::compress_data(const dealii::hp::FECollection<1,2> &);
+    template void DoFLevel::compress_data(const dealii::hp::FECollection<1,3> &);
+    template void DoFLevel::compress_data(const dealii::hp::FECollection<2,2> &);
+    template void DoFLevel::compress_data(const dealii::hp::FECollection<2,3> &);
+    template void DoFLevel::compress_data(const dealii::hp::FECollection<3,3> &);
+
+    template void DoFLevel::uncompress_data(const dealii::hp::FECollection<1,1> &);
+    template void DoFLevel::uncompress_data(const dealii::hp::FECollection<1,2> &);
+    template void DoFLevel::uncompress_data(const dealii::hp::FECollection<1,3> &);
+    template void DoFLevel::uncompress_data(const dealii::hp::FECollection<2,2> &);
+    template void DoFLevel::uncompress_data(const dealii::hp::FECollection<2,3> &);
+    template void DoFLevel::uncompress_data(const dealii::hp::FECollection<3,3> &);
   }
 }
 
-- 
2.39.5