]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Make use of the new FE_Poly class. Move some of the renumbering business to the Tenso...
authorhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 16 Mar 2004 16:39:39 +0000 (16:39 +0000)
committerhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 16 Mar 2004 16:39:39 +0000 (16:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@8778 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_q_hierarchical.h
deal.II/deal.II/source/fe/fe_q_hierarchical.cc

index 540fb2a363c00e09d49d6695e26cadd893572029..2ddad0c47cc7b8673d8ab935ee3815baa13ac505 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2002, 2003 by the deal.II authors
+//    Copyright (C) 2002, 2003, 2004 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
 #define __deal2__fe_q_hierarchical_h
 
 #include <base/config.h>
-#include <base/polynomial.h>
 #include <base/tensor_product_polynomials.h>
+#include <fe/fe_poly.h>
 #include <fe/fe.h>
-#include <lac/full_matrix.h>
+//#include <lac/full_matrix.h>
 
 template <int dim> class TensorProductPolynomials;
 template <int dim> class MappingQ;
@@ -231,15 +231,15 @@ template <int dim> class MappingQ;
  * Note the reverse ordering of degrees of freedom on the left and upper
  * line.
  *
- * @author Brian Carnes, 2002
+ * @author Brian Carnes, 2002, Ralf Hartmann 2004
  */
 template <int dim>
-class FE_Q_Hierarchical : public FiniteElement<dim>
+class FE_Q_Hierarchical : public FE_Poly<TensorProductPolynomials<dim>,dim>
 {
   public:
                                     /**
                                      * Constructor for tensor product
-                                     * polynomials of degree @p{p}.
+                                     * polynomials of degree @p p.
                                      */
     FE_Q_Hierarchical (const unsigned int p);
     
@@ -254,106 +254,6 @@ class FE_Q_Hierarchical : public FiniteElement<dim>
                                      */
     virtual std::string get_name () const;
 
-                                    /**
-                                     * Return the value of the
-                                     * @p{i}th shape function at the
-                                     * point @p{p}.  @p{p} is a point
-                                     * on the reference element.
-                                     */
-    virtual double shape_value (const unsigned int i,
-                               const Point<dim> &p) const;
-    
-                                    /**
-                                     * Return the value of the
-                                     * @p{component}th vector
-                                     * component of the @p{i}th shape
-                                     * function at the point
-                                     * @p{p}. See the
-                                     * @ref{FiniteElementBase} base
-                                     * class for more information
-                                     * about the semantics of this
-                                     * function.
-                                     *
-                                     * Since this element is scalar,
-                                     * the returned value is the same
-                                     * as if the function without the
-                                     * @p{_component} suffix were
-                                     * called, provided that the
-                                     * specified component is zero.
-                                     */
-    virtual double shape_value_component (const unsigned int i,
-                                         const Point<dim> &p,
-                                         const unsigned int component) const;
-
-                                    /**
-                                     * Return the gradient of the
-                                     * @p{i}th shape function at the
-                                     * point @p{p}. @p{p} is a point
-                                     * on the reference element, and
-                                     * likewise the gradient is the
-                                     * gradient on the unit cell with
-                                     * respect to unit cell
-                                     * coordinates.
-                                     */
-    virtual Tensor<1,dim> shape_grad (const unsigned int  i,
-                                     const Point<dim>   &p) const;
-
-                                    /**
-                                     * Return the gradient of the
-                                     * @p{component}th vector
-                                     * component of the @p{i}th shape
-                                     * function at the point
-                                     * @p{p}. See the
-                                     * @ref{FiniteElementBase} base
-                                     * class for more information
-                                     * about the semantics of this
-                                     * function.
-                                     *
-                                     * Since this element is scalar,
-                                     * the returned value is the same
-                                     * as if the function without the
-                                     * @p{_component} suffix were
-                                     * called, provided that the
-                                     * specified component is zero.
-                                     */
-    virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
-                                               const Point<dim> &p,
-                                               const unsigned int component) const;
-
-                                    /**
-                                     * Return the tensor of second
-                                     * derivatives of the @p{i}th
-                                     * shape function at point @p{p}
-                                     * on the unit cell. The
-                                     * derivatives are derivatives on
-                                     * the unit cell with respect to
-                                     * unit cell coordinates.
-                                     */
-    virtual Tensor<2,dim> shape_grad_grad (const unsigned int  i,
-                                          const Point<dim> &p) const;
-
-                                    /**
-                                     * Return the second derivative
-                                     * of the @p{component}th vector
-                                     * component of the @p{i}th shape
-                                     * function at the point
-                                     * @p{p}. See the
-                                     * @ref{FiniteElementBase} base
-                                     * class for more information
-                                     * about the semantics of this
-                                     * function.
-                                     *
-                                     * Since this element is scalar,
-                                     * the returned value is the same
-                                     * as if the function without the
-                                     * @p{_component} suffix were
-                                     * called, provided that the
-                                     * specified component is zero.
-                                     */
-    virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
-                                                    const Point<dim> &p,
-                                                    const unsigned int component) const;
-
                                     /**
                                      * Return the polynomial degree
                                      * of this finite element,
@@ -362,34 +262,6 @@ class FE_Q_Hierarchical : public FiniteElement<dim>
                                      */
     unsigned int get_degree () const;
     
-                                     /**
-                                     * Number of base elements in a
-                                     * mixed discretization. Since
-                                     * this is a scalar element,
-                                     * return one.
-                                     */
-    virtual unsigned int n_base_elements () const;
-    
-                                    /**
-                                     * Access to base element
-                                     * objects. Since this element is
-                                     * scalar, @p{base_element(0)} is
-                                     * @p{this}, and all other
-                                     * indices throw an error.
-                                     */
-    virtual const FiniteElement<dim> &
-    base_element (const unsigned int index) const;
-    
-                                     /**
-                                      * Multiplicity of base element
-                                      * @p{index}. Since this is a
-                                      * scalar element,
-                                      * @p{element_multiplicity(0)}
-                                      * returns one, and all other
-                                      * indices will throw an error.
-                                      */
-    virtual unsigned int element_multiplicity (const unsigned int index) const;
-    
                                     /**
                                      * Check for non-zero values on a face.
                                      *
@@ -438,59 +310,6 @@ class FE_Q_Hierarchical : public FiniteElement<dim>
                                      * constructors of @p{FESystem}.
                                      */
     virtual FiniteElement<dim> * clone() const;
-  
-                                    /**
-                                     * Prepare internal data
-                                     * structures and fill in values
-                                     * independent of the cell.
-                                     */
-    virtual
-    typename Mapping<dim>::InternalDataBase *
-    get_data (const UpdateFlags,
-             const Mapping<dim>& mapping,
-             const Quadrature<dim>& quadrature) const ;
-
-                                    /**
-                                     * Implementation of the same
-                                     * function in
-                                     * @ref{FiniteElement}.
-                                     */
-    virtual void
-    fill_fe_values (const Mapping<dim> &mapping,
-                   const typename DoFHandler<dim>::cell_iterator &cell,
-                   const Quadrature<dim>                         &quadrature,
-                   typename Mapping<dim>::InternalDataBase       &mapping_internal,
-                   typename Mapping<dim>::InternalDataBase       &fe_internal,
-                   FEValuesData<dim>& data) const;
-    
-                                    /**
-                                     * Implementation of the same
-                                     * function in
-                                     * @ref{FiniteElement}.
-                                     */
-    virtual void
-    fill_fe_face_values (const Mapping<dim> &mapping,
-                        const typename DoFHandler<dim>::cell_iterator &cell,
-                        const unsigned int                            face_no,
-                        const Quadrature<dim-1>                       &quadrature,
-                        typename Mapping<dim>::InternalDataBase       &mapping_internal,
-                        typename Mapping<dim>::InternalDataBase       &fe_internal,
-                        FEValuesData<dim>& data) const ;
-    
-                                    /**
-                                     * Implementation of the same
-                                     * function in
-                                     * @ref{FiniteElement}.
-                                     */
-    virtual void
-    fill_fe_subface_values (const Mapping<dim> &mapping,
-                           const typename DoFHandler<dim>::cell_iterator &cell,
-                           const unsigned int                            face_no,
-                           const unsigned int                            sub_no,
-                           const Quadrature<dim-1>                       &quadrature,
-                           typename Mapping<dim>::InternalDataBase       &mapping_internal,
-                           typename Mapping<dim>::InternalDataBase       &fe_internal,
-                           FEValuesData<dim>& data) const ;
 
   private:
 
@@ -628,159 +447,16 @@ class FE_Q_Hierarchical : public FiniteElement<dim>
                                      */
     void initialize_unit_face_support_points ();
     
-                                    /**
-                                     * Determine the values that need
-                                     * to be computed on the unit
-                                     * cell to be able to compute all
-                                     * values required by @p{flags}.
-                                     *
-                                     * For the purpuse of this
-                                     * function, refer to the
-                                     * documentation in
-                                     * @p{FiniteElement}.
-                                     *
-                                     * The effect in this element is
-                                     * as follows: if
-                                     * @p{update_values} is set in
-                                     * @p{flags}, copy it to the
-                                     * result. All other flags of the
-                                     * result are cleared, since
-                                     * everything else must be
-                                     * computed for each cell.
-                                     */
-    virtual UpdateFlags update_once (const UpdateFlags flags) const;
-  
-                                    /**
-                                     * Determine the values that need
-                                     * to be computed on every
-                                     * cell to be able to compute all
-                                     * values required by @p{flags}.
-                                     *
-                                     * For the purpuse of this
-                                     * function, refer to the
-                                     * documentation in
-                                     * @p{FiniteElement}.
-                                     *
-                                     * The effect in this element is
-                                     * as follows:
-                                     * @begin{itemize}
-                                     * @item if @p{update_gradients}
-                                     * is set, the result will
-                                     * contain @p{update_gradients}
-                                     * and
-                                     * @p{update_covariant_transformation}.
-                                     * The latter is required to
-                                     * transform the gradient on the
-                                     * unit cell to the real
-                                     * cell. Remark, that the action
-                                     * required by
-                                     * @p{update_covariant_transformation}
-                                     * is actually performed by the
-                                     * @p{Mapping} object used in
-                                     * conjunction with this finite
-                                     * element.
-                                     * @item if
-                                     * @p{update_second_derivatives}
-                                     * is set, the result will
-                                     * contain
-                                     * @p{update_second_derivatives}
-                                     * and
-                                     * @p{update_covariant_transformation}.
-                                     * The rationale is the same as
-                                     * above and no higher
-                                     * derivatives of the
-                                     * transformation are required,
-                                     * since we use difference
-                                     * quotients for the actual
-                                     * computation.
-                                     * @end{itemize}
-                                     */
-    virtual UpdateFlags update_each (const UpdateFlags flags) const;
-    
                                     /**
                                      * Degree of the polynomials.
                                      */  
     const unsigned int degree;
-
-                                    /**
-                                     * Mapping from lexicographic to
-                                     * shape function numbering.
-                                     */
-    const std::vector<unsigned int> renumber;
-
-                                    /**
-                                     * Inverse renumber
-                                     * vector. i.e. mapping from
-                                     * shape function numbering to
-                                     * lexicographic numbering.
-                                     */
-    const std::vector<unsigned int> renumber_inverse;
              
                                     /**
                                      * Mapping from lexicographic to
                                      * shape function numbering on first face.
                                      */
     const std::vector<unsigned int> face_renumber;
-
-                                    /**
-                                     * Pointer to the tensor
-                                     * product polynomials.
-                                     */
-    const TensorProductPolynomials<dim> polynomial_space;
-
-                                    /**
-                                     * Fields of cell-independent data.
-                                     *
-                                     * For information about the
-                                     * general purpose of this class,
-                                     * see the documentation of the
-                                     * base class.
-                                     */
-    class InternalData : public FiniteElementBase<dim>::InternalDataBase
-    {
-      public:
-                                        /**
-                                         * Array with shape function
-                                         * values in quadrature
-                                         * points. There is one
-                                         * row for each shape
-                                         * function, containing
-                                         * values for each quadrature
-                                         * point.
-                                         *
-                                         * In this array, we store
-                                         * the values of the shape
-                                         * function in the quadrature
-                                         * points on the unit
-                                         * cell. Since these values
-                                         * do not change under
-                                         * transformation to the real
-                                         * cell, we only need to copy
-                                         * them over when visiting a
-                                         * concrete cell.
-                                         */
-       Table<2,double> shape_values;
-
-                                        /**
-                                         * Array with shape function
-                                         * gradients in quadrature
-                                         * points. There is one
-                                         * row for each shape
-                                         * function, containing
-                                         * values for each quadrature
-                                         * point.
-                                         *
-                                         * We store the gradients in
-                                         * the quadrature points on
-                                         * the unit cell. We then
-                                         * only have to apply the
-                                         * transformation (which is a
-                                         * matrix-vector
-                                         * multiplication) when
-                                         * visiting an actual cell.
-                                         */      
-       Table<2,Tensor<1,dim> > shape_gradients;
-    };
     
                                     /**
                                      * Allow access from other
index edf384f4091e782b0fe5e31c2637165d7f0c2c8f..7b9908a9d5c525a0d2e2006d9b1a25b43888a268 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2002, 2003 by the deal.II authors
+//    Copyright (C) 2002, 2003, 2004 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -47,17 +47,19 @@ namespace
 template <int dim>
 FE_Q_Hierarchical<dim>::FE_Q_Hierarchical (const unsigned int degree)
                :
-               FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree),1, degree),
-                                   std::vector<bool> (FiniteElementData<dim>(get_dpo_vector(degree),1, degree).dofs_per_cell,
-                                                       false),
-                                   std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),1, degree).dofs_per_cell,
-                                                                   std::vector<bool>(1,true))),
-               degree(degree),
-               renumber(lexicographic_to_hierarchic_numbering (*this, degree)),
-               renumber_inverse(invert_numbering(renumber)),
-               face_renumber(face_lexicographic_to_hierarchic_numbering (degree)),
-               polynomial_space(Polynomials::Hierarchical::generate_complete_basis(degree))
+               FE_Poly<TensorProductPolynomials<dim>, dim> (
+                 Polynomials::Hierarchical::generate_complete_basis(degree),
+                 FiniteElementData<dim>(get_dpo_vector(degree),1, degree),
+                 std::vector<bool> (FiniteElementData<dim>(
+                   get_dpo_vector(degree),1, degree).dofs_per_cell, false),
+                 std::vector<std::vector<bool> >(FiniteElementData<dim>(
+                   get_dpo_vector(degree),1, degree).dofs_per_cell, std::vector<bool>(1,true))),
+                                                   degree(degree),
+                                                   face_renumber(face_lexicographic_to_hierarchic_numbering (degree))
 {
+  this->poly_space.set_numbering(invert_numbering(
+    lexicographic_to_hierarchic_numbering (*this, degree)));
+  
                                   // The matrix @p{dofs_cell} contains the 
                                   // values of the linear functionals of 
                                   // the master 1d cell applied to the 
@@ -133,78 +135,6 @@ FE_Q_Hierarchical<dim>::clone() const
 }
 
 
-
-template <int dim>
-double
-FE_Q_Hierarchical<dim>::shape_value (const unsigned int i,
-                                    const Point<dim> &p) const
-{
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  return polynomial_space.compute_value(renumber_inverse[i], p);
-}
-
-
-
-template <int dim>
-double
-FE_Q_Hierarchical<dim>::shape_value_component (const unsigned int i,
-                                              const Point<dim> &p,
-                                              const unsigned int component) const
-{
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  Assert (component == 0, ExcIndexRange (component, 0, 1));
-  return polynomial_space.compute_value(renumber_inverse[i], p);
-}
-
-
-
-template <int dim>
-Tensor<1,dim>
-FE_Q_Hierarchical<dim>::shape_grad (const unsigned int i,
-                                   const Point<dim> &p) const
-{
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  return polynomial_space.compute_grad(renumber_inverse[i], p);
-}
-
-
-
-template <int dim>
-Tensor<1,dim>
-FE_Q_Hierarchical<dim>::shape_grad_component (const unsigned int i,
-                                             const Point<dim> &p,
-                                             const unsigned int component) const
-{
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  Assert (component == 0, ExcIndexRange (component, 0, 1));
-  return polynomial_space.compute_grad(renumber_inverse[i], p);
-}
-
-
-
-template <int dim>
-Tensor<2,dim>
-FE_Q_Hierarchical<dim>::shape_grad_grad (const unsigned int i,
-                                        const Point<dim> &p) const
-{
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  return polynomial_space.compute_grad_grad(renumber_inverse[i], p);
-}
-
-
-
-template <int dim>
-Tensor<2,dim>
-FE_Q_Hierarchical<dim>::shape_grad_grad_component (const unsigned int i,
-                                                  const Point<dim> &p,
-                                                  const unsigned int component) const
-{
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  Assert (component == 0, ExcIndexRange (component, 0, 1));
-  return polynomial_space.compute_grad_grad(renumber_inverse[i], p);
-}
-
-
 //----------------------------------------------------------------------
 // Auxiliary functions
 //----------------------------------------------------------------------
@@ -414,6 +344,9 @@ initialize_embedding_and_restriction (const std::vector<FullMatrix<double> > &do
                                      const std::vector<FullMatrix<double> > &dofs_subcell)
 {
   const unsigned int dofs_1d = 2*this->dofs_per_vertex + this->dofs_per_line;
+  const std::vector<unsigned int> &renumber=
+    this->poly_space.get_numbering();
+  
 
   for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
     {
@@ -455,16 +388,16 @@ initialize_embedding_and_restriction (const std::vector<FullMatrix<double> > &do
                unsigned int c1 = ((c==2) || (c==3)) ? 1 : 0;
 
                this->prolongation[c](j,i) = 
-                 dofs_subcell[c0](renumber_inverse[j] % dofs_1d,
-                                  renumber_inverse[i] % dofs_1d) *
-                 dofs_subcell[c1]((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d,
-                                  (renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d);
+                 dofs_subcell[c0](renumber[j] % dofs_1d,
+                                  renumber[i] % dofs_1d) *
+                 dofs_subcell[c1]((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d,
+                                  (renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d);
 
                this->restriction[c](j,i) = 
-                 dofs_cell[c0](renumber_inverse[j] % dofs_1d,
-                               renumber_inverse[i] % dofs_1d) *
-                 dofs_cell[c1]((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d,
-                               (renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d);
+                 dofs_cell[c0](renumber[j] % dofs_1d,
+                               renumber[i] % dofs_1d) *
+                 dofs_cell[c1]((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d,
+                               (renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d);
              }
            break;
          }
@@ -478,20 +411,20 @@ initialize_embedding_and_restriction (const std::vector<FullMatrix<double> > &do
                unsigned int c2 = ((c==2) || (c==3) || (c==6) || (c==7)) ? 1 : 0;
 
                this->prolongation[c](j,i) = 
-                 dofs_subcell[c0](renumber_inverse[j] % dofs_1d,
-                                  renumber_inverse[i] % dofs_1d) *
-                 dofs_subcell[c1](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d) % dofs_1d,
-                                  ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d) % dofs_1d) *
-                 dofs_subcell[c2](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d - (((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d,
-                                  ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d - (((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d);
+                 dofs_subcell[c0](renumber[j] % dofs_1d,
+                                  renumber[i] % dofs_1d) *
+                 dofs_subcell[c1](((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) % dofs_1d,
+                                  ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) % dofs_1d) *
+                 dofs_subcell[c2](((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d - (((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d,
+                                  ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d - (((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d);
 
                this->restriction[c](j,i) = 
-                 dofs_cell[c0](renumber_inverse[j] % dofs_1d,
-                               renumber_inverse[i] % dofs_1d) *
-                 dofs_cell[c1](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d) % dofs_1d,
-                               ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d) % dofs_1d) *
-                 dofs_cell[c2](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d - (((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d,
-                               ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d - (((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d);
+                 dofs_cell[c0](renumber[j] % dofs_1d,
+                               renumber[i] % dofs_1d) *
+                 dofs_cell[c1](((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) % dofs_1d,
+                               ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) % dofs_1d) *
+                 dofs_cell[c2](((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d - (((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d,
+                               ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d - (((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d);
              }
            break;
          }
@@ -512,6 +445,9 @@ void FE_Q_Hierarchical<dim>::initialize_unit_support_points ()
     n *= degree+1;
   
   this->unit_support_points.resize(n);
+
+  const std::vector<unsigned int> &index_map_inverse=
+    this->poly_space.get_numbering_inverse();
   
   Point<dim> p;
                                    // the method of numbering allows
@@ -567,7 +503,7 @@ void FE_Q_Hierarchical<dim>::initialize_unit_support_points ()
              else
                p(2) = .5;
            }
-         this->unit_support_points[renumber[k++]] = p;
+         this->unit_support_points[index_map_inverse[k++]] = p;
        };
 }
 
@@ -907,296 +843,6 @@ FE_Q_Hierarchical<1>::face_lexicographic_to_hierarchic_numbering (const unsigned
 #endif
 
 
-template <int dim>
-UpdateFlags
-FE_Q_Hierarchical<dim>::update_once (const UpdateFlags flags) const
-{
-                                  // for this kind of elements, only
-                                  // the values can be precomputed
-                                  // once and for all. set this flag
-                                  // if the values are requested at
-                                  // all
-  return (update_default | (flags & update_values));
-}
-
-
-
-template <int dim>
-UpdateFlags
-FE_Q_Hierarchical<dim>::update_each (const UpdateFlags flags) const
-{
-  UpdateFlags out = update_default;
-
-  if (flags & update_gradients)
-    out |= update_gradients | update_covariant_transformation;
-  if (flags & update_second_derivatives)
-    out |= update_second_derivatives | update_covariant_transformation;
-
-  return out;
-}
-
-
-
-//----------------------------------------------------------------------
-// Data field initialization
-//----------------------------------------------------------------------
-
-template <int dim>
-typename Mapping<dim>::InternalDataBase *
-FE_Q_Hierarchical<dim>::get_data (const UpdateFlags      update_flags,
-                                 const Mapping<dim>    &mapping,
-                                 const Quadrature<dim> &quadrature) const
-{
-                                  // generate a new data object and
-                                  // initialize some fields
-  InternalData* data = new InternalData;
-
-                                  // check what needs to be
-                                  // initialized only once and what
-                                  // on every cell/face/subface we
-                                  // visit
-  data->update_once = update_once(update_flags);
-  data->update_each = update_each(update_flags);
-  data->update_flags = data->update_once | data->update_each;
-
-  const UpdateFlags flags(data->update_flags);
-  const unsigned int n_q_points = quadrature.n_quadrature_points;
-
-                                  // some scratch arrays
-  std::vector<double> values(0);
-  std::vector<Tensor<1,dim> > grads(0);
-  std::vector<Tensor<2,dim> > grad_grads(0);
-
-                                  // initialize fields only if really
-                                  // necessary. otherwise, don't
-                                  // allocate memory
-  if (flags & update_values)
-    {
-      values.resize (this->dofs_per_cell);
-      data->shape_values.reinit (this->dofs_per_cell,
-                                n_q_points);
-    }
-
-  if (flags & update_gradients)
-    {
-      grads.resize (this->dofs_per_cell);
-      data->shape_gradients.reinit (this->dofs_per_cell,
-                                   n_q_points);
-    }
-
-                                  // if second derivatives through
-                                  // finite differencing is required,
-                                  // then initialize some objects for
-                                  // that
-  if (flags & update_second_derivatives)
-    data->initialize_2nd (this, mapping, quadrature);
-
-                                  // next already fill those fields
-                                  // of which we have information by
-                                  // now. note that the shape
-                                  // gradients are only those on the
-                                  // unit cell, and need to be
-                                  // transformed when visiting an
-                                  // actual cell
-  if (flags & (update_values | update_gradients))
-    for (unsigned int i=0; i<n_q_points; ++i)
-      {
-       polynomial_space.compute(quadrature.point(i),
-                                values, grads, grad_grads);
-       
-       if (flags & update_values)
-         for (unsigned int k=0; k<this->dofs_per_cell; ++k)
-           data->shape_values[renumber[k]][i] = values[k];
-       
-       if (flags & update_gradients)
-         for (unsigned int k=0; k<this->dofs_per_cell; ++k)
-           data->shape_gradients[renumber[k]][i] = grads[k];
-      }
-  return data;
-}
-
-
-
-
-//----------------------------------------------------------------------
-// Fill data of FEValues
-//----------------------------------------------------------------------
-
-template <int dim>
-void
-FE_Q_Hierarchical<dim>::fill_fe_values (
-  const Mapping<dim>                            &mapping,
-  const typename DoFHandler<dim>::cell_iterator &cell,
-  const Quadrature<dim>                         &quadrature,
-  typename Mapping<dim>::InternalDataBase       &mapping_data,
-  typename Mapping<dim>::InternalDataBase       &fedata,
-  FEValuesData<dim>                             &data) const
-{
-                                  // convert data object to internal
-                                  // data for this class. fails with
-                                  // an exception if that is not
-                                  // possible
-  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-  
-  const UpdateFlags flags(fe_data.current_update_flags());
-
-  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
-    {
-      if (flags & update_values)
-       for (unsigned int i=0; i<quadrature.n_quadrature_points; ++i)
-         data.shape_values(k,i) = fe_data.shape_values[k][i];
-      
-      if (flags & update_gradients)
-       {
-         Assert (data.shape_gradients[k].size() <=
-                 fe_data.shape_gradients[k].size(),
-                 ExcInternalError());    
-         mapping.transform_covariant(data.shape_gradients[k].begin(),
-                                     data.shape_gradients[k].end(),
-                                     fe_data.shape_gradients[k].begin(),
-                                     mapping_data);
-       };
-    }
-
-  if (flags & update_second_derivatives)
-    this->compute_2nd (mapping, cell,
-                       QProjector<dim>::DataSetDescriptor::cell(),
-                       mapping_data, fe_data, data);
-}
-
-
-
-template <int dim>
-void
-FE_Q_Hierarchical<dim>::fill_fe_face_values (
-  const Mapping<dim>                            &mapping,
-  const typename DoFHandler<dim>::cell_iterator &cell,
-  const unsigned int                             face,
-  const Quadrature<dim-1>                       &quadrature,
-  typename Mapping<dim>::InternalDataBase       &mapping_data,
-  typename Mapping<dim>::InternalDataBase       &fedata,
-  FEValuesData<dim>                             &data) const
-{
-                                  // convert data object to internal
-                                  // data for this class. fails with
-                                  // an exception if that is not
-                                  // possible
-  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-
-                                  // offset determines which data set
-                                  // to take (all data sets for all
-                                  // faces are stored contiguously)
-  const typename QProjector<dim>::DataSetDescriptor offset
-    = (QProjector<dim>::DataSetDescriptor::
-       face (face, cell->face_orientation(face),
-             quadrature.n_quadrature_points));
-  
-  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
-
-  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
-    {
-      if (flags & update_values)
-        for (unsigned int i=0; i<quadrature.n_quadrature_points; ++i)
-         data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
-      
-      if (flags & update_gradients)
-       {
-         Assert (data.shape_gradients[k].size() + offset <=
-                 fe_data.shape_gradients[k].size(),
-                 ExcInternalError());
-         mapping.transform_covariant(data.shape_gradients[k].begin(),
-                                     data.shape_gradients[k].end(),
-                                     fe_data.shape_gradients[k].begin()+offset,
-                                     mapping_data);
-       };
-    }
-
-  if (flags & update_second_derivatives)
-    this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
-}
-
-
-
-template <int dim>
-void
-FE_Q_Hierarchical<dim>::fill_fe_subface_values (
-  const Mapping<dim>                            &mapping,
-  const typename DoFHandler<dim>::cell_iterator &cell,
-  const unsigned int                             face,
-  const unsigned int                             subface,
-  const Quadrature<dim-1>                       &quadrature,
-  typename Mapping<dim>::InternalDataBase       &mapping_data,
-  typename Mapping<dim>::InternalDataBase       &fedata,
-  FEValuesData<dim>                             &data) const
-{
-                                  // convert data object to internal
-                                  // data for this class. fails with
-                                  // an exception if that is not
-                                  // possible
-  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-
-                                  // offset determines which data set
-                                  // to take (all data sets for all
-                                  // sub-faces are stored contiguously)
-  const typename QProjector<dim>::DataSetDescriptor offset
-    = (QProjector<dim>::DataSetDescriptor::
-       sub_face (face, subface, cell->face_orientation(face),
-                 quadrature.n_quadrature_points));
-
-  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
-
-  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
-    {
-      if (flags & update_values)
-        for (unsigned int i=0; i<quadrature.n_quadrature_points; ++i)
-         data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
-      
-      if (flags & update_gradients)
-       {
-         Assert (data.shape_gradients[k].size() + offset <=
-                 fe_data.shape_gradients[k].size(),
-                 ExcInternalError());
-         mapping.transform_covariant(data.shape_gradients[k].begin(),
-                                     data.shape_gradients[k].end(),
-                                     fe_data.shape_gradients[k].begin()+offset,
-                                     mapping_data);
-       };
-    }
-  
-  if (flags & update_second_derivatives)
-    this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
-}
-
-
-
-template <int dim>
-unsigned int
-FE_Q_Hierarchical<dim>::n_base_elements () const
-{
-  return 1;
-}
-
-
-
-template <int dim>
-const FiniteElement<dim> &
-FE_Q_Hierarchical<dim>::base_element (const unsigned int index) const
-{
-  Assert (index==0, ExcIndexRange(index, 0, 1));
-  return *this;
-}
-
-
-
-template <int dim>
-unsigned int
-FE_Q_Hierarchical<dim>::element_multiplicity (const unsigned int index) const
-{
-  Assert (index==0, ExcIndexRange(index, 0, 1));
-  return 1;
-}
-
-
 
 template <int dim>
 bool

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.