From f391a18092e2aea179abe63049bd2b890724b565 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 4 May 2006 03:51:29 +0000 Subject: [PATCH] Allow introspection into which fe indices are actually used on what objects. We'll need this for hanging node constraints git-svn-id: https://svn.dealii.org/trunk@13038 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/dofs/dof_accessor.h | 47 +- .../include/dofs/dof_accessor.templates.h | 38 ++ deal.II/deal.II/include/dofs/dof_levels.h | 68 +++ deal.II/deal.II/include/dofs/hp_dof_levels.h | 176 +++++- deal.II/deal.II/source/dofs/hp_dof_levels.cc | 561 ++++++++++++++++-- tests/hp/Makefile | 3 +- tests/hp/n_active_fe_indices.cc | 158 +++++ 7 files changed, 997 insertions(+), 54 deletions(-) create mode 100644 tests/hp/n_active_fe_indices.cc diff --git a/deal.II/deal.II/include/dofs/dof_accessor.h b/deal.II/deal.II/include/dofs/dof_accessor.h index deeb3fd039..282212579f 100644 --- a/deal.II/deal.II/include/dofs/dof_accessor.h +++ b/deal.II/deal.II/include/dofs/dof_accessor.h @@ -390,7 +390,52 @@ class DoFAccessor : public DoFObjectAccessor_Inheritance::set_dof_index (const unsigned int i, +template +inline +unsigned int +DoFAccessor::n_active_fe_indices () const +{ + // access the respective DoF. the + // last argument disambiguates + // between the functions of same + // name in the DoFLevel hierarchy, + // to make sure we get at lines, + // quads, or hexes + return this->dof_handler->levels[this->present_level] + ->template n_active_fe_indices (*this->dof_handler, + this->present_index, + internal::StructuralDimension()); +} + + + +template +inline +bool +DoFAccessor::fe_index_is_active (const unsigned int fe_index) const +{ + // access the respective DoF. the + // last argument disambiguates + // between the functions of same + // name in the DoFLevel hierarchy, + // to make sure we get at lines, + // quads, or hexes + return this->dof_handler->levels[this->present_level] + ->template fe_index_is_active (*this->dof_handler, + this->present_index, + fe_index, + internal::StructuralDimension()); +} + + /*------------------------- Functions: DoFObjectAccessor<1,dim> -----------------------*/ diff --git a/deal.II/deal.II/include/dofs/dof_levels.h b/deal.II/deal.II/include/dofs/dof_levels.h index 8ec89df36c..24cb1ad589 100644 --- a/deal.II/deal.II/include/dofs/dof_levels.h +++ b/deal.II/deal.II/include/dofs/dof_levels.h @@ -173,6 +173,46 @@ namespace internal */ std::vector cell_dof_indices_cache; + /** + * Return the value 1. The + * meaning of this function + * becomes clear by looking + * at what the corresponding + * functions in the classes + * internal::hp::DoFLevel do. + * Since we always have to + * return 1, there is no + * point in implementing this + * function multiple times in + * the derived classes, and + * we just do it once as a + * template in the common + * base class. + */ + template + unsigned int + n_active_fe_indices (const ::DoFHandler &dof_handler, + const unsigned int index, + const StructuralDimension) const; + + /** + * Similar to the function + * above. Assert that the + * given index is zero, and + * the return true. + * + * See the + * internal::hp::DoFLevel + * hierarchy for more + * information. + */ + template + bool + fe_index_is_active (const ::DoFHandler &dof_handler, + const unsigned int index, + const unsigned int fe_index, + const StructuralDimension) const; + /** * Determine an estimate for the * memory consumption (in bytes) @@ -564,6 +604,34 @@ namespace internal template friend class ::DoFHandler; template friend class ::MGDoFHandler; }; + + +// ----------- inline functions --------------------- + + template + inline + unsigned int + DoFLevel<0>::n_active_fe_indices (const ::DoFHandler &, + const unsigned int, + const StructuralDimension) const + { + return 1; + } + + + template + inline + bool + DoFLevel<0>::fe_index_is_active (const ::DoFHandler &, + const unsigned int, + const unsigned int fe_index, + const StructuralDimension) const + { + Assert (fe_index == 0, + ExcMessage ("Only zero fe_index values are allowed for " + "non-hp DoFHandlers.")); + return true; + } } diff --git a/deal.II/deal.II/include/dofs/hp_dof_levels.h b/deal.II/deal.II/include/dofs/hp_dof_levels.h index 8252aaf43c..c280692ab7 100644 --- a/deal.II/deal.II/include/dofs/hp_dof_levels.h +++ b/deal.II/deal.II/include/dofs/hp_dof_levels.h @@ -263,6 +263,62 @@ namespace internal const unsigned int local_index, internal::StructuralDimension<1> dummy) const; + /** + * Return the number of + * finite elements that are + * active on a given + * object. If this is a cell, + * the answer is of course + * one. If it is a face, the + * answer may be one or two, + * depending on whether the + * two adjacent cells use the + * same finite element or + * not. If it is an edge in + * 3d, the possible return + * value may be one or any + * other value larger than + * that. + * + * The last argument is used + * in the same way as for the + * other functions, i.e. it + * disambiguates between the + * different functions in the + * class hierarchy, whereas + * the template argument + * denotes the space + * dimension we operate in. + */ + template + unsigned int + n_active_fe_indices (const ::hp::DoFHandler &dof_handler, + const unsigned int line_index, + const StructuralDimension<1>) const; + + /** + * Check whether a given + * finite element index is + * used on the present + * object or not. + * + * The last argument is used + * in the same way as for the + * other functions, i.e. it + * disambiguates between the + * different functions in the + * class hierarchy, whereas + * the template argument + * denotes the space + * dimension we operate in. + */ + template + bool + fe_index_is_active (const ::hp::DoFHandler &dof_handler, + const unsigned int line_index, + const unsigned int fe_index, + const StructuralDimension<1>) const; + /** * Determine an estimate for the * memory consumption (in bytes) @@ -388,13 +444,71 @@ namespace internal const unsigned int local_index, internal::StructuralDimension<2> dummy) const; - /** + /** + * Return the number of + * finite elements that are + * active on a given + * object. If this is a cell, + * the answer is of course + * one. If it is a face, the + * answer may be one or two, + * depending on whether the + * two adjacent cells use the + * same finite element or + * not. If it is an edge in + * 3d, the possible return + * value may be one or any + * other value larger than + * that. + * + * The last argument is used + * in the same way as for the + * other functions, i.e. it + * disambiguates between the + * different functions in the + * class hierarchy, whereas + * the template argument + * denotes the space + * dimension we operate in. + */ + template + unsigned int + n_active_fe_indices (const ::hp::DoFHandler &dof_handler, + const unsigned int quad_index, + const StructuralDimension<2>) const; + + /** + * Check whether a given + * finite element index is + * used on the present + * object or not. + * + * The last argument is used + * in the same way as for the + * other functions, i.e. it + * disambiguates between the + * different functions in the + * class hierarchy, whereas + * the template argument + * denotes the space + * dimension we operate in. + */ + template + bool + fe_index_is_active (const ::hp::DoFHandler &dof_handler, + const unsigned int quad_index, + const unsigned int fe_index, + const StructuralDimension<2>) const; + + /** * Import the respective * functions from the base * class. */ using DoFLevel<1>::set_dof_index; using DoFLevel<1>::get_dof_index; + using DoFLevel<1>::n_active_fe_indices; + using DoFLevel<1>::fe_index_is_active; /** * Determine an estimate for the @@ -522,13 +636,71 @@ namespace internal const unsigned int local_index, internal::StructuralDimension<3> dummy) const; - /** + /** + * Return the number of + * finite elements that are + * active on a given + * object. If this is a cell, + * the answer is of course + * one. If it is a face, the + * answer may be one or two, + * depending on whether the + * two adjacent cells use the + * same finite element or + * not. If it is an edge in + * 3d, the possible return + * value may be one or any + * other value larger than + * that. + * + * The last argument is used + * in the same way as for the + * other functions, i.e. it + * disambiguates between the + * different functions in the + * class hierarchy, whereas + * the template argument + * denotes the space + * dimension we operate in. + */ + template + unsigned int + n_active_fe_indices (const ::hp::DoFHandler &dof_handler, + const unsigned int hex_index, + const StructuralDimension<3>) const; + + /** + * Check whether a given + * finite element index is + * used on the present + * object or not. + * + * The last argument is used + * in the same way as for the + * other functions, i.e. it + * disambiguates between the + * different functions in the + * class hierarchy, whereas + * the template argument + * denotes the space + * dimension we operate in. + */ + template + bool + fe_index_is_active (const ::hp::DoFHandler &dof_handler, + const unsigned int hex_index, + const unsigned int fe_index, + const StructuralDimension<3>) const; + + /** * Import the respective * functions from the base * classes. */ using DoFLevel<2>::set_dof_index; using DoFLevel<2>::get_dof_index; + using DoFLevel<2>::n_active_fe_indices; + using DoFLevel<2>::fe_index_is_active; /** * Determine an estimate for the diff --git a/deal.II/deal.II/source/dofs/hp_dof_levels.cc b/deal.II/deal.II/source/dofs/hp_dof_levels.cc index dc04c1b02b..cdd8d96b75 100644 --- a/deal.II/deal.II/source/dofs/hp_dof_levels.cc +++ b/deal.II/deal.II/source/dofs/hp_dof_levels.cc @@ -44,6 +44,8 @@ namespace internal Assert (local_index < dof_handler.get_fe()[fe_index].dofs_per_line, ExcIndexRange(local_index, 0, dof_handler.get_fe()[fe_index].dofs_per_line)); + Assert (line_index < line_dof_offsets.size(), + ExcIndexRange (line_index, 0, line_dof_offsets.size())); // make sure we are on an // object for which DoFs have @@ -100,11 +102,11 @@ namespace internal void DoFLevel<1>:: set_dof_index (const ::hp::DoFHandler &dof_handler, - const unsigned int line_index, - const unsigned int fe_index, - const unsigned int local_index, - const unsigned int global_index, - internal::StructuralDimension<1> /*dummy*/) + const unsigned int line_index, + const unsigned int fe_index, + const unsigned int local_index, + const unsigned int global_index, + internal::StructuralDimension<1> /*dummy*/) { Assert (fe_index != ::hp::DoFHandler::default_fe_index, ExcMessage ("You need to specify a FE index when working " @@ -117,6 +119,8 @@ namespace internal Assert (local_index < dof_handler.get_fe()[fe_index].dofs_per_line, ExcIndexRange(local_index, 0, dof_handler.get_fe()[fe_index].dofs_per_line)); + Assert (line_index < line_dof_offsets.size(), + ExcIndexRange (line_index, 0, line_dof_offsets.size())); // make sure we are on an // object for which DoFs have @@ -172,14 +176,149 @@ namespace internal + template + unsigned int + DoFLevel<1>:: + n_active_fe_indices (const ::hp::DoFHandler &dof_handler, + const unsigned int line_index, + internal::StructuralDimension<1> /*dummy*/) const + { + Assert (&dof_handler != 0, + ExcMessage ("No DoFHandler is specified for this iterator")); + Assert (&dof_handler.get_fe() != 0, + ExcMessage ("No finite element collection is associated with " + "this DoFHandler")); + Assert (line_index < line_dof_offsets.size(), + ExcIndexRange (line_index, 0, line_dof_offsets.size())); + + // make sure we are on an + // object for which DoFs have + // been allocated at all + Assert (line_dof_offsets[line_index] != deal_II_numbers::invalid_unsigned_int, + ExcMessage ("You are trying to access degree of freedom " + "information for an object on which no such " + "information is available")); + + // if we are in 1d, then the + // only set of indices we store + // is the one for the cell, + // which is unique + if (dim == 1) + { + return 1; + } + else + { + // we are in higher space + // dimensions, so there may + // be multiple finite + // elements associated with + // this object. hop along + // the list of index sets + // until we find the one + // with the correct + // fe_index, and then poke + // into that part. trigger + // an exception if we can't + // find a set for this + // particular fe_index + const unsigned int starting_offset = line_dof_offsets[line_index]; + const unsigned int *pointer = &line_dofs[starting_offset]; + unsigned int counter = 0; + while (true) + { + if (*pointer == deal_II_numbers::invalid_unsigned_int) + // end of list reached + return counter; + else + { + ++counter; + pointer += dof_handler.get_fe()[*pointer].dofs_per_line + 1; + } + } + } + } + + + + template + bool + DoFLevel<1>:: + fe_index_is_active (const ::hp::DoFHandler &dof_handler, + const unsigned int line_index, + const unsigned int fe_index, + internal::StructuralDimension<1> /*dummy*/) const + { + Assert (&dof_handler != 0, + ExcMessage ("No DoFHandler is specified for this iterator")); + Assert (&dof_handler.get_fe() != 0, + ExcMessage ("No finite element collection is associated with " + "this DoFHandler")); + Assert (line_index < line_dof_offsets.size(), + ExcIndexRange (line_index, 0, line_dof_offsets.size())); + Assert (fe_index != ::hp::DoFHandler::default_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 (line_dof_offsets[line_index] != deal_II_numbers::invalid_unsigned_int, + ExcMessage ("You are trying to access degree of freedom " + "information for an object on which no such " + "information is available")); + + // if we are in 1d, then the + // only set of indices we store + // is the one for the cell, + // which is unique + if (dim == 1) + { + Assert (line_index < active_fe_indices.size(), + ExcInternalError()); + return (fe_index == active_fe_indices[line_index]); + } + else + { + // we are in higher space + // dimensions, so there may + // be multiple finite + // elements associated with + // this object. hop along + // the list of index sets + // until we find the one + // with the correct + // fe_index, and then poke + // into that part. trigger + // an exception if we can't + // find a set for this + // particular fe_index + const unsigned int starting_offset = line_dof_offsets[line_index]; + const unsigned int *pointer = &line_dofs[starting_offset]; + while (true) + { + if (*pointer == deal_II_numbers::invalid_unsigned_int) + // end of list reached + return false; + else + if (*pointer == fe_index) + return true; + else + pointer += dof_handler.get_fe()[*pointer].dofs_per_line + 1; + } + } + } + + + template unsigned int DoFLevel<2>:: get_dof_index (const ::hp::DoFHandler &dof_handler, - const unsigned int quad_index, - const unsigned int fe_index, - const unsigned int local_index, - internal::StructuralDimension<2> /*dummy*/) const + const unsigned int quad_index, + const unsigned int fe_index, + const unsigned int local_index, + internal::StructuralDimension<2> /*dummy*/) const { Assert (dim >= 2, ExcMessage ("You can only access quads in 2d or higher")); Assert (fe_index != ::hp::DoFHandler::default_fe_index, @@ -193,6 +332,8 @@ namespace internal Assert (local_index < dof_handler.get_fe()[fe_index].dofs_per_quad, ExcIndexRange(local_index, 0, dof_handler.get_fe()[fe_index].dofs_per_quad)); + Assert (quad_index < quad_dof_offsets.size(), + ExcIndexRange (quad_index, 0, quad_dof_offsets.size())); // make sure we are on an // object for which DoFs have @@ -249,11 +390,11 @@ namespace internal void DoFLevel<2>:: set_dof_index (const ::hp::DoFHandler &dof_handler, - const unsigned int quad_index, - const unsigned int fe_index, - const unsigned int local_index, - const unsigned int global_index, - internal::StructuralDimension<2> /*dummy*/) + const unsigned int quad_index, + const unsigned int fe_index, + const unsigned int local_index, + const unsigned int global_index, + internal::StructuralDimension<2> /*dummy*/) { Assert (dim >= 2, ExcMessage ("You can only access quads in 2d or higher")); Assert (fe_index != ::hp::DoFHandler::default_fe_index, @@ -267,6 +408,8 @@ namespace internal Assert (local_index < dof_handler.get_fe()[fe_index].dofs_per_quad, ExcIndexRange(local_index, 0, dof_handler.get_fe()[fe_index].dofs_per_quad)); + Assert (quad_index < quad_dof_offsets.size(), + ExcIndexRange (quad_index, 0, quad_dof_offsets.size())); // make sure we are on an // object for which DoFs have @@ -322,14 +465,149 @@ namespace internal + template + unsigned int + DoFLevel<2>:: + n_active_fe_indices (const ::hp::DoFHandler &dof_handler, + const unsigned int quad_index, + internal::StructuralDimension<2> /*dummy*/) const + { + Assert (&dof_handler != 0, + ExcMessage ("No DoFHandler is specified for this iterator")); + Assert (&dof_handler.get_fe() != 0, + ExcMessage ("No finite element collection is associated with " + "this DoFHandler")); + Assert (quad_index < quad_dof_offsets.size(), + ExcIndexRange (quad_index, 0, quad_dof_offsets.size())); + + // make sure we are on an + // object for which DoFs have + // been allocated at all + Assert (quad_dof_offsets[quad_index] != deal_II_numbers::invalid_unsigned_int, + ExcMessage ("You are trying to access degree of freedom " + "information for an object on which no such " + "information is available")); + + // if we are in 2d, then the + // only set of indices we store + // is the one for the cell, + // which is unique + if (dim == 2) + { + return 1; + } + else + { + // we are in higher space + // dimensions, so there may + // be multiple finite + // elements associated with + // this object. hop along + // the list of index sets + // until we find the one + // with the correct + // fe_index, and then poke + // into that part. trigger + // an exception if we can't + // find a set for this + // particular fe_index + const unsigned int starting_offset = quad_dof_offsets[quad_index]; + const unsigned int *pointer = &quad_dofs[starting_offset]; + unsigned int counter = 0; + while (true) + { + if (*pointer == deal_II_numbers::invalid_unsigned_int) + // end of list reached + return counter; + else + { + ++counter; + pointer += dof_handler.get_fe()[*pointer].dofs_per_quad + 1; + } + } + } + } + + + + template + bool + DoFLevel<2>:: + fe_index_is_active (const ::hp::DoFHandler &dof_handler, + const unsigned int quad_index, + const unsigned int fe_index, + internal::StructuralDimension<2> /*dummy*/) const + { + Assert (&dof_handler != 0, + ExcMessage ("No DoFHandler is specified for this iterator")); + Assert (&dof_handler.get_fe() != 0, + ExcMessage ("No finite element collection is associated with " + "this DoFHandler")); + Assert (quad_index < quad_dof_offsets.size(), + ExcIndexRange (quad_index, 0, quad_dof_offsets.size())); + Assert (fe_index != ::hp::DoFHandler::default_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 (quad_dof_offsets[quad_index] != deal_II_numbers::invalid_unsigned_int, + ExcMessage ("You are trying to access degree of freedom " + "information for an object on which no such " + "information is available")); + + // if we are in 2d, then the + // only set of indices we store + // is the one for the cell, + // which is unique + if (dim == 2) + { + Assert (quad_index < active_fe_indices.size(), + ExcInternalError()); + return (fe_index == active_fe_indices[quad_index]); + } + else + { + // we are in higher space + // dimensions, so there may + // be multiple finite + // elements associated with + // this object. hop along + // the list of index sets + // until we find the one + // with the correct + // fe_index, and then poke + // into that part. trigger + // an exception if we can't + // find a set for this + // particular fe_index + const unsigned int starting_offset = quad_dof_offsets[quad_index]; + const unsigned int *pointer = &quad_dofs[starting_offset]; + while (true) + { + if (*pointer == deal_II_numbers::invalid_unsigned_int) + // end of list reached + return false; + else + if (*pointer == fe_index) + return true; + else + pointer += dof_handler.get_fe()[*pointer].dofs_per_quad + 1; + } + } + } + + + template unsigned int DoFLevel<3>:: get_dof_index (const ::hp::DoFHandler &dof_handler, - const unsigned int hex_index, - const unsigned int fe_index, - const unsigned int local_index, - internal::StructuralDimension<3> /*dummy*/) const + const unsigned int hex_index, + const unsigned int fe_index, + const unsigned int local_index, + internal::StructuralDimension<3> /*dummy*/) const { Assert (dim >= 3, ExcMessage ("You can only access hexs in 3d or higher")); Assert (fe_index != ::hp::DoFHandler::default_fe_index, @@ -343,6 +621,8 @@ namespace internal Assert (local_index < dof_handler.get_fe()[fe_index].dofs_per_hex, ExcIndexRange(local_index, 0, dof_handler.get_fe()[fe_index].dofs_per_hex)); + Assert (hex_index < hex_dof_offsets.size(), + ExcIndexRange (hex_index, 0, hex_dof_offsets.size())); // make sure we are on an // object for which DoFs have @@ -399,11 +679,11 @@ namespace internal void DoFLevel<3>:: set_dof_index (const ::hp::DoFHandler &dof_handler, - const unsigned int hex_index, - const unsigned int fe_index, - const unsigned int local_index, - const unsigned int global_index, - internal::StructuralDimension<3> /*dummy*/) + const unsigned int hex_index, + const unsigned int fe_index, + const unsigned int local_index, + const unsigned int global_index, + internal::StructuralDimension<3> /*dummy*/) { Assert (dim >= 3, ExcMessage ("You can only access hexs in 3d or higher")); Assert (fe_index != ::hp::DoFHandler::default_fe_index, @@ -417,6 +697,8 @@ namespace internal Assert (local_index < dof_handler.get_fe()[fe_index].dofs_per_hex, ExcIndexRange(local_index, 0, dof_handler.get_fe()[fe_index].dofs_per_hex)); + Assert (hex_index < hex_dof_offsets.size(), + ExcIndexRange (hex_index, 0, hex_dof_offsets.size())); // make sure we are on an // object for which DoFs have @@ -469,7 +751,142 @@ namespace internal } } } + + + + template + unsigned int + DoFLevel<3>:: + n_active_fe_indices (const ::hp::DoFHandler &dof_handler, + const unsigned int hex_index, + internal::StructuralDimension<3> /*dummy*/) const + { + Assert (&dof_handler != 0, + ExcMessage ("No DoFHandler is specified for this iterator")); + Assert (&dof_handler.get_fe() != 0, + ExcMessage ("No finite element collection is associated with " + "this DoFHandler")); + Assert (hex_index < hex_dof_offsets.size(), + ExcIndexRange (hex_index, 0, hex_dof_offsets.size())); + + // make sure we are on an + // object for which DoFs have + // been allocated at all + Assert (hex_dof_offsets[hex_index] != deal_II_numbers::invalid_unsigned_int, + ExcMessage ("You are trying to access degree of freedom " + "information for an object on which no such " + "information is available")); + + // if we are in 3d, then the + // only set of indices we store + // is the one for the cell, + // which is unique + if (dim == 3) + { + return 1; + } + else + { + // we are in higher space + // dimensions, so there may + // be multiple finite + // elements associated with + // this object. hop along + // the list of index sets + // until we find the one + // with the correct + // fe_index, and then poke + // into that part. trigger + // an exception if we can't + // find a set for this + // particular fe_index + const unsigned int starting_offset = hex_dof_offsets[hex_index]; + const unsigned int *pointer = &hex_dofs[starting_offset]; + unsigned int counter = 0; + while (true) + { + if (*pointer == deal_II_numbers::invalid_unsigned_int) + // end of list reached + return counter; + else + { + ++counter; + pointer += dof_handler.get_fe()[*pointer].dofs_per_hex + 1; + } + } + } + } + + + template + bool + DoFLevel<3>:: + fe_index_is_active (const ::hp::DoFHandler &dof_handler, + const unsigned int hex_index, + const unsigned int fe_index, + internal::StructuralDimension<3> /*dummy*/) const + { + Assert (&dof_handler != 0, + ExcMessage ("No DoFHandler is specified for this iterator")); + Assert (&dof_handler.get_fe() != 0, + ExcMessage ("No finite element collection is associated with " + "this DoFHandler")); + Assert (hex_index < hex_dof_offsets.size(), + ExcIndexRange (hex_index, 0, hex_dof_offsets.size())); + Assert (fe_index != ::hp::DoFHandler::default_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 (hex_dof_offsets[hex_index] != deal_II_numbers::invalid_unsigned_int, + ExcMessage ("You are trying to access degree of freedom " + "information for an object on which no such " + "information is available")); + + // if we are in 2d, then the + // only set of indices we store + // is the one for the cell, + // which is unique + if (dim == 3) + { + Assert (hex_index < active_fe_indices.size(), + ExcInternalError()); + return (fe_index == active_fe_indices[hex_index]); + } + else + { + // we are in higher space + // dimensions, so there may + // be multiple finite + // elements associated with + // this object. hop along + // the list of index sets + // until we find the one + // with the correct + // fe_index, and then poke + // into that part. trigger + // an exception if we can't + // find a set for this + // particular fe_index + const unsigned int starting_offset = hex_dof_offsets[hex_index]; + const unsigned int *pointer = &hex_dofs[starting_offset]; + while (true) + { + if (*pointer == deal_II_numbers::invalid_unsigned_int) + // end of list reached + return false; + else + if (*pointer == fe_index) + return true; + else + pointer += dof_handler.get_fe()[*pointer].dofs_per_hex + 1; + } + } + } + // explicit instantiations @@ -477,42 +894,72 @@ namespace internal unsigned int DoFLevel<1>:: get_dof_index (const ::hp::DoFHandler &dof_handler, - const unsigned int line_index, - const unsigned int fe_index, - const unsigned int local_index, - internal::StructuralDimension<1> dummy) const; + const unsigned int line_index, + const unsigned int fe_index, + const unsigned int local_index, + internal::StructuralDimension<1> dummy) const; template void DoFLevel<1>:: set_dof_index (const ::hp::DoFHandler &dof_handler, - const unsigned int line_index, - const unsigned int fe_index, - const unsigned int local_index, - const unsigned int global_index, - internal::StructuralDimension<1> dummy); + const unsigned int line_index, + const unsigned int fe_index, + const unsigned int local_index, + const unsigned int global_index, + internal::StructuralDimension<1> dummy); + template + unsigned int + DoFLevel<1>:: + n_active_fe_indices (const ::hp::DoFHandler &dof_handler, + const unsigned int line_index, + internal::StructuralDimension<1> dummy) const; + + template + bool + DoFLevel<1>:: + fe_index_is_active (const ::hp::DoFHandler &dof_handler, + const unsigned int line_index, + const unsigned int fe_index, + internal::StructuralDimension<1> dummy) const; + #if deal_II_dimension >= 2 template unsigned int DoFLevel<2>:: get_dof_index (const ::hp::DoFHandler &dof_handler, - const unsigned int quad_index, - const unsigned int fe_index, - const unsigned int local_index, - internal::StructuralDimension<2> dummy) const; + const unsigned int quad_index, + const unsigned int fe_index, + const unsigned int local_index, + internal::StructuralDimension<2> dummy) const; template void DoFLevel<2>:: set_dof_index (const ::hp::DoFHandler &dof_handler, - const unsigned int quad_index, - const unsigned int fe_index, - const unsigned int local_index, - const unsigned int global_index, - internal::StructuralDimension<2> dummy); + const unsigned int quad_index, + const unsigned int fe_index, + const unsigned int local_index, + const unsigned int global_index, + internal::StructuralDimension<2> dummy); + template + unsigned int + DoFLevel<2>:: + n_active_fe_indices (const ::hp::DoFHandler &dof_handler, + const unsigned int quad_index, + internal::StructuralDimension<2> dummy) const; + + + template + bool + DoFLevel<2>:: + fe_index_is_active (const ::hp::DoFHandler &dof_handler, + const unsigned int line_index, + const unsigned int fe_index, + internal::StructuralDimension<2> dummy) const; #endif #if deal_II_dimension >= 3 @@ -521,21 +968,35 @@ namespace internal unsigned int DoFLevel<3>:: get_dof_index (const ::hp::DoFHandler &dof_handler, - const unsigned int hex_index, - const unsigned int fe_index, - const unsigned int local_index, - internal::StructuralDimension<3> dummy) const; + const unsigned int hex_index, + const unsigned int fe_index, + const unsigned int local_index, + internal::StructuralDimension<3> dummy) const; template void DoFLevel<3>:: set_dof_index (const ::hp::DoFHandler &dof_handler, - const unsigned int hex_index, - const unsigned int fe_index, - const unsigned int local_index, - const unsigned int global_index, - internal::StructuralDimension<3> dummy); + const unsigned int hex_index, + const unsigned int fe_index, + const unsigned int local_index, + const unsigned int global_index, + internal::StructuralDimension<3> dummy); + template + unsigned int + DoFLevel<3>:: + n_active_fe_indices (const ::hp::DoFHandler &dof_handler, + const unsigned int hex_index, + internal::StructuralDimension<3> dummy) const; + + template + bool + DoFLevel<3>:: + fe_index_is_active (const ::hp::DoFHandler &dof_handler, + const unsigned int line_index, + const unsigned int fe_index, + internal::StructuralDimension<3> dummy) const; #endif } } diff --git a/tests/hp/Makefile b/tests/hp/Makefile index 6ef799bf48..59aee90e20 100644 --- a/tests/hp/Makefile +++ b/tests/hp/Makefile @@ -29,7 +29,8 @@ tests_x = hp_dof_handler \ continuous_* \ n_dofs \ random \ - step-* + step-* \ + n_active_fe_indices # from above list of regular expressions, generate the real set of # tests diff --git a/tests/hp/n_active_fe_indices.cc b/tests/hp/n_active_fe_indices.cc new file mode 100644 index 0000000000..9b3085d720 --- /dev/null +++ b/tests/hp/n_active_fe_indices.cc @@ -0,0 +1,158 @@ +//---------------------------- n_active_fe_indices.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005, 2006 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- n_active_fe_indices.cc --------------------------- + + +// record the output of n_active_fe_indices for cells, faces, and +// edges, as well as fe_index_is_active + + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +template +void check_cells (const hp::DoFHandler &dof_handler) +{ + for (typename hp::DoFHandler::active_cell_iterator + cell=dof_handler.begin_active(); + cell!=dof_handler.end(); ++cell) + { + deallog << "cell=" << cell << std::endl; + deallog << "n=" + << cell->n_active_fe_indices () + << std::endl; + deallog << "x="; + for (unsigned int i=0; ife_index_is_active(i); + deallog << std::endl; + + Assert (cell->n_active_fe_indices () == 1, ExcInternalError()); + } +} + + +void check_faces (const hp::DoFHandler<1> &) +{} + + +template +void check_faces (const hp::DoFHandler &dof_handler) +{ + for (typename hp::DoFHandler::active_cell_iterator + cell=dof_handler.begin_active(); + cell!=dof_handler.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + { + deallog << "face=" << cell->face(f) << std::endl; + deallog << "n=" + << cell->face(f)->n_active_fe_indices () + << std::endl; + deallog << "x="; + for (unsigned int i=0; iface(f)->fe_index_is_active(i); + deallog << std::endl; + + Assert (cell->face(f)->n_active_fe_indices () >= 1, ExcInternalError()); + Assert (cell->face(f)->n_active_fe_indices () <= 2, ExcInternalError()); + } +} + + +void check_edges (const hp::DoFHandler<1> &) +{} + + +void check_edges (const hp::DoFHandler<2> &) +{} + + +template +void check_edges (const hp::DoFHandler &dof_handler) +{ + for (typename hp::DoFHandler::active_cell_iterator + cell=dof_handler.begin_active(); + cell!=dof_handler.end(); ++cell) + for (unsigned int e=0; e::lines_per_cell; ++e) + { + deallog << "edge=" << cell->line(e) << std::endl; + deallog << "n=" + << cell->line(e)->n_active_fe_indices () + << std::endl; + deallog << "x="; + for (unsigned int i=0; iline(e)->fe_index_is_active(i); + deallog << std::endl; + + Assert (cell->line(e)->n_active_fe_indices () >= 1, ExcInternalError()); + Assert (cell->line(e)->n_active_fe_indices () <= dof_handler.get_fe().size(), + ExcInternalError()); + } +} + + + +template +void test () +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global (2); + tria.begin_active()->set_refine_flag (); + tria.execute_coarsening_and_refinement (); + tria.refine_global (2); + + hp::FECollection fe_collection; + fe_collection.push_back(FE_Q (1)); + fe_collection.push_back(FE_Q (2)); + fe_collection.push_back(FE_Q (3)); + fe_collection.push_back(FE_Q (4)); + + hp::DoFHandler dof_handler(tria); + + for (typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + cell->set_active_fe_index (rand() % fe_collection.size()); + + dof_handler.distribute_dofs(fe_collection); + + check_cells (dof_handler); + check_faces (dof_handler); + check_edges (dof_handler); +} + + +int main () +{ + std::ofstream logfile("n_active_fe_indices/output"); + logfile.precision(2); + + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + test<1> (); + test<2> (); + test<3> (); + + deallog << "OK" << std::endl; +} -- 2.39.5