Allow introspection into which fe indices are actually used on what objects. We'll...

git-svn-id: https://svn.dealii.org/trunk@13038 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/deal.II/include/dofs/dof_accessor.h b/deal.II/deal.II/include/dofs/dof_accessor.h
index deeb3fd03926dbfd0095fa05abde361344ca574f..282212579f85f6eaf33ea6e6050ec7714a7f3fa0 100644 (file)
--- 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<structdim, DH::dimensio
                        const unsigned int index,
                        const unsigned int fe_index = DH::default_fe_index) const;
 
-                                    /**
+                                     /**
+                                      * Return the number of finite
+                                      * elements that are active on a
+                                      * given object.
+                                      *
+                                      * For non-hp DoFHandler objects,
+                                      * the answer is of course always
+                                      * one. However, for
+                                      * hp::DoFHandler objects, this
+                                      * isn't the case: 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.
+                                      */
+    unsigned int
+    n_active_fe_indices () const;
+
+                                     /**
+                                      * Check whether a given finite
+                                      * element index is used on the
+                                      * present object or not.
+                                      *
+                                      * For non-hp DoFHandler objects,
+                                      * only a zero index may be used,
+                                      * since only a single finite
+                                      * element can ever be active,
+                                      * and the answer is then @p
+                                      * true. On the other hand, for
+                                      * hp::DoFHandler objects,
+                                      * several fe indices may be used
+                                      * on the same object (except for
+                                      * a cell, where only a single
+                                      * finite element can be active
+                                      * at any given time).
+                                      */
+    bool
+    fe_index_is_active (const unsigned int fe_index) const;
+
+
+                                     /**
                                      * Exceptions for child classes
                                      *
                                      * @ingroup Exceptions

diff --git a/deal.II/deal.II/include/dofs/dof_accessor.templates.h b/deal.II/deal.II/include/dofs/dof_accessor.templates.h
index 6ed5d2e6b904a3fabcea280025172c63e7939ed9..6086d32c87105eccbaa689c5f4f0b07f1de76a27 100644 (file)
--- a/deal.II/deal.II/include/dofs/dof_accessor.templates.h
+++ b/deal.II/deal.II/include/dofs/dof_accessor.templates.h
@@ -223,6 +223,44 @@ DoFAccessor<dim,DH>::set_dof_index (const unsigned int i,
 
 
 
+template <int dim, class DH>
+inline
+unsigned int
+DoFAccessor<dim,DH>::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<DH::dimension> (*this->dof_handler,
+                                                   this->present_index,
+                                                   internal::StructuralDimension<dim>());
+}
+
+
+
+template <int dim, class DH>
+inline
+bool
+DoFAccessor<dim,DH>::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<DH::dimension> (*this->dof_handler,
+                                                  this->present_index,
+                                                  fe_index,
+                                                  internal::StructuralDimension<dim>());
+}
+
+
 
 /*------------------------- 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 8ec89df36c762811b9e497e6e864f2659d59b1be..24cb1ad589c7558e9c44189a7033eb3c39699a1a 100644 (file)
--- 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<unsigned int> 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 <int dim, int structdim>
+        unsigned int
+        n_active_fe_indices (const ::DoFHandler<dim> &dof_handler,
+                             const unsigned int       index,
+                             const StructuralDimension<structdim>) 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 <int dim, int structdim>
+        bool
+        fe_index_is_active (const ::DoFHandler<dim> &dof_handler,
+                            const unsigned int       index,
+                            const unsigned int       fe_index,
+                            const StructuralDimension<structdim>) const;
+        
                                          /**
                                           * Determine an estimate for the
                                           * memory consumption (in bytes)
@@ -564,6 +604,34 @@ namespace internal
         template <int> friend class ::DoFHandler;
         template <int> friend class ::MGDoFHandler;
     };
+
+
+// ----------- inline functions ---------------------
+
+    template <int dim, int structdim>
+    inline
+    unsigned int
+    DoFLevel<0>::n_active_fe_indices (const ::DoFHandler<dim> &,
+                                      const unsigned int,
+                                      const StructuralDimension<structdim>) const
+    {
+      return 1;
+    }
+
+
+    template <int dim, int structdim>
+    inline
+    bool
+    DoFLevel<0>::fe_index_is_active (const ::DoFHandler<dim> &,
+                                     const unsigned int,
+                                     const unsigned int fe_index,
+                                     const StructuralDimension<structdim>) 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 8252aaf43c492bf3afceaa804a88d28b2d6421b1..c280692ab77d1e0636dedcd2091680e34a257622 100644 (file)
--- 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 <int dim>
+        unsigned int
+        n_active_fe_indices (const ::hp::DoFHandler<dim> &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 <int dim>
+        bool
+        fe_index_is_active (const ::hp::DoFHandler<dim> &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 <int dim>
+        unsigned int
+        n_active_fe_indices (const ::hp::DoFHandler<dim> &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 <int dim>
+        bool
+        fe_index_is_active (const ::hp::DoFHandler<dim> &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 <int dim>
+        unsigned int
+        n_active_fe_indices (const ::hp::DoFHandler<dim> &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 <int dim>
+        bool
+        fe_index_is_active (const ::hp::DoFHandler<dim> &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 dc04c1b02b899ef193f4349d306f8428f61575fe..cdd8d96b75fc9b237331c8970db8aa33550dfe9d 100644 (file)
--- 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<dim> &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<dim>::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 <int dim>
+    unsigned int
+    DoFLevel<1>::
+    n_active_fe_indices (const ::hp::DoFHandler<dim> &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 <int dim>
+    bool
+    DoFLevel<1>::
+    fe_index_is_active (const ::hp::DoFHandler<dim> &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<dim>::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 <int dim>
     unsigned int
     DoFLevel<2>::
     get_dof_index (const ::hp::DoFHandler<dim> &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<dim>::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<dim> &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<dim>::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 <int dim>
+    unsigned int
+    DoFLevel<2>::
+    n_active_fe_indices (const ::hp::DoFHandler<dim> &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 <int dim>
+    bool
+    DoFLevel<2>::
+    fe_index_is_active (const ::hp::DoFHandler<dim> &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<dim>::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 <int dim>
     unsigned int
     DoFLevel<3>::
     get_dof_index (const ::hp::DoFHandler<dim> &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<dim>::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<dim> &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<dim>::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 <int dim>
+    unsigned int
+    DoFLevel<3>::
+    n_active_fe_indices (const ::hp::DoFHandler<dim> &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 <int dim>
+    bool
+    DoFLevel<3>::
+    fe_index_is_active (const ::hp::DoFHandler<dim> &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<dim>::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<deal_II_dimension> &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<deal_II_dimension> &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<deal_II_dimension> &dof_handler,
+                         const unsigned int           line_index,
+                         internal::StructuralDimension<1> dummy) const;
+
+    template
+    bool
+    DoFLevel<1>::
+    fe_index_is_active (const ::hp::DoFHandler<deal_II_dimension> &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<deal_II_dimension> &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<deal_II_dimension> &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<deal_II_dimension> &dof_handler,
+                         const unsigned int           quad_index,
+                         internal::StructuralDimension<2> dummy) const;
+    
+
+    template
+    bool
+    DoFLevel<2>::
+    fe_index_is_active (const ::hp::DoFHandler<deal_II_dimension> &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<deal_II_dimension> &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<deal_II_dimension> &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<deal_II_dimension> &dof_handler,
+                         const unsigned int           hex_index,
+                         internal::StructuralDimension<3> dummy) const;
+
+    template
+    bool
+    DoFLevel<3>::
+    fe_index_is_active (const ::hp::DoFHandler<deal_II_dimension> &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 6ef799bf48f18a72a9e8b96799c8ae9a41db5575..59aee90e20c3e386e1fba1331aa0d25119cc1c3c 100644 (file)
--- 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 (file)
index 0000000..9b3085d
--- /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 <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_out.h>
+#include <grid/tria_iterator.h>
+#include <dofs/hp_dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe_q.h>
+
+#include <fstream>
+
+
+template <int dim>
+void check_cells (const hp::DoFHandler<dim> &dof_handler)
+{
+  for (typename hp::DoFHandler<dim>::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; i<dof_handler.get_fe().size(); ++i)
+        deallog << cell->fe_index_is_active(i);
+      deallog << std::endl;
+
+      Assert (cell->n_active_fe_indices () == 1, ExcInternalError());
+    }
+}
+
+
+void check_faces (const hp::DoFHandler<1> &)
+{}
+
+
+template <int dim>
+void check_faces (const hp::DoFHandler<dim> &dof_handler)
+{
+  for (typename hp::DoFHandler<dim>::active_cell_iterator
+         cell=dof_handler.begin_active();
+       cell!=dof_handler.end(); ++cell)
+    for (unsigned int f=0; f<GeometryInfo<dim>::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; i<dof_handler.get_fe().size(); ++i)
+          deallog << cell->face(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 <int dim>
+void check_edges (const hp::DoFHandler<dim> &dof_handler)
+{
+  for (typename hp::DoFHandler<dim>::active_cell_iterator
+         cell=dof_handler.begin_active();
+       cell!=dof_handler.end(); ++cell)
+    for (unsigned int e=0; e<GeometryInfo<dim>::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; i<dof_handler.get_fe().size(); ++i)
+          deallog << cell->line(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 <int dim>
+void test ()
+{
+  Triangulation<dim> 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<dim> fe_collection;
+  fe_collection.push_back(FE_Q<dim> (1));
+  fe_collection.push_back(FE_Q<dim> (2));
+  fe_collection.push_back(FE_Q<dim> (3));
+  fe_collection.push_back(FE_Q<dim> (4));
+
+  hp::DoFHandler<dim> dof_handler(tria);
+
+  for (typename hp::DoFHandler<dim>::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;
+}