From: Wolfgang Bangerth Date: Mon, 9 Sep 2013 19:10:18 +0000 (+0000) Subject: Implement the compression scheme for DoF indices stored on cells. X-Git-Tag: v8.1.0~847 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=80fd6fb73f81e8021a37003b20b9522bd17a320a;p=dealii.git Implement the compression scheme for DoF indices stored on cells. git-svn-id: https://svn.dealii.org/trunk@30655 0785d39b-7218-0410-832d-ea1e28bc413d --- 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 class DoFLevel; + class DoFLevel; template class DoFFaces; template class DoFObjects; @@ -790,7 +790,7 @@ namespace hp * levels[] tree of * the Triangulation objects. */ - std::vector*> levels; + std::vector levels; /** * Space to store the DoF @@ -900,7 +900,6 @@ namespace hp * the functions that set and * retrieve vertex dof indices. */ - template friend class dealii::internal::hp::DoFLevel; template 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 class DoFHandler; + template 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 interior 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. + * + *

Compression

+ * + * 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 12, 13, 14, 15. 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 12, 14, 13, 15, which + * can not be compressed using run-length encoding. */ - template 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 active_fe_indices; + std::vector 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::size_type> dof_offsets; + std::vector 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 + void compress_data (const dealii::hp::FECollection &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 + void uncompress_data (const dealii::hp::FECollection &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 inline types::global_dof_index - DoFLevel:: + 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 inline void - DoFLevel:: + 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 inline unsigned int - DoFLevel:: + 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 inline bool - DoFLevel:: + 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 inline void - DoFLevel:: + 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::size_type> dof_offsets; + std::vector 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 > active_fe_backup(dof_handler.levels.size ()); for (unsigned int level = 0; levelactive_fe_indices, @@ -553,7 +553,7 @@ namespace internal for (unsigned int level=0; leveln_levels(); ++level) { - dof_handler.levels.push_back (new internal::hp::DoFLevel); + 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; leveln_levels(); ++level) { dof_handler.levels[level]->dof_offsets - = std::vector::size_type> ( + = std::vector ( dof_handler.tria->n_raw_lines(level), - (std::vector::size_type)(-1)); + (DoFLevel::offset_type)(-1)); types::global_dof_index next_free_dof = 0; for (typename DoFHandler::active_cell_iterator @@ -619,7 +619,7 @@ namespace internal Assert (static_cast (std::count (dof_handler.levels[level]->dof_offsets.begin(), dof_handler.levels[level]->dof_offsets.end(), - (std::vector::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 > active_fe_backup(dof_handler.levels.size ()); for (unsigned int level = 0; levelactive_fe_indices, @@ -668,7 +668,7 @@ namespace internal for (unsigned int level=0; leveln_levels(); ++level) { - dof_handler.levels.push_back (new internal::hp::DoFLevel); + 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; leveln_levels(); ++level) { dof_handler.levels[level]->dof_offsets - = std::vector::size_type> ( + = std::vector ( dof_handler.tria->n_raw_quads(level), - (std::vector::size_type)(-1)); + (DoFLevel::offset_type)(-1)); types::global_dof_index next_free_dof = 0; for (typename DoFHandler::active_cell_iterator @@ -736,7 +736,7 @@ namespace internal Assert (static_cast (std::count (dof_handler.levels[level]->dof_offsets.begin(), dof_handler.levels[level]->dof_offsets.end(), - (std::vector::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::size_type> - (dof_handler.tria->n_raw_lines(), - (std::vector::size_type)(-1)); + = std::vector (dof_handler.tria->n_raw_lines(), + (unsigned int)(-1)); dof_handler.faces->lines.dofs = std::vector (n_line_slots, DoFHandler::invalid_dof_index); @@ -1020,7 +1019,7 @@ namespace internal // active_fe_indices field which // we have to backup before { - std::vector > + std::vector > active_fe_backup(dof_handler.levels.size ()); for (unsigned int level = 0; levelactive_fe_indices, @@ -1034,7 +1033,7 @@ namespace internal for (unsigned int level=0; leveln_levels(); ++level) { - dof_handler.levels.push_back (new internal::hp::DoFLevel); + 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; leveln_levels(); ++level) { dof_handler.levels[level]->dof_offsets - = std::vector::size_type> + = std::vector (dof_handler.tria->n_raw_hexs(level), - (std::vector::size_type)(-1)); + (DoFLevel::offset_type)(-1)); types::global_dof_index next_free_dof = 0; for (typename DoFHandler::active_cell_iterator @@ -1102,7 +1101,7 @@ namespace internal Assert (static_cast (std::count (dof_handler.levels[level]->dof_offsets.begin(), dof_handler.levels[level]->dof_offsets.end(), - (std::vector::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::size_type> + = std::vector (dof_handler.tria->n_raw_quads(), - (std::vector::size_type)(-1)); + (unsigned int)(-1)); dof_handler.faces->quads.dofs = std::vector (n_quad_slots, DoFHandler::invalid_dof_index); @@ -2691,6 +2690,9 @@ namespace hp = std::vector (1, number_cache.locally_owned_dofs); + for (unsigned int level=0; levelcompress_data (*finite_elements); + // finally restore the user flags const_cast &>(*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; leveluncompress_data (*finite_elements); + + // do the renumbering renumber_dofs_internal (new_numbers, dealii::internal::int2type()); + + // now re-compress the dof indices + for (unsigned int level=0; levelcompress_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); + 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); + 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 #include +#include +#include DEAL_II_NAMESPACE_OPEN @@ -23,9 +25,196 @@ namespace internal { namespace hp { - template + template + void + DoFLevel::compress_data (const dealii::hp::FECollection &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(), + 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 new_dof_indices; + new_dof_indices.reserve(new_size); + for (unsigned int cell=0; cell(), + 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= 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 + void + DoFLevel::uncompress_data(const dealii::hp::FECollection &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(); + } + + // 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=0) + { + // apparently not. simply copy them + Assert (next_offset-dof_offsets[cell] == fe_collection[active_fe_indices[cell]].template n_dofs_per_object(), + ExcInternalError()); + for (unsigned int i=dof_offsets[cell]; i(); ++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::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> &); } }