From 8114b52716247069e0a34ccf4ba2091bda1ca123 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 6 Jan 2017 08:14:03 -0700 Subject: [PATCH] Fix and enable compression of hp DoF indices. 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 | 97 ++++++++++++++++++++++++++++++---- source/hp/dof_handler.cc | 10 +++- source/hp/dof_level.cc | 36 ++++++++++--- 3 files changed, 124 insertions(+), 19 deletions(-) diff --git a/include/deal.II/hp/dof_level.h b/include/deal.II/hp/dof_level.h index 8c6316f7a3..0904bf58eb 100644 --- a/include/deal.II/hp/dof_level.h +++ b/include/deal.II/hp/dof_level.h @@ -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 void uncompress_data (const dealii::hp::FECollection &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 inline void diff --git a/source/hp/dof_handler.cc b/source/hp/dof_handler.cc index 18d7c11403..65833da8bb 100644 --- a/source/hp/dof_handler.cc +++ b/source/hp/dof_handler.cc @@ -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 (); } } diff --git a/source/hp/dof_level.cc b/source/hp/dof_level.cc index 595d16f33d..b8506563fe 100644 --- a/source/hp/dof_level.cc +++ b/source/hp/dof_level.cc @@ -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 &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 new_dof_indices; new_dof_indices.reserve(new_size); + std::vector new_dof_offsets (dof_offsets.size(), (offset_type)(-1)); for (unsigned int cell=0; cell(), 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 &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(), @@ -196,8 +198,14 @@ namespace internal // apparently so. uncompress Assert (next_offset-dof_offsets[cell] == 1, ExcInternalError()); - for (unsigned int i=0; i(); ++i) + const unsigned int dofs_per_object + = fe_collection[get_toggled_compression_state(active_fe_indices[cell])] + .template n_dofs_per_object(); + for (unsigned int i=0; i &); template void DoFLevel::compress_data(const dealii::hp::FECollection<1,2> &); -- 2.39.5