]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix and enable compression of hp DoF indices.
authorWolfgang Bangerth <bangerth@colostate.edu>
Fri, 6 Jan 2017 15:14:03 +0000 (08:14 -0700)
committerWolfgang Bangerth <bangerth@colostate.edu>
Sun, 8 Jan 2017 01:50:44 +0000 (18:50 -0700)
A previous patch, a long time ago, implemented a scheme to compress the data structure
that stores the DoF indices in the hp case. For reasons no longer clear, the functions
that do the compression and decompression just returned immediately, without actually
doing anything. Maybe not surprisingly, the algorithm -- not executed after all -- did
not work.

This patch fixes the algorithm at various locations and enables the compression.

include/deal.II/hp/dof_level.h
source/hp/dof_handler.cc
source/hp/dof_level.cc

index 8c6316f7a3df1a7396dbc064acb0e3ed0864f8a3..0904bf58ebaed6b0a385b526fba49c11fcfdb2e7 100644 (file)
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 2003 - 2016 by the deal.II authors
+// Copyright (C) 2003 - 2017 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -119,6 +119,25 @@ namespace internal
        */
       typedef signed short int signed_active_fe_index_type;
 
+      /**
+       * Given an active_fe_index, return whether the corresponding
+       * set of DoF indices are compressed. See the general documentation
+       * of this class for a description of when this is the case.
+       */
+      static
+      bool
+      is_compressed_entry (const active_fe_index_type active_fe_index);
+
+      /**
+       * Given an active_fe_index (either corresponding to an uncompressed
+       * or compressed state), return the active_fe_index that corresponds
+       * to the respectively other state. See the general documentation
+       * of this class for a description of how compression is indicated.
+       */
+      static
+      active_fe_index_type
+      get_toggled_compression_state (const active_fe_index_type active_fe_index);
+
       /**
        * Indices specifying the finite element of hp::FECollection to use for
        * the different cells on the current level. The vector stores one
@@ -273,6 +292,27 @@ namespace internal
       template <int dim, int spacedim>
       void uncompress_data (const dealii::hp::FECollection<dim,spacedim> &fe_collection);
 
+
+      /**
+       * Restore the active fe indices stored by the current object to
+       * an uncompressed state. At the same time, do not touch any of
+       * the other data structures. This leaves the active_fe_indices
+       * array out-of-synch from the other member variables, and the
+       * function can consequently only be used when the remaining
+       * data structures are going to be rebuilt next. This is
+       * specifically the case in hp::DoFHandler::distribute_dofs()
+       * where we throw away all data in the DoFLevels objects *except
+       * for the active_fe_indices* array. This function therefore
+       * simply makes sure that that one array is in uncompressed
+       * format so that we can use the information about active
+       * fe_indices for all cells that were used in the previous mesh
+       * refinement cycle (or the previous time distribute_dofs() was
+       * called) without having to care about any of the other data
+       * fields.
+       */
+      void normalize_active_fe_indices ();
+
+
       /**
        * Make hp::DoFHandler and its auxiliary class a friend since it is the
        * class that needs to create these data structures.
@@ -285,6 +325,26 @@ namespace internal
 
     // -------------------- template functions --------------------------------
 
+
+    inline
+    bool DoFLevel::is_compressed_entry (const active_fe_index_type active_fe_index)
+    {
+      return ((signed_active_fe_index_type)active_fe_index < 0);
+    }
+
+
+
+    inline
+    DoFLevel::active_fe_index_type
+    DoFLevel::get_toggled_compression_state (const active_fe_index_type active_fe_index)
+    {
+      // convert the active_fe_index into a signed type, flip all
+      // bits, and get the unsigned representation back
+      return (active_fe_index_type)~(signed_active_fe_index_type)active_fe_index;
+    }
+
+
+
     inline
     types::global_dof_index
     DoFLevel::
@@ -296,20 +356,21 @@ namespace internal
       Assert (obj_index < dof_offsets.size(),
               ExcIndexRange (obj_index, 0, dof_offsets.size()));
 
-      // make sure we are on an
-      // object for which DoFs have
-      // been allocated at all
+      // make sure we are on an object for which DoFs have been
+      // allocated at all
       Assert (dof_offsets[obj_index] != (offset_type)(-1),
               ExcMessage ("You are trying to access degree of freedom "
                           "information for an object on which no such "
                           "information is available"));
 
-      Assert (fe_index == active_fe_indices[obj_index],
+      Assert (fe_index == (is_compressed_entry(active_fe_indices[obj_index]) == false ?
+                           active_fe_indices[obj_index] :
+                           get_toggled_compression_state(active_fe_indices[obj_index])),
               ExcMessage ("FE index does not match that of the present cell"));
 
       // 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)
+      if (is_compressed_entry(active_fe_indices[obj_index]) == false)
         return dof_indices[dof_offsets[obj_index]+local_index];
       else
         return dof_indices[dof_offsets[obj_index]]+local_index;
@@ -336,7 +397,7 @@ 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,
+      Assert (is_compressed_entry(active_fe_indices[obj_index]) == false,
               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"));
@@ -353,10 +414,10 @@ namespace internal
       Assert (obj_index < active_fe_indices.size(),
               ExcIndexRange (obj_index, 0, active_fe_indices.size()));
 
-      if (((signed_active_fe_index_type)active_fe_indices[obj_index]) >= 0)
+      if (is_compressed_entry(active_fe_indices[obj_index]) == false)
         return active_fe_indices[obj_index];
       else
-        return (active_fe_index_type)~(signed_active_fe_index_type)active_fe_indices[obj_index];
+        return get_toggled_compression_state(active_fe_indices[obj_index]);
     }
 
 
@@ -367,10 +428,14 @@ namespace internal
     fe_index_is_active (const unsigned int                obj_index,
                         const unsigned int                fe_index) const
     {
-      return (fe_index == active_fe_index(obj_index));
+      if (is_compressed_entry(active_fe_indices[obj_index]) == false)
+        return (fe_index == active_fe_index(obj_index));
+      else
+        return (fe_index == get_toggled_compression_state(active_fe_index(obj_index)));
     }
 
 
+
     inline
     void
     DoFLevel::
@@ -380,6 +445,16 @@ namespace internal
       Assert (obj_index < active_fe_indices.size(),
               ExcIndexRange (obj_index, 0, active_fe_indices.size()));
 
+      // check whether the given fe_index is within the range of
+      // values that we interpret as "not compressed". if not, then
+      // the index is so large that we cannot accept it. (but this
+      // will not likely happen because it requires someone using an
+      // FECollection that has more than 32k entries.)
+      Assert (is_compressed_entry (fe_index) == false,
+              ExcMessage ("You are using an active_fe_index that is larger than an "
+                          "internal limitation for these objects. Try to work with "
+                          "hp::FECollection objects that have a more modest size."));
+
       active_fe_indices[obj_index] = fe_index;
     }
 
@@ -405,6 +480,8 @@ namespace internal
       return &cell_dof_indices_cache[cell_cache_offsets[obj_index]];
     }
 
+
+
     template <class Archive>
     inline
     void
index 18d7c11403658a550bd1b8d278c328a00b666601..65833da8bbc574df44db5b57a851724548e01e88 100644 (file)
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 2003 - 2016 by the deal.II authors
+// Copyright (C) 2003 - 2017 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -3175,6 +3175,14 @@ namespace hp
                     tria->n_raw_cells(level),
                     ExcInternalError ());
           }
+
+        // it may be that the previous table was compressed; in that
+        // case, restore the correct active_fe_index. the fact that
+        // this no longer matches the indices in the table is of no
+        // importance because the current function is called at a
+        // point where we have to recreate the dof_indices tables in
+        // the levels anyway
+        levels[level]->normalize_active_fe_indices ();
       }
   }
 
index 595d16f33d248d802fab9af69dde5a2cd3aa7adc..b8506563fe05ca577448190934bf88e26e5efc35 100644 (file)
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 2003 - 2015 by the deal.II authors
+// Copyright (C) 2003 - 2017 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -29,7 +29,6 @@ namespace internal
     DoFLevel::compress_data (const dealii::hp::FECollection<dim,spacedim> &fe_collection)
     {
       (void)fe_collection;
-      return;
 
       if (dof_offsets.size() == 0 || dof_indices.size()==0)
         return;
@@ -81,6 +80,7 @@ namespace internal
       // 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))
@@ -98,6 +98,8 @@ namespace internal
             Assert (next_offset-dof_offsets[cell] == fe_collection[active_fe_indices[cell]].template n_dofs_per_object<dim>(),
                     ExcInternalError());
 
+            new_dof_offsets[cell] = new_dof_indices.size();
+
             // 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])
@@ -118,10 +120,11 @@ namespace internal
 
                     // 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());
+                    Assert (is_compressed_entry(active_fe_indices[cell]) == false,
+                            ExcInternalError());
 
                     // then mark the compression
-                    active_fe_indices[cell] = (active_fe_index_type)~(signed_active_fe_index_type)active_fe_indices[cell];
+                    active_fe_indices[cell] = get_toggled_compression_state(active_fe_indices[cell]);
                   }
                 else
                   for (unsigned int i=dof_offsets[cell]; i<next_offset; ++i)
@@ -137,6 +140,7 @@ namespace internal
       // finally swap old and new content
       Assert (new_dof_indices.size() == new_size, ExcInternalError());
       dof_indices.swap (new_dof_indices);
+      dof_offsets.swap (new_dof_offsets);
     }
 
 
@@ -145,8 +149,6 @@ namespace internal
     void
     DoFLevel::uncompress_data(const dealii::hp::FECollection<dim,spacedim> &fe_collection)
     {
-      return;
-
       if (dof_offsets.size() == 0 || dof_indices.size()==0)
         return;
 
@@ -183,7 +185,7 @@ namespace internal
             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)
+            if (is_compressed_entry(active_fe_indices[cell]) == false)
               {
                 // apparently not. simply copy them
                 Assert (next_offset-dof_offsets[cell] == fe_collection[active_fe_indices[cell]].template n_dofs_per_object<dim>(),
@@ -196,8 +198,14 @@ namespace internal
                 // 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)
+                const unsigned int dofs_per_object
+                  = fe_collection[get_toggled_compression_state(active_fe_indices[cell])]
+                    .template n_dofs_per_object<dim>();
+                for (unsigned int i=0; i<dofs_per_object; ++i)
                   new_dof_indices.push_back (dof_indices[dof_offsets[cell]]+i);
+
+                // then mark the uncompression
+                active_fe_indices[cell] = get_toggled_compression_state(active_fe_indices[cell]);
               }
 
             // then move on to the next cell
@@ -213,6 +221,7 @@ namespace internal
     }
 
 
+
     std::size_t
     DoFLevel::memory_consumption () const
     {
@@ -224,6 +233,17 @@ namespace internal
     }
 
 
+
+    void DoFLevel::normalize_active_fe_indices ()
+    {
+      for (unsigned int i=0; i<active_fe_indices.size(); ++i)
+        if (is_compressed_entry(active_fe_indices[i]))
+          active_fe_indices[i]
+            = get_toggled_compression_state(active_fe_indices[i]);
+    }
+
+
+
     // explicit instantiations
     template void DoFLevel::compress_data(const dealii::hp::FECollection<1,1> &);
     template void DoFLevel::compress_data(const dealii::hp::FECollection<1,2> &);

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.