]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Completely rewrite the base class for the hp::FE*Values classes to allow for arbitrar...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 22 Feb 2006 07:00:03 +0000 (07:00 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 22 Feb 2006 07:00:03 +0000 (07:00 +0000)
git-svn-id: https://svn.dealii.org/trunk@12455 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/hp_fe_values.h
deal.II/deal.II/source/fe/hp_fe_values.cc
deal.II/examples/step-21/step-21.cc

index e2d6871c3960b0c76d78adc1a95a98fae4fdb3b4..19d86c94cf44abe6e8f9fe1965aab39d525cc3fc 100644 (file)
@@ -32,130 +32,23 @@ namespace internal
   namespace hp
   {
 /**
- * Map between finite element objects and @p FEValues (the second
- * template parameter, which might be <tt>FEValues<dim></tt>,
- * <tt>FEFaceValues<dim></tt>, ...) objects. The <tt>hpFE*Values</tt> classes
- * use this to hold an <tt>FE*Values</tt> object for each finite element
- * that is used in the triangulation that it integrates on.
- * 
- * @ingroup hp
- */
-    template <int dim, class FEValues>
-    class FEValuesMap
-    {
-      public:
-                                         /**
-                                          * Make destructor virtual,
-                                          * since we have virtual
-                                          * functions in this class.
-                                          */
-        virtual ~FEValuesMap ();
-      
-                                         /**
-                                          * Return a reference to the
-                                          * @p FEValues object selected
-                                          * by the last call to
-                                          * @p select_fe_values. The
-                                          * returned value is a constant
-                                          * reference since the only
-                                          * non-const function in the
-                                          * underlying class would be
-                                          * @p reinit, which you must
-                                          * not call directly any way;
-                                          * rather, use the @p reinit
-                                          * function of the
-                                          * <tt>hpFE*Values</tt> class.
-                                          */
-        const FEValues & get_present_fe_values () const;
-
-      protected:
-                                         /**
-                                          * Create an object of type
-                                          * @p FEValues for this
-                                          * particular finite
-                                          * element. Since the type of
-                                          * @p FEValues is not known
-                                          * here, we have to leave this
-                                          * to derived classes.
-                                          */
-        virtual
-        FEValues *
-        create_fe_values (const FiniteElement<dim> &fe,
-                          const unsigned int active_fe_index) const = 0;
-      
-                                         /**
-                                          * Select the @p FEValues
-                                          * object corresponding to the
-                                          * finite element object given
-                                          * as argument. If there is no
-                                          * such @p FEValues object
-                                          * yet, then one is created by
-                                          * calling @p create_fe_values
-                                          * with this finite element
-                                          * object.
-                                          *
-                                          * A non-const reference to
-                                          * this object is returned.
-                                          */
-        FEValues &
-        select_fe_values (const FiniteElement<dim> &fe,
-                          const unsigned int active_fe_index);
-
-
-                                         /**
-                                          * This field remembers the
-                                          * @p{active_fe_index} of the
-                                          * cell, which was used for
-                                          * the last call to @p{reinit}.
-                                          */
-        unsigned int present_fe_index;
-
-      private:
-                                         /**
-                                          * Have a map from pointers to
-                                          * finite elements to objects
-                                          * of type @p FEValues.
-                                          *
-                                          * Note that the shared pointer
-                                          * will make sure that the
-                                          * memory previously allocated
-                                          * is released when the
-                                          * destructor of this class is
-                                          * run.
-                                          */
-        std::map<std::pair<SmartPointer<const FiniteElement<dim> >, unsigned int>,
-                 boost::shared_ptr<FEValues> > fe_to_fe_values_map;
-
-                                         /**
-                                          * Pointer to the @p FEValues
-                                          * object that is used for the
-                                          * present cell. This always
-                                          * points to one of the objects
-                                          * in the map above, and to
-                                          * which one is determined by
-                                          * the last call to
-                                          * @p select_fe_values.
-                                          */
-        boost::shared_ptr<FEValues> present_fe_values;
-
-    };
-
-
-/**
- * Base class for the <tt>hp::FE*Values</tt> classes, storing the data that
- * is common to them. The first template parameter denotes the space
- * dimension we are in, the second the dimensionality of the object
- * that we integrate on, i.e. for usual @p hp::FEValues it is equal to
- * the first one, while for face integration it is one less.
+ * Base class for the <tt>hp::FE*Values</tt> classes, storing the data
+ * that is common to them. The first template parameter denotes the
+ * space dimension we are in, the second the dimensionality of the
+ * object that we integrate on, i.e. for usual @p hp::FEValues it is
+ * equal to the first one, while for face integration it is one
+ * less. The third template parameter indicates the type of underlying
+ * non-hp FE*Values base type, i.e. it could either be ::FEValues,
+ * ::FEFaceValues, or ::FESubfaceValues.
  * 
  * @ingroup hp
  *
  * @author Wolfgang Bangerth, 2003
  */
-    template <int dim, int q_dim>
+    template <int dim, int q_dim, class FEValues>
     class FEValuesBase 
     {
-      public:
+      public:        
                                          /**
                                           * Constructor. Set the fields
                                           * of this class to the values
@@ -163,10 +56,9 @@ namespace internal
                                           * to the constructor.
                                           */
         FEValuesBase (const ::hp::MappingCollection<dim> &mapping_collection,
+                      const ::hp::FECollection<dim>      &fe_collection,
                       const ::hp::QCollection<q_dim>     &q_collection,
                       const UpdateFlags             update_flags);
-
-
                                          /**
                                           * Constructor. Set the fields
                                           * of this class to the values
@@ -176,18 +68,59 @@ namespace internal
                                           * object for the mapping
                                           * object.
                                           */
-        FEValuesBase (const ::hp::QCollection<q_dim> &q_collection,
+        FEValuesBase (const ::hp::FECollection<dim> &fe_collection,
+                      const ::hp::QCollection<q_dim> &q_collection,
                       const UpdateFlags         update_flags);
+        
+                                         /**
+                                          * Return a reference to the
+                                          * @p FEValues object
+                                          * selected by the last call
+                                          * to
+                                          * select_fe_values(). select_fe_values()
+                                          * in turn is called when you
+                                          * called the @p reinit
+                                          * function of the
+                                          * <tt>hp::FE*Values</tt>
+                                          * class the last time.
+                                          */
+        const FEValues & get_present_fe_values () const;
 
+      protected:
+
+                                         /**
+                                          * Select a FEValues object
+                                          * suitable for the given FE,
+                                          * quadrature, and mapping
+                                          * indices. If such an object
+                                          * doesn't yet exist, create
+                                          * one.
+                                          *
+                                          * The function returns a
+                                          * writable reference so that
+                                          * derived classes can also
+                                          * reinit() the selected
+                                          * FEValues object.
+                                          */
+        FEValues &
+        select_fe_values (const unsigned int fe_index,
+                          const unsigned int mapping_index,
+                          const unsigned int q_index);
 
       protected:
                                          /**
-                                          * A copy of the hp::MappingCollection
-                                          * object, which was specified
-                                          * upon construction of the
-                                          * object.
+                                          * A pointer to the
+                                          * collection of finite
+                                          * elements to be used.
+                                          */
+        const SmartPointer<const ::hp::FECollection<dim> > fe_collection;
+        
+                                         /**
+                                          * A pointer to the
+                                          * collection of mappings to
+                                          * be used.
                                           */
-        SmartPointer<const ::hp::MappingCollection<dim> > mapping_collection;
+        const SmartPointer<const ::hp::MappingCollection<dim> > mapping_collection;
     
                                          /**
                                           * Copy of the quadrature
@@ -197,6 +130,42 @@ namespace internal
                                           */
         const ::hp::QCollection<q_dim> q_collection;
 
+      private:
+                                         /**
+                                          * A table in which we store
+                                          * pointers to fe_values
+                                          * objects for different
+                                          * finite element, mapping,
+                                          * and quadrature objects
+                                          * from our collection. The
+                                          * first index indicates the
+                                          * index of the finite
+                                          * element within the
+                                          * fe_collection, the second
+                                          * the index of the mapping
+                                          * within the mapping
+                                          * collection, and the last
+                                          * one the index of the
+                                          * quadrature formula within
+                                          * the q_collection.
+                                          *
+                                          * Initially, all entries
+                                          * have zero pointers, and we
+                                          * will allocate them lazily
+                                          * as needed in
+                                          * select_fe_values().
+                                          */
+        Table<3,boost::shared_ptr<FEValues> > fe_values_table;
+
+                                         /**
+                                          * Set of indices pointing at
+                                          * the fe_values object
+                                          * selected last time the
+                                          * select_fe_value() function
+                                          * was called.
+                                          */
+        TableIndices<3> present_fe_values_index;
+        
                                          /**
                                           * Values of the update flags as
                                           * given to the constructor.
@@ -217,8 +186,7 @@ namespace hp
  * @ingroup hp
  */  
   template <int dim>
-  class FEValues : public internal::hp::FEValuesMap<dim,::FEValues<dim> >,
-                   protected internal::hp::FEValuesBase<dim,dim>
+  class FEValues : public internal::hp::FEValuesBase<dim,dim,::FEValues<dim> >
   {
     public:
                                        /**
@@ -276,26 +244,108 @@ namespace hp
 
                                        /**
                                         * Reinitialize the object for
-                                        * the given cell. This selects
-                                        * the right @p FEValues object
-                                        * for the finite element in use
-                                        * by the cell given, and calling
-                                        * the @p reinit function on
-                                        * this object.
+                                        * the given cell.
+                                        *
+                                        * After the call, you can get
+                                        * an FEValues object using the
+                                        * get_present_fe_values()
+                                        * function that corresponds to
+                                        * the present cell. For this
+                                        * FEValues object, we use the
+                                        * additional arguments
+                                        * described below to determine
+                                        * which finite element,
+                                        * mapping, and quadrature
+                                        * formula to use. They are
+                                        * order in such a way that the
+                                        * arguments one may want to
+                                        * change most frequently come
+                                        * first. The rules for these
+                                        * arguments are as follows:
+                                        *
+                                        * If the @p fe_index argument
+                                        * to this function is left at
+                                        * its default value, then we
+                                        * use that finite element
+                                        * within the hp::FECollection
+                                        * passed to the constructor of
+                                        * this class with index given
+                                        * by
+                                        * <code>cell-@>active_fe_index()</code>. Consequently,
+                                        * the hp::FECollection
+                                        * argument given to this
+                                        * object should really be the
+                                        * same as that used in the
+                                        * construction of the
+                                        * hp::DofHandler associated
+                                        * with the present cell. On
+                                        * the other hand, if a value
+                                        * is given for this argument,
+                                        * it overrides the choice of
+                                        * <code>cell-@>active_fe_index()</code>.
+                                        *
+                                        * If the @p q_index argument
+                                        * is left at its default
+                                        * value, then we use that
+                                        * quadrature formula within
+                                        * the hp::QCollection passed
+                                        * to the constructor of this
+                                        * class with index given by
+                                        * <code>cell-@>active_fe_index()</code>,
+                                        * i.e. the same index as that
+                                        * of the finite element. In
+                                        * this case, there should be a
+                                        * corresponding quadrature
+                                        * formula for each finite
+                                        * element in the
+                                        * hp::FECollection. As a
+                                        * special case, if the
+                                        * quadrature collection
+                                        * contains only a single
+                                        * element (a frequent case if
+                                        * one wants to use the same
+                                        * quadrature object for all
+                                        * finite elements in an hp
+                                        * discretization, even if that
+                                        * may not be the most
+                                        * efficient), then this single
+                                        * quadrature is used unless a
+                                        * different value for this
+                                        * argument is specified. On
+                                        * the other hand, if a value
+                                        * is given for this argument,
+                                        * it overrides the choice of
+                                        * <code>cell-@>active_fe_index()</code>
+                                        * or the choice for the single
+                                        * quadrature.
+                                        *
+                                        * If the @p mapping_index
+                                        * argument is left at its
+                                        * default value, then we use
+                                        * that mapping object within
+                                        * the hp::MappingCollection
+                                        * passed to the constructor of
+                                        * this class with index given
+                                        * by
+                                        * <code>cell-@>active_fe_index()</code>,
+                                        * i.e. the same index as that
+                                        * of the finite
+                                        * element. As above, if the
+                                        * mapping collection contains
+                                        * only a single element (a
+                                        * frequent case if one wants
+                                        * to use a MappingQ1 object
+                                        * for all finite elements in
+                                        * an hp discretization), then
+                                        * this single mapping is used
+                                        * unless a different value for
+                                        * this argument is specified.
                                         */
       void
-      reinit (const typename hp::DoFHandler<dim>::cell_iterator &cell);
-
-    protected:
-                                       /**
-                                        * Create an object of type
-                                        * @p FEValues for this
-                                        * particular finite element.
-                                        */
-      virtual
-      ::FEValues<dim> *
-      create_fe_values (const FiniteElement<dim> &fe,
-                        const unsigned int active_fe_index) const;
+      reinit (const typename hp::DoFHandler<dim>::cell_iterator &cell,
+              const unsigned int q_index = deal_II_numbers::invalid_unsigned_int,
+              const unsigned int mapping_index = deal_II_numbers::invalid_unsigned_int,
+              const unsigned int fe_index = deal_II_numbers::invalid_unsigned_int);
   };
 
 
@@ -305,8 +355,7 @@ namespace hp
  * @ingroup hp
  */  
   template <int dim>
-  class FEFaceValues : public internal::hp::FEValuesMap<dim,::FEFaceValues<dim> >,
-                       protected internal::hp::FEValuesBase<dim,dim-1>
+  class FEFaceValues : public internal::hp::FEValuesBase<dim,dim-1,::FEFaceValues<dim> >
   {
     public:
                                        /**
@@ -358,57 +407,114 @@ namespace hp
                                         * function.
                                         */
       FEFaceValues (const hp::FECollection<dim>  &fe_collection,
-                    const hp::QCollection<dim-1>     &q_collection,
+                    const hp::QCollection<dim-1> &q_collection,
                     const UpdateFlags             update_flags);
 
-
-                                       /**
-                                        * Reinitialize the object for
-                                        * the given cell. This selects
-                                        * the right @p FEValues object
-                                        * for the finite element in use
-                                        * by the cell given, and calling
-                                        * the @p reinit function on
-                                        * this object.
-                                        */
-      void
-      reinit (const typename hp::DoFHandler<dim>::cell_iterator &cell,
-              const unsigned int face_no);
-
-
                                        /**
                                         * Reinitialize the object for
-                                        * the given cell. In this
-                                        * method, the user can specify
-                                        * which active_fe_index
-                                        * should be used to initialise
-                                        * the underlying FEValues
-                                        * object. This functionality
-                                        * is required, if the face terms
-                                        * between two cells with different
-                                        * polynomial degree should be
-                                        * assembled. In this case the
-                                        * values on the face of the
-                                        * lower order element have to
-                                        * be evaluated at the quadratrure
-                                        * points of the higher order
-                                        * element.
+                                        * the given cell and face.
+                                        *
+                                        * After the call, you can get
+                                        * an FEFaceValues object using the
+                                        * get_present_fe_values()
+                                        * function that corresponds to
+                                        * the present cell. For this
+                                        * FEFaceValues object, we use the
+                                        * additional arguments
+                                        * described below to determine
+                                        * which finite element,
+                                        * mapping, and quadrature
+                                        * formula to use. They are
+                                        * order in such a way that the
+                                        * arguments one may want to
+                                        * change most frequently come
+                                        * first. The rules for these
+                                        * arguments are as follows:
+                                        *
+                                        * If the @p fe_index argument
+                                        * to this function is left at
+                                        * its default value, then we
+                                        * use that finite element
+                                        * within the hp::FECollection
+                                        * passed to the constructor of
+                                        * this class with index given
+                                        * by
+                                        * <code>cell-@>active_fe_index()</code>. Consequently,
+                                        * the hp::FECollection
+                                        * argument given to this
+                                        * object should really be the
+                                        * same as that used in the
+                                        * construction of the
+                                        * hp::DofHandler associated
+                                        * with the present cell. On
+                                        * the other hand, if a value
+                                        * is given for this argument,
+                                        * it overrides the choice of
+                                        * <code>cell-@>active_fe_index()</code>.
+                                        *
+                                        * If the @p q_index argument
+                                        * is left at its default
+                                        * value, then we use that
+                                        * quadrature formula within
+                                        * the hp::QCollection passed
+                                        * to the constructor of this
+                                        * class with index given by
+                                        * <code>cell-@>active_fe_index()</code>,
+                                        * i.e. the same index as that
+                                        * of the finite element. In
+                                        * this case, there should be a
+                                        * corresponding quadrature
+                                        * formula for each finite
+                                        * element in the
+                                        * hp::FECollection. As a
+                                        * special case, if the
+                                        * quadrature collection
+                                        * contains only a single
+                                        * element (a frequent case if
+                                        * one wants to use the same
+                                        * quadrature object for all
+                                        * finite elements in an hp
+                                        * discretization, even if that
+                                        * may not be the most
+                                        * efficient), then this single
+                                        * quadrature is used unless a
+                                        * different value for this
+                                        * argument is specified. On
+                                        * the other hand, if a value
+                                        * is given for this argument,
+                                        * it overrides the choice of
+                                        * <code>cell-@>active_fe_index()</code>
+                                        * or the choice for the single
+                                        * quadrature.
+                                        *
+                                        * If the @p mapping_index
+                                        * argument is left at its
+                                        * default value, then we use
+                                        * that mapping object within
+                                        * the hp::MappingCollection
+                                        * passed to the constructor of
+                                        * this class with index given
+                                        * by
+                                        * <code>cell-@>active_fe_index()</code>,
+                                        * i.e. the same index as that
+                                        * of the finite
+                                        * element. As above, if the
+                                        * mapping collection contains
+                                        * only a single element (a
+                                        * frequent case if one wants
+                                        * to use a MappingQ1 object
+                                        * for all finite elements in
+                                        * an hp discretization), then
+                                        * this single mapping is used
+                                        * unless a different value for
+                                        * this argument is specified.
                                         */
       void
       reinit (const typename hp::DoFHandler<dim>::cell_iterator &cell,
               const unsigned int face_no,
-              const unsigned int active_fe_index);
-
-    protected:
-                                       /**
-                                        * Create an object of type
-                                        * @p FEValues for this
-                                        * particular finite element.
-                                        */
-      virtual
-      ::FEFaceValues<dim> *
-      create_fe_values (const FiniteElement<dim> &fe,
-                        const unsigned int active_fe_index) const;
+              const unsigned int q_index = deal_II_numbers::invalid_unsigned_int,
+              const unsigned int mapping_index = deal_II_numbers::invalid_unsigned_int,
+              const unsigned int fe_index = deal_II_numbers::invalid_unsigned_int);
   };
 
 
@@ -418,8 +524,7 @@ namespace hp
  * @ingroup hp
  */  
   template <int dim>
-  class FESubfaceValues : public internal::hp::FEValuesMap<dim,::FESubfaceValues<dim> >,
-                          protected internal::hp::FEValuesBase<dim,dim-1>
+  class FESubfaceValues : public internal::hp::FEValuesBase<dim,dim-1,::FESubfaceValues<dim> >
   {
     public:
                                        /**
@@ -474,57 +579,91 @@ namespace hp
                        const hp::QCollection<dim-1>    &q_collection,
                        const UpdateFlags            update_flags);
 
-
                                        /**
                                         * Reinitialize the object for
-                                        * the given cell. This selects
-                                        * the right @p FEValues object
-                                        * for the finite element in use
-                                        * by the cell given, and calling
-                                        * the @p reinit function on
-                                        * this object.
-                                        */
-      void
-      reinit (const typename hp::DoFHandler<dim>::cell_iterator &cell,
-              const unsigned int face_no,
-              const unsigned int subface_no);
-
-
-                                       /**
-                                        * Reinitialize the object for
-                                        * the given cell. In this
-                                        * method, the user can specify
-                                        * which active_fe_index
-                                        * should be used to initialise
-                                        * the underlying FEValues
-                                        * object. This functionality
-                                        * is required, if the face terms
-                                        * between two cells with different
-                                        * polynomial degree should be
-                                        * assembled. In this case the
-                                        * values on the face of the
-                                        * lower order element have to
-                                        * be evaluated at the quadratrure
-                                        * points of the higher order
-                                        * element.
+                                        * the given cell, face, and subface.
+                                        *
+                                        * After the call, you can get
+                                        * an FESubfaceValues object using the
+                                        * get_present_fe_values()
+                                        * function that corresponds to
+                                        * the present cell. For this
+                                        * FESubfaceValues object, we use the
+                                        * additional arguments
+                                        * described below to determine
+                                        * which finite element,
+                                        * mapping, and quadrature
+                                        * formula to use. They are
+                                        * order in such a way that the
+                                        * arguments one may want to
+                                        * change most frequently come
+                                        * first. The rules for these
+                                        * arguments are as follows:
+                                        *
+                                        * If the @p q_index argument
+                                        * is left at its default
+                                        * value, then we use that
+                                        * quadrature formula within
+                                        * the hp::QCollection passed
+                                        * to the constructor of this
+                                        * class with index given by
+                                        * <code>cell-@>active_fe_index()</code>,
+                                        * i.e. the same index as that
+                                        * of the finite element. In
+                                        * this case, there should be a
+                                        * corresponding quadrature
+                                        * formula for each finite
+                                        * element in the
+                                        * hp::FECollection. As a
+                                        * special case, if the
+                                        * quadrature collection
+                                        * contains only a single
+                                        * element (a frequent case if
+                                        * one wants to use the same
+                                        * quadrature object for all
+                                        * finite elements in an hp
+                                        * discretization, even if that
+                                        * may not be the most
+                                        * efficient), then this single
+                                        * quadrature is used unless a
+                                        * different value for this
+                                        * argument is specified. On
+                                        * the other hand, if a value
+                                        * is given for this argument,
+                                        * it overrides the choice of
+                                        * <code>cell-@>active_fe_index()</code>
+                                        * or the choice for the single
+                                        * quadrature.
+                                        *
+                                        * If the @p mapping_index
+                                        * argument is left at its
+                                        * default value, then we use
+                                        * that mapping object within
+                                        * the hp::MappingCollection
+                                        * passed to the constructor of
+                                        * this class with index given
+                                        * by
+                                        * <code>cell-@>active_fe_index()</code>,
+                                        * i.e. the same index as that
+                                        * of the finite
+                                        * element. As above, if the
+                                        * mapping collection contains
+                                        * only a single element (a
+                                        * frequent case if one wants
+                                        * to use a MappingQ1 object
+                                        * for all finite elements in
+                                        * an hp discretization), then
+                                        * this single mapping is used
+                                        * unless a different value for
+                                        * this argument is specified.
                                         */
       void
       reinit (const typename hp::DoFHandler<dim>::cell_iterator &cell,
               const unsigned int face_no,
               const unsigned int subface_no,
-              const unsigned int active_fe_index);
-
-
-    protected:
-                                       /**
-                                        * Create an object of type
-                                        * @p FEValues for this
-                                        * particular finite element.
-                                        */
-      virtual
-      ::FESubfaceValues<dim> *
-      create_fe_values (const FiniteElement<dim> &fe,
-                        const unsigned int active_fe_index) const;
+              const unsigned int q_index = deal_II_numbers::invalid_unsigned_int,
+              const unsigned int mapping_index = deal_II_numbers::invalid_unsigned_int,
+              const unsigned int fe_index = deal_II_numbers::invalid_unsigned_int);
   };
   
 }
@@ -536,12 +675,12 @@ namespace internal
 {
   namespace hp
   {
-    template <int dim, class FEValues>
+    template <int dim, int q_dim, class FEValues>
     inline
     const FEValues &
-    FEValuesMap<dim,FEValues>::get_present_fe_values () const
+    FEValuesBase<dim,q_dim,FEValues>::get_present_fe_values () const
     {
-      return *present_fe_values;
+      return *fe_values_table(present_fe_values_index);
     }
   }
   
index ad8fc8b421b2fc6c6aecd4cf3bc339b2cf685cd3..9d38944f0e9ed65830885a406838217de9fd2e05 100644 (file)
 #include <fe/mapping_q1.h>
 
 
-// -------------------------- FEValuesMap -------------------------
-
 namespace internal
 {
   
   namespace hp
   {
-    template <int dim, class FEValues>
-    FEValuesMap<dim,FEValues>::~FEValuesMap () 
-    {}
-  
-  
-    template <int dim, class FEValues>
-    FEValues &
-    FEValuesMap<dim,FEValues>::select_fe_values (const FiniteElement<dim> &fe,
-                                                 const unsigned int active_fe_index)
-    {
-      std::pair<SmartPointer<const FiniteElement<dim> >, unsigned int> fe_pair=
-       std::make_pair(&fe, active_fe_index);
-                                       // check if the finite element
-                                       // does not exist as a key in the
-                                       // map
-      if (fe_to_fe_values_map.find (fe_pair) == fe_to_fe_values_map.end())
-                                         // a-ha! doesn't yet, so let's
-                                         // make it up
-        fe_to_fe_values_map[fe_pair]
-          = boost::shared_ptr<FEValues> (create_fe_values (fe, active_fe_index));
-
-
-                                       // now there definitely is one!
-      present_fe_values = fe_to_fe_values_map[fe_pair];
-
-      return *present_fe_values;
-    }
-
-
 // -------------------------- FEValuesBase -------------------------
 
-    template <int dim, int q_dim>
-    FEValuesBase<dim,q_dim>::FEValuesBase (
-      const ::hp::MappingCollection<dim> &mapping_collection,
-      const ::hp::QCollection<q_dim>     &q_collection,
-      const UpdateFlags             update_flags)
+    template <int dim, int q_dim, class FEValues>
+    FEValuesBase<dim,q_dim,FEValues>::
+    FEValuesBase (const ::hp::MappingCollection<dim> &mapping_collection,
+                  const ::hp::FECollection<dim>      &fe_collection,
+                  const ::hp::QCollection<q_dim>     &q_collection,
+                  const UpdateFlags                   update_flags)
                     :
+                    fe_collection (&fe_collection),
                     mapping_collection (&mapping_collection),
                     q_collection (q_collection),
+                    fe_values_table (fe_collection.size(),
+                                     mapping_collection.size(),
+                                     q_collection.size()),
+                    present_fe_values_index (deal_II_numbers::invalid_unsigned_int,
+                                             deal_II_numbers::invalid_unsigned_int,
+                                             deal_II_numbers::invalid_unsigned_int),
                     update_flags (update_flags)
     {}
 
 
-    template <int dim, int q_dim>
-    FEValuesBase<dim,q_dim>::FEValuesBase (const ::hp::QCollection<q_dim> &q_collection,
-                                           const UpdateFlags         update_flags)
+    template <int dim, int q_dim, class FEValues>
+    FEValuesBase<dim,q_dim,FEValues>::
+    FEValuesBase (const ::hp::FECollection<dim>      &fe_collection,
+                  const ::hp::QCollection<q_dim> &q_collection,
+                  const UpdateFlags         update_flags)
                     :
+                    fe_collection (&fe_collection),
                     mapping_collection (&::hp::StaticMappingQ1<dim>::mapping_collection),
                     q_collection (q_collection),
+                    fe_values_table (fe_collection.size(),
+                                     1,
+                                     q_collection.size()),
+                    present_fe_values_index (deal_II_numbers::invalid_unsigned_int,
+                                             deal_II_numbers::invalid_unsigned_int,
+                                             deal_II_numbers::invalid_unsigned_int),
                     update_flags (update_flags)
     {}
-    
+
+
+
+    template <int dim, int q_dim, class FEValues>
+    FEValues &
+    FEValuesBase<dim,q_dim,FEValues>::
+    select_fe_values (const unsigned int fe_index,
+                      const unsigned int mapping_index,
+                      const unsigned int q_index)
+    {
+      Assert (fe_index < fe_collection->size(),
+              ExcIndexRange (fe_index, 0, fe_collection->size()));
+      Assert (mapping_index < mapping_collection->size(),
+              ExcIndexRange (mapping_index, 0, mapping_collection->size()));
+      Assert (q_index < q_collection.size(),
+              ExcIndexRange (q_index, 0, q_collection.size()));
+
+
+                                       // set the triple of indices
+                                       // that we want to work with
+      present_fe_values_index = TableIndices<3> (fe_index,
+                                                 mapping_index,
+                                                 q_index);
+                                       
+                                       // first check whether we
+                                       // already have an object for
+                                       // this particular combination
+                                       // of indices
+      if (fe_values_table(present_fe_values_index) == 0)
+        fe_values_table(present_fe_values_index)
+          =
+          boost::shared_ptr<FEValues>
+          (new FEValues ((*mapping_collection)[mapping_index],
+                         (*fe_collection)[fe_index],
+                         q_collection[q_index],
+                         update_flags));
+
+                                       // now there definitely is one!
+      return *fe_values_table(present_fe_values_index);
+    }    
   }
 }
 
@@ -88,44 +113,70 @@ namespace hp
 
   template <int dim>
   FEValues<dim>::FEValues (const hp::MappingCollection<dim> &mapping,
-                           const hp::FECollection<dim>  &/*fe_collection*/,
+                           const hp::FECollection<dim>      &fe_collection,
                            const hp::QCollection<dim>       &q_collection,
-                           const UpdateFlags             update_flags)
+                           const UpdateFlags                 update_flags)
                   :
-                  internal::hp::FEValuesBase<dim,dim> (mapping,
-                                                       q_collection,
-                                                       update_flags)
+                  internal::hp::FEValuesBase<dim,dim,::FEValues<dim> > (mapping,
+                                                                  fe_collection,
+                                                                  q_collection,
+                                                                  update_flags)
   {}
 
 
   template <int dim>
-  FEValues<dim>::FEValues (const hp::FECollection<dim> &/*fe_collection*/,
+  FEValues<dim>::FEValues (const hp::FECollection<dim> &fe_collection,
                            const hp::QCollection<dim>      &q_collection,
                            const UpdateFlags            update_flags)
                   :
-                  internal::hp::FEValuesBase<dim,dim> (q_collection,
+                  internal::hp::FEValuesBase<dim,dim,::FEValues<dim> > (fe_collection,
+                                                       q_collection,
                                                        update_flags)
   {}
 
 
   template <int dim>
   void
-  FEValues<dim>::reinit (const typename hp::DoFHandler<dim>::cell_iterator &cell)
-  {
-    this->present_fe_index = cell->active_fe_index ();
-    this->select_fe_values (cell->get_fe(), this->present_fe_index).reinit (cell);
-  }
-
-
-
-  template <int dim>
-  ::FEValues<dim> *
-  FEValues<dim>::create_fe_values (const FiniteElement<dim> &fe,
-                                   const unsigned int active_fe_index) const
+  FEValues<dim>::reinit (const typename hp::DoFHandler<dim>::cell_iterator &cell,
+                         const unsigned int q_index,
+                         const unsigned int mapping_index,
+                         const unsigned int fe_index)
   {
-    return new ::FEValues<dim> (
-      (*this->mapping_collection)[active_fe_index], fe,
-      this->q_collection[active_fe_index], this->update_flags);
+                                     // determine which indices we
+                                     // should actually use
+    unsigned int real_q_index       = q_index,
+                 real_mapping_index = mapping_index,
+                 real_fe_index      = fe_index;
+
+    if (real_q_index == deal_II_numbers::invalid_unsigned_int)
+      if (q_collection.size() > 1)
+        real_q_index = cell->active_fe_index();
+      else
+        real_q_index = 0;
+    
+    if (real_mapping_index == deal_II_numbers::invalid_unsigned_int)
+      if (mapping_collection->size() > 1)
+        real_mapping_index = cell->active_fe_index();
+      else
+        real_mapping_index = 0;
+
+    if (real_fe_index == deal_II_numbers::invalid_unsigned_int)
+      real_fe_index = cell->active_fe_index();
+
+                                     // some checks
+    Assert (real_q_index < q_collection.size(),
+            ExcIndexRange (real_q_index, 0, q_collection.size()));
+    Assert (real_mapping_index < mapping_collection->size(),
+            ExcIndexRange (real_mapping_index, 0, mapping_collection->size()));
+    Assert (real_fe_index < fe_collection->size(),
+            ExcIndexRange (real_fe_index, 0, fe_collection->size()));
+    
+                                     // now finally actually get the
+                                     // corresponding object and
+                                     // initialize it
+    this->select_fe_values (real_fe_index,
+                            real_mapping_index,
+                            real_q_index).reinit (cell);
   }
 
 
@@ -134,55 +185,71 @@ namespace hp
 
   template <int dim>
   FEFaceValues<dim>::FEFaceValues (const hp::MappingCollection<dim> &mapping,
-                                   const hp::FECollection<dim>  &/*fe_collection*/,
+                                   const hp::FECollection<dim>  &fe_collection,
                                    const hp::QCollection<dim-1> &q_collection,
                                    const UpdateFlags         update_flags)
                   :
-                  internal::hp::FEValuesBase<dim,dim-1> (mapping,
+                  internal::hp::FEValuesBase<dim,dim-1,::FEFaceValues<dim> > (mapping,
+                                                         fe_collection,
                                                          q_collection,
                                                          update_flags)
   {}
 
 
   template <int dim>
-  FEFaceValues<dim>::FEFaceValues (const hp::FECollection<dim>  &/*fe_collection*/,
+  FEFaceValues<dim>::FEFaceValues (const hp::FECollection<dim>  &fe_collection,
                                    const hp::QCollection<dim-1> &q_collection,
                                    const UpdateFlags         update_flags)
                   :
-                  internal::hp::FEValuesBase<dim,dim-1> (q_collection,
+                  internal::hp::FEValuesBase<dim,dim-1,::FEFaceValues<dim> > (fe_collection,
+                                                         q_collection,
                                                          update_flags)
   {}
 
 
-  template <int dim>
-  void
-  FEFaceValues<dim>::reinit (const typename hp::DoFHandler<dim>::cell_iterator &cell,
-                             const unsigned int face_no)
-  {
-    this->present_fe_index = cell->active_fe_index ();
-    this->select_fe_values (cell->get_fe(), this->present_fe_index).reinit (cell, face_no);
-  }
-
-
   template <int dim>
   void
   FEFaceValues<dim>::reinit (const typename hp::DoFHandler<dim>::cell_iterator &cell,
                              const unsigned int face_no,
-                             const unsigned int active_fe_index)
-  {
-    this->present_fe_index = active_fe_index;
-    this->select_fe_values (cell->get_fe(), active_fe_index).reinit (cell, face_no);
-  }
-
-
-  template <int dim>
-  ::FEFaceValues<dim> *
-  FEFaceValues<dim>::create_fe_values (const FiniteElement<dim> &fe,
-                                       const unsigned int active_fe_index) const
+                             const unsigned int q_index,
+                             const unsigned int mapping_index,
+                             const unsigned int fe_index)
   {
-    return new ::FEFaceValues<dim> (
-      (*this->mapping_collection)[active_fe_index], fe,
-      this->q_collection[active_fe_index], this->update_flags);
+                                     // determine which indices we
+                                     // should actually use
+    unsigned int real_q_index       = q_index,
+                 real_mapping_index = mapping_index,
+                 real_fe_index      = fe_index;
+
+    if (real_q_index == deal_II_numbers::invalid_unsigned_int)
+      if (q_collection.size() > 1)
+        real_q_index = cell->active_fe_index();
+      else
+        real_q_index = 0;
+    
+    if (real_mapping_index == deal_II_numbers::invalid_unsigned_int)
+      if (mapping_collection->size() > 1)
+        real_mapping_index = cell->active_fe_index();
+      else
+        real_mapping_index = 0;
+
+    if (real_fe_index == deal_II_numbers::invalid_unsigned_int)
+      real_fe_index = cell->active_fe_index();
+
+                                     // some checks
+    Assert (real_q_index < q_collection.size(),
+            ExcIndexRange (real_q_index, 0, q_collection.size()));
+    Assert (real_mapping_index < mapping_collection->size(),
+            ExcIndexRange (real_mapping_index, 0, mapping_collection->size()));
+    Assert (real_fe_index < fe_collection->size(),
+            ExcIndexRange (real_fe_index, 0, fe_collection->size()));
+    
+                                     // now finally actually get the
+                                     // corresponding object and
+                                     // initialize it
+    this->select_fe_values (real_fe_index,
+                            real_mapping_index,
+                            real_q_index).reinit (cell, face_no);
   }
 
 
@@ -191,57 +258,72 @@ namespace hp
 
   template <int dim>
   FESubfaceValues<dim>::FESubfaceValues (const hp::MappingCollection<dim> &mapping,
-                                         const hp::FECollection<dim>  &/*fe_collection*/,
+                                         const hp::FECollection<dim>  &fe_collection,
                                          const hp::QCollection<dim-1> &q_collection,
                                          const UpdateFlags         update_flags)
                   :
-                  internal::hp::FEValuesBase<dim,dim-1> (mapping,
+                  internal::hp::FEValuesBase<dim,dim-1,::FESubfaceValues<dim> > (mapping,
+                                                         fe_collection,
                                                          q_collection,
                                                          update_flags)
   {}
 
 
   template <int dim>
-  FESubfaceValues<dim>::FESubfaceValues (const hp::FECollection<dim>  &/*fe_collection*/,
+  FESubfaceValues<dim>::FESubfaceValues (const hp::FECollection<dim>  &fe_collection,
                                          const hp::QCollection<dim-1> &q_collection,
                                          const UpdateFlags         update_flags)
                   :
-                  internal::hp::FEValuesBase<dim,dim-1> (q_collection,
+                  internal::hp::FEValuesBase<dim,dim-1,::FESubfaceValues<dim> > (fe_collection,
+                                                         q_collection,
                                                          update_flags)
   {}
 
 
-  template <int dim>
-  void
-  FESubfaceValues<dim>::reinit (const typename hp::DoFHandler<dim>::cell_iterator &cell,
-                                const unsigned int face_no,
-                                const unsigned int subface_no)
-  {
-    this->present_fe_index = cell->active_fe_index ();
-    this->select_fe_values (cell->get_fe(), this->present_fe_index).reinit (cell, face_no, subface_no);
-  }
-
-
   template <int dim>
   void
   FESubfaceValues<dim>::reinit (const typename hp::DoFHandler<dim>::cell_iterator &cell,
                                 const unsigned int face_no,
                                 const unsigned int subface_no,
-                                const unsigned int active_fe_index)
+                                const unsigned int q_index,
+                                const unsigned int mapping_index,
+                                const unsigned int fe_index)
   {
-    this->present_fe_index = active_fe_index;
-    this->select_fe_values (cell->get_fe(), active_fe_index).reinit (cell, face_no, subface_no);
-  }
-
-
-  template <int dim>
-  ::FESubfaceValues<dim> *
-  FESubfaceValues<dim>::create_fe_values (const FiniteElement<dim> &fe,
-                                          const unsigned int active_fe_index) const
-  {
-    return new ::FESubfaceValues<dim> (
-      (*this->mapping_collection)[active_fe_index], fe,
-      this->q_collection[active_fe_index], this->update_flags);
+                                     // determine which indices we
+                                     // should actually use
+    unsigned int real_q_index       = q_index,
+                 real_mapping_index = mapping_index,
+                 real_fe_index      = fe_index;
+
+    if (real_q_index == deal_II_numbers::invalid_unsigned_int)
+      if (q_collection.size() > 1)
+        real_q_index = cell->active_fe_index();
+      else
+        real_q_index = 0;
+    
+    if (real_mapping_index == deal_II_numbers::invalid_unsigned_int)
+      if (mapping_collection->size() > 1)
+        real_mapping_index = cell->active_fe_index();
+      else
+        real_mapping_index = 0;
+
+    if (real_fe_index == deal_II_numbers::invalid_unsigned_int)
+      real_fe_index = cell->active_fe_index();
+
+                                     // some checks
+    Assert (real_q_index < q_collection.size(),
+            ExcIndexRange (real_q_index, 0, q_collection.size()));
+    Assert (real_mapping_index < mapping_collection->size(),
+            ExcIndexRange (real_mapping_index, 0, mapping_collection->size()));
+    Assert (real_fe_index < fe_collection->size(),
+            ExcIndexRange (real_fe_index, 0, fe_collection->size()));
+    
+                                     // now finally actually get the
+                                     // corresponding object and
+                                     // initialize it
+    this->select_fe_values (real_fe_index,
+                            real_mapping_index,
+                            real_q_index).reinit (cell, face_no, subface_no);
   }
 }
 
@@ -251,9 +333,13 @@ namespace internal
 {
   namespace hp
   {
-    template class FEValuesBase<deal_II_dimension,deal_II_dimension>;
+    template class FEValuesBase<deal_II_dimension,deal_II_dimension,
+                                ::FEValues<deal_II_dimension> >;
 #if deal_II_dimension >= 2
-    template class FEValuesBase<deal_II_dimension,deal_II_dimension-1>;
+    template class FEValuesBase<deal_II_dimension,deal_II_dimension-1,
+                                ::FEFaceValues<deal_II_dimension> >;
+    template class FEValuesBase<deal_II_dimension,deal_II_dimension-1,
+                                ::FESubfaceValues<deal_II_dimension> >;
 #endif
   }
 }
@@ -266,10 +352,3 @@ namespace hp
   template class FESubfaceValues<deal_II_dimension>;
 #endif
 }
-
-// Putting the following explicit instantiations into the brackets 
-// of the appropriate namespace somehow causes problems with the 
-// Apple gcc3.3. Therefore these are separated.
-template class internal::hp::FEValuesMap<deal_II_dimension,FEValues<deal_II_dimension> >;
-template class internal::hp::FEValuesMap<deal_II_dimension,FEFaceValues<deal_II_dimension> >;
-template class internal::hp::FEValuesMap<deal_II_dimension,FESubfaceValues<deal_II_dimension> >;
index 564f3d787bc3bc3a25c1d82c1cd18d4b243959cb..0e4a21fd0dfa4e366749efe87d5dcf3c79a136da 100644 (file)
@@ -1010,9 +1010,9 @@ void DGMethod<dim>::assemble_system1 ()
                                                       // the face quadrature rule of the
                                                       // higher order element
                                                       // will be used.
-                     const unsigned int use_fe_index = 
-                       neighbor_child->active_fe_index () > cell->active_fe_index () ?
-                       neighbor_child->active_fe_index () : cell->active_fe_index ();
+                     const unsigned int quadrature_index = 
+                       std::max (neighbor_child->active_fe_index (),
+                                  cell->active_fe_index ());
                        
 
                                                       // As these are
@@ -1067,8 +1067,8 @@ void DGMethod<dim>::assemble_system1 ()
                                                       // of the
                                                       // neighboring
                                                       // child cell.
-                     fe_v_subface_x.reinit (cell, face_no, subface_no, use_fe_index);
-                     fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2, use_fe_index);
+                     fe_v_subface_x.reinit (cell, face_no, subface_no, quadrature_index);
+                     fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2, quadrature_index);
 
                      dg.assemble_face_term1(fe_v_subface_x.get_present_fe_values (),
                                             fe_v_face_neighbor_x.get_present_fe_values (),
@@ -1124,9 +1124,9 @@ void DGMethod<dim>::assemble_system1 ()
                                                       // quadrature rule
                                                        // of higher order
                                                        // cell.
-                     const unsigned int use_fe_index = 
-                       neighbor->active_fe_index () > cell->active_fe_index () ?
-                       neighbor->active_fe_index () : cell->active_fe_index ();
+                     const unsigned int quadrature_index = 
+                       std::max (neighbor->active_fe_index (),
+                                  cell->active_fe_index ());
 
                                                       // We reinit
                                                       // the
@@ -1140,8 +1140,8 @@ void DGMethod<dim>::assemble_system1 ()
                                                       // the
                                                       // corresponding
                                                       // face terms.
-                     fe_v_face_x.reinit (cell, face_no, use_fe_index);
-                     fe_v_face_neighbor_x.reinit (neighbor, neighbor2, use_fe_index);
+                     fe_v_face_x.reinit (cell, face_no, quadrature_index);
+                     fe_v_face_neighbor_x.reinit (neighbor, neighbor2, quadrature_index);
                      
                      dg.assemble_face_term1(fe_v_face_x.get_present_fe_values (),
                                             fe_v_face_neighbor_x.get_present_fe_values (),
@@ -1193,9 +1193,9 @@ void DGMethod<dim>::assemble_system1 ()
                                                       // quadrature rule
                                                        // of higher order
                                                        // cell.
-                     const unsigned int use_fe_index = 
-                       neighbor->active_fe_index () > cell->active_fe_index () ?
-                       neighbor->active_fe_index () : cell->active_fe_index ();
+                     const unsigned int quadrature_index = 
+                       std::max (neighbor->active_fe_index (),
+                                  cell->active_fe_index ());
 
                                                       // Reinit the
                                                       // appropriate
@@ -1203,9 +1203,9 @@ void DGMethod<dim>::assemble_system1 ()
                                                       // and assemble
                                                       // the face
                                                       // terms.
-                     fe_v_face_x.reinit (cell, face_no, use_fe_index);
+                     fe_v_face_x.reinit (cell, face_no, quadrature_index);
                      fe_v_subface_neighbor_x.reinit (neighbor, neighbor_face_no,
-                                                     neighbor_subface_no, use_fe_index);
+                                                     neighbor_subface_no, quadrature_index);
                      
                      dg.assemble_face_term1(fe_v_face_x.get_present_fe_values (),
                                             fe_v_subface_neighbor_x.get_present_fe_values (),
@@ -1390,9 +1390,9 @@ void DGMethod<dim>::assemble_system2 ()
                          neighbor2,subface_no));
                      const unsigned int dofs_on_neighbor_child = neighbor_child->get_fe().dofs_per_cell;
 
-                     const unsigned int use_fe_index = 
-                       neighbor_child->active_fe_index () > cell->active_fe_index () ?
-                       neighbor_child->active_fe_index () : cell->active_fe_index ();
+                     const unsigned int quadrature_index = 
+                       std::max (neighbor_child->active_fe_index (),
+                                  cell->active_fe_index ());
 
                      Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
                              ExcInternalError());
@@ -1402,8 +1402,8 @@ void DGMethod<dim>::assemble_system2 ()
                      u_vn_matrix = 0;
                      un_vn_matrix = 0;
                      
-                     fe_v_subface_x.reinit (cell, face_no, subface_no, use_fe_index);
-                     fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2, use_fe_index);
+                     fe_v_subface_x.reinit (cell, face_no, subface_no, quadrature_index);
+                     fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2, quadrature_index);
 
                      dg.assemble_face_term2(fe_v_subface_x.get_present_fe_values (),
                                             fe_v_face_neighbor_x.get_present_fe_values (),
@@ -1441,16 +1441,16 @@ void DGMethod<dim>::assemble_system2 ()
                    {
                      const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
                      
-                     const unsigned int use_fe_index = 
-                       neighbor->active_fe_index () > cell->active_fe_index () ?
-                       neighbor->active_fe_index () : cell->active_fe_index ();
+                     const unsigned int quadrature_index = 
+                       std::max (neighbor->active_fe_index (),
+                                  cell->active_fe_index ());
 
                      un_v_matrix = 0;
                      u_vn_matrix = 0;
                      un_vn_matrix = 0;
                      
-                     fe_v_face_x.reinit (cell, face_no, use_fe_index);
-                     fe_v_face_neighbor_x.reinit (neighbor, neighbor2, use_fe_index);
+                     fe_v_face_x.reinit (cell, face_no, quadrature_index);
+                     fe_v_face_neighbor_x.reinit (neighbor, neighbor2, quadrature_index);
 
                      dg.assemble_face_term2(fe_v_face_x.get_present_fe_values (),
                                             fe_v_face_neighbor_x.get_present_fe_values (),
@@ -1658,7 +1658,7 @@ void DGMethod<dim>::output_results (const unsigned int cycle) const
   Assert (cycle < 10, ExcInternalError());
   
   //  filename += ".gnuplot";
-  filename += ".vtk";
+  filename += ".gmv";
   deallog << "Writing solution to <" << filename << ">..."
            << std::endl << std::endl;
   std::ofstream gnuplot_output (filename.c_str());
@@ -1670,7 +1670,7 @@ void DGMethod<dim>::output_results (const unsigned int cycle) const
   data_out.build_patches (5);
   
 //  data_out.write_gnuplot(gnuplot_output);
-  data_out.write_vtk(gnuplot_output);
+  data_out.write_gmv(gnuplot_output);
 }
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.