]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
initial version of hierarchical continuous fe.
authorbrian <brian@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 12 Dec 2002 16:47:00 +0000 (16:47 +0000)
committerbrian <brian@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 12 Dec 2002 16:47:00 +0000 (16:47 +0000)
git-svn-id: https://svn.dealii.org/trunk@6811 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_q_hierarchical.h [new file with mode: 0644]
deal.II/deal.II/source/fe/fe_q_hierarchical.cc [new file with mode: 0644]

diff --git a/deal.II/deal.II/include/fe/fe_q_hierarchical.h b/deal.II/deal.II/include/fe/fe_q_hierarchical.h
new file mode 100644 (file)
index 0000000..c3cdbe7
--- /dev/null
@@ -0,0 +1,799 @@
+//---------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2002 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------
+#ifndef __deal2__fe_q_hierarchical_h
+#define __deal2__fe_q_hierarchical_h
+
+#include <base/config.h>
+#include <base/polynomial.h>
+#include <base/tensor_product_polynomials.h>
+#include <fe/fe.h>
+#include <lac/full_matrix.h>
+
+template <int dim> class TensorProductPolynomials;
+template <int dim> class MappingQ;
+
+
+
+/**
+ * Implementation of Hierarchical finite elements @p{Qp} that yield the
+ * finite element space of continuous, piecewise polynomials of degree
+ * @p{p}. This class is realized using tensor product polynomials
+ * based on a hierarchical basis of the interval @p{[-1,1]} suitable
+ * for building an @p{hp} tensor product finite element. There are not
+ * many differences between @p{FE_Q_Hierarchical} and @p{FE_Q}, except 
+ * that we now allow the degree to be nonconstant for @p{p}-refinement.
+ *
+ * The constructor of this class takes the degree @p{p} of this finite
+ * element.
+ *
+ * @sect3{Implementation}
+ *
+ * The constructor creates a @ref{TensorProductPolynomials} object
+ * that includes the tensor product of @p{Hierarchical}
+ * polynomials of degree @p{p}. This @p{TensorProductPolynomials}
+ * object provides all values and derivatives of the shape functions.
+ *
+ * Furthermore the constructor filles the @p{interface_constraints},
+ * the @p{prolongation} (embedding) and the @p{restriction}
+ * matrices. These are implemented only up to a certain degree, that
+ * is listed in the following: (fix this eventually......)
+ *
+ * @begin{itemize}
+ * @item @p{dim==1}
+ *   @begin{itemize}
+ *   @item the @p{interface_constraints} are not needed
+ *   @item the @p{prolongation} matrices up to degree 4, and
+ *   @item the @p{restriction} matrices up to degree 4.
+ *   @end{itemize}
+ * @item @p{dim==2}
+ *   @begin{itemize}
+ *   @item the @p{interface_constraints} up to degree 4,
+ *   @item the @p{prolongation} matrices up to degree 3, and
+ *   @item the @p{restriction} matrices up to degree 4.
+ *   @end{itemize}
+ * @item @p{dim==3}
+ *   @begin{itemize}
+ *   @item the @p{interface_constraints} up to degree 2,
+ *   @item the @p{prolongation} matrices up to degree 2, and
+ *   @item the @p{restriction} matrices up to degree 4.
+ *   @end{itemize}
+ * @end{itemize}
+ *
+ * @sect3{Numbering of the degrees of freedom (DoFs)}
+ *
+ * The original ordering of the shape functions represented by the
+ * @ref{TensorProductPolynomials} is a tensor product
+ * numbering. However, the shape functions on a cell are renumbered
+ * beginning with the shape functions whose support points are at the
+ * vertices, then on the line, on the quads, and finally (for 3d) on
+ * the hexes. To be explicit, these numberings are listed in the
+ * following: (support points for @p{hp}??)
+ *
+ * @sect4{Q1 elements}
+ * @begin{itemize}
+ * @item 1D case:
+ *   @begin{verbatim}
+ *      0-------1
+ *   @end{verbatim}
+ *
+ * @item 2D case:
+ *   @begin{verbatim}
+ *      3-------2
+ *      |       |
+ *      |       |
+ *      |       |
+ *      0-------1
+ *   @end{verbatim}
+ *
+ * @item 3D case:
+ *   @begin{verbatim}
+ *         7-------6        7-------6
+ *        /|       |       /       /|
+ *       / |       |      /       / |
+ *      /  |       |     /       /  |
+ *     3   |       |    3-------2   |
+ *     |   4-------5    |       |   5
+ *     |  /       /     |       |  /
+ *     | /       /      |       | /
+ *     |/       /       |       |/
+ *     0-------1        0-------1
+ *
+ *   The respective coordinate values of the support points of the degrees
+ *   of freedom are as follows:
+ *   @begin{itemize}
+ *   @item Index 0: @p{[0, 0, 0]};
+ *   @item Index 1: @p{[1, 0, 0]};
+ *   @item Index 2: @p{[1, 0, 1]};
+ *   @item Index 3: @p{[0, 0, 1]};
+ *   @item Index 4: @p{[0, 1, 0]};
+ *   @item Index 5: @p{[1, 1, 0]};
+ *   @item Index 6: @p{[1, 1, 1]};
+ *   @item Index 7: @p{[0, 1, 1]};
+ *   @end{itemize}
+ * @end{itemize}
+ * @sect4{Q2 elements}
+ * @begin{itemize}
+ * @item 1D case:
+ *   @begin{verbatim}
+ *      0---2---1
+ *   @end{verbatim}
+ *
+ * @item 2D case:
+ *   @begin{verbatim}
+ *      3---6---2
+ *      |       |
+ *      7   8   5
+ *      |       |
+ *      0---4---1
+ *   @end{verbatim}
+ *
+ * @item 3D case:
+ *   @begin{verbatim}
+ *         7--14---6        7--14---6
+ *        /|       |       /       /|
+ *      19 |       13     19      1813
+ *      /  15      |     /       /  |
+ *     3   |       |    3---10--2   |
+ *     |   4--12---5    |       |   5
+ *     |  /       /     |       9  /
+ *    11 16      17     11      | 17
+ *     |/       /       |       |/
+ *     0---8---1        0---8---1
+ *
+ *         *-------*        *-------*
+ *        /|       |       /       /|
+ *       / |  21   |      /  24   / |
+ *      /  |       |     /       /  |
+ *     *   |       |    *-------*   |
+ *     |25 *-------*    |       |23 *
+ *     |  /       /     |   20  |  /
+ *     | /  22   /      |       | /
+ *     |/       /       |       |/
+ *     *-------*        *-------* 
+ *   @end{verbatim}
+ *   The center vertex has number 26.
+ *
+ *   The respective coordinate values of the support points of the degrees
+ *   of freedom are as follows:
+ *   @begin{itemize}
+ *   @item Index 0: @p{[0, 0, 0]};
+ *   @item Index 1: @p{[1, 0, 0]};
+ *   @item Index 2: @p{[1, 0, 1]};
+ *   @item Index 3: @p{[0, 0, 1]};
+ *   @item Index 4: @p{[0, 1, 0]};
+ *   @item Index 5: @p{[1, 1, 0]};
+ *   @item Index 6: @p{[1, 1, 1]};
+ *   @item Index 7: @p{[0, 1, 1]};
+ *   @item Index 8: @p{[1/2, 0, 0]};
+ *   @item Index 9: @p{[1, 0, 1/2]};
+ *   @item Index 10: @p{[1/2, 0, 1]};
+ *   @item Index 11: @p{[0, 0, 1/2]};
+ *   @item Index 12: @p{[1/2, 1, 0]};
+ *   @item Index 13: @p{[1, 1, 1/2]};
+ *   @item Index 14: @p{[1/2, 1, 1]};
+ *   @item Index 15: @p{[0, 1, 1/2]};
+ *   @item Index 16: @p{[0, 1/2, 0]};
+ *   @item Index 17: @p{[1, 1/2, 0]};
+ *   @item Index 18: @p{[1, 1/2, 1]};
+ *   @item Index 19: @p{[0, 1/2, 1]};
+ *   @item Index 20: @p{[1/2, 0, 1/2]};
+ *   @item Index 21: @p{[1/2, 1, 1/2]};
+ *   @item Index 22: @p{[1/2, 1/2, 0]};
+ *   @item Index 23: @p{[1, 1/2, 1/2]};
+ *   @item Index 24: @p{[1/2, 1/2, 1]};
+ *   @item Index 25: @p{[0, 1/2, 1/2]};
+ *   @item Index 26: @p{[1/2, 1/2, 1/2]}; 
+ *   @end{itemize}
+ * @end{itemize}
+ * @sect4{Q3 elements}
+ * @begin{itemize}
+ * @item 1D case:
+ *   @begin{verbatim}
+ *      0--2--3--1
+ *   @end{verbatim}
+ *
+ * @item 2D case:
+ *   @begin{verbatim}
+ *      3--8--9--2
+ *      |        |
+ *      11 14 15 7
+ *      |        |
+ *      10 12 13 6
+ *      |        |
+ *      0--4--5--1
+ *   @end{verbatim}
+ *   Note the reverse ordering of degrees of freedom on the left and
+ *   upper line.
+ * @end{itemize}
+ * @sect4{Q4 elements}
+ * @begin{itemize}
+ * @item 1D case:
+ *   @begin{verbatim}
+ *      0--2--3--4--1
+ *   @end{verbatim}
+ *
+ * @item 2D case:
+ *   @begin{verbatim}
+ *      3--10-11-12-2
+ *      |           |
+ *      15 22 23 24 9
+ *      |           |
+ *      14 19 20 21 8
+ *      |           |
+ *      13 16 17 18 7
+ *      |           |
+ *      0--4--5--6--1
+ *   @end{verbatim}
+ * @end{itemize}
+ * Note the reverse ordering of degrees of freedom on the left and upper
+ * line.
+ *
+ * @author Brian Carnes, 2002
+ */
+template <int dim>
+class FE_Q_Hierarchical : public FiniteElement<dim>
+{
+  public:
+                                    /**
+                                     * Constructor for tensor product
+                                     * polynomials of degree @p{p}.
+                                     */
+    FE_Q_Hierarchical (const unsigned int p);
+    
+                                    /**
+                                     * 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,
+                                     * i.e. the value passed to the
+                                     * constructor.
+                                     */
+    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.
+                                     *
+                                     * This function returns
+                                     * @p{true}, if the shape
+                                     * function @p{shape_index} has
+                                     * non-zero values on the face
+                                     * @p{face_index}.
+                                     *
+                                     * Implementation of the
+                                     * interface in
+                                     * @ref{FiniteElement}
+                                     */
+    virtual bool has_support_on_face (const unsigned int shape_index,
+                                     const unsigned int face_index) const;
+
+                                    /**
+                                     * Determine an estimate for the
+                                     * memory consumption (in bytes)
+                                     * of this object.
+                                     *
+                                     * This function is made virtual,
+                                     * since finite element objects
+                                     * are usually accessed through
+                                     * pointers to their base class,
+                                     * rather than the class itself.
+                                     */
+    virtual unsigned int memory_consumption () const;
+
+                                     /**
+                                     * For a finite element of degree
+                                     * @p{sub_degree} < @p{degree}, we 
+                                     * return a vector which maps the 
+                                     * numbering on an FE
+                                     * of degree @p{sub_degree} into the 
+                                     * numbering on this element.
+                                     */
+    std::vector<unsigned int> get_embedding_dofs (const unsigned int sub_degree) const;
+
+  protected:    
+                                    /**
+                                     * @p{clone} function instead of
+                                     * a copy constructor.
+                                     *
+                                     * This function is needed by the
+                                     * 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:
+
+                                    /**
+                                     * Only for internal use. Its
+                                     * full name is
+                                     * @p{get_dofs_per_object_vector}
+                                     * function and it creates the
+                                     * @p{dofs_per_object} vector that is
+                                     * needed within the constructor to
+                                     * be passed to the constructor of
+                                     * @p{FiniteElementData}.
+                                     */
+    static std::vector<unsigned int> get_dpo_vector(const unsigned int degree);
+    
+                                    /**
+                                     * Map tensor product data to
+                                     * shape function numbering. This
+                                     * function is actually an alike
+                                     * replica of the respective
+                                     * function in the @ref{FETools}
+                                     * class, but is kept for three
+                                     * reasons:
+                                     *
+                                     * 1. It only operates on a
+                                     * @ref{FiniteElementData}
+                                     * structure. This is ok in the
+                                     * present context, since we can
+                                     * control which types of
+                                     * arguments it is called with
+                                     * because this is a private
+                                     * function. However, the
+                                     * publicly visible function in
+                                     * the @ref{FETools} class needs
+                                     * to make sure that the
+                                     * @ref{FiniteElementData} object
+                                     * it works on actually
+                                     * represents a continuous finite
+                                     * element, which we found too
+                                     * difficult if we do not pass an
+                                     * object of type @ref{FE_Q}
+                                     * directly.
+                                     *
+                                     * 2. If we would call the
+                                     * publicly available version of
+                                     * this function instead of this
+                                     * one, we would have to pass a
+                                     * finite element
+                                     * object. However, since the
+                                     * construction of an entire
+                                     * finite element object can be
+                                     * costly, we rather chose to
+                                     * retain this function.
+                                     *
+                                     * 3. Third reason is that we
+                                     * want to call this function for
+                                     * faces as well, by just calling
+                                     * this function for the finite
+                                     * element of one dimension
+                                     * less. If we would call the
+                                     * global function instead, this
+                                     * would require us to construct
+                                     * a second finite element object
+                                     * of one dimension less, just to
+                                     * call this function. Since that
+                                     * function does not make use of
+                                     * hanging nodes constraints,
+                                     * interpolation and restriction
+                                     * matrices, etc, this would have
+                                     * been a waste. Furthermore, it
+                                     * would have posed problems with
+                                     * template instantiations.
+                                     *
+                                     * To sum up, the existence of
+                                     * this function is a compromise
+                                     * between simplicity and proper
+                                     * library design, where we have
+                                     * chosen to weigh the simplicity
+                                     * aspect a little more than
+                                     * proper design.
+                                     */
+    static
+    void
+    lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe_data,
+                                          const unsigned int            degree,
+                                          std::vector<unsigned int>    &numbering);
+
+                                    /**
+                                     * This is an analogon to the
+                                     * previous function, but working
+                                     * on faces.
+                                     */
+    static
+    void
+    face_lexicographic_to_hierarchic_numbering (const unsigned int         degree,
+                                               std::vector<unsigned int> &numbering);
+
+
+  // not sure if needed
+                                    /**
+                                     * Initialize the
+                                     * @p{unit_support_points} field
+                                     * of the @ref{FiniteElementBase}
+                                     * class. Called from the
+                                     * constructor.
+                                     */
+    void initialize_unit_support_points ();
+
+  // not sure if needed
+                                    /**
+                                     * Initialize the
+                                     * @p{unit_face_support_points} field
+                                     * of the @ref{FiniteElementBase}
+                                     * class. Called from the
+                                     * constructor.
+                                     */
+    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.
+                                     */
+    std::vector<unsigned int> renumber;
+
+                                    /**
+                                     * Inverse renumber
+                                     * vector. i.e. mapping from
+                                     * shape function numbering to
+                                     * lexicographic numbering.
+                                     */
+    std::vector<unsigned int> renumber_inverse;
+             
+                                    /**
+                                     * Mapping from lexicographic to
+                                     * shape function numbering on first face.
+                                     */
+    std::vector<unsigned int> face_renumber;
+
+                                     /**
+                                      * The matrix @p{dofs_cell} contains the 
+                                     * values of the linear functionals of 
+                                      * the master 1d cell applied to the 
+                                     * shape functions of the two 1d subcells.
+                                     * The matrix @p{dofs_subcell} constains
+                                     * the values of the linear functionals 
+                                     * on each 1d subcell applied to the 
+                                     * shape functions on the master 1d 
+                                     * subcell. 
+                                     * We use @p{dofs_cell} and 
+                                     * @p{dofs_subcell} to compute the 
+                                     * @p{prolongation}, @p{restriction} and 
+                                     * @p{interface_constraints} matrices 
+                                     * for all dimensions.
+                                     */
+    std::vector<FullMatrix<double> > dofs_cell;
+    std::vector<FullMatrix<double> > dofs_subcell;
+
+                                    /**
+                                     * 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
+                                     * dimensions. We need this since
+                                     * we want to call the functions
+                                     * @p{get_dpo_vector} and
+                                     * @p{lexicographic_to_hierarchic_numbering}
+                                     * for the faces of the finite
+                                     * element of dimension dim+1.
+                                     */
+    template <int dim1> friend class FE_Q_Hierarchical;
+};
+
+
+/* -------------- declaration of explicit specializations ------------- */
+
+template <> void FE_Q_Hierarchical<1>::initialize_unit_face_support_points ();
+template <> void FE_Q_Hierarchical<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int,
+                                                                                  std::vector<unsigned int>&);
+
+#endif
diff --git a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc
new file mode 100644 (file)
index 0000000..bc6a907
--- /dev/null
@@ -0,0 +1,1338 @@
+//----------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2002 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------------------------------------------
+
+#include <base/quadrature.h>
+#include <base/polynomial.h>
+#include <base/tensor_product_polynomials.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe.h>
+#include <fe/mapping.h>
+#include <fe/fe_q_hierarchical.h>
+#include <fe/fe_values.h>
+
+
+
+template <int dim>
+FE_Q_Hierarchical<dim>::FE_Q_Hierarchical (const unsigned int degree)
+               :
+               FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree),1),
+                                   std::vector<bool> (FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
+                                                       false),
+                                   std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
+                                                                   std::vector<bool>(1,true))),
+               degree(degree),
+               renumber(this->dofs_per_cell, 0),
+               renumber_inverse(this->dofs_per_cell, 0),
+               face_renumber(this->dofs_per_face, 0),
+               polynomial_space(Polynomials::Hierarchical<double>::generate_complete_basis(degree))
+{
+                                  // do some internal book-keeping on
+                                  // cells and faces. if in 1d, the
+                                  // face function is empty
+  lexicographic_to_hierarchic_numbering (*this, degree, renumber);
+  face_lexicographic_to_hierarchic_numbering (degree, face_renumber);
+  
+                                  // build inverse of renumbering
+                                  // vector
+  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+    renumber_inverse[renumber[i]]=i;
+
+                                  // build the dofs_subcell, dofs_cell 
+                                  // matrices for use in building prolongation,                                   // restriction, and constraint matrices.
+  const unsigned int dofs_1d = 2*this->dofs_per_vertex + this->dofs_per_line;
+
+  for (unsigned int c=0; c<GeometryInfo<1>::children_per_cell; ++c)
+  {
+    dofs_cell.push_back (FullMatrix<double> (dofs_1d,dofs_1d) );
+    dofs_subcell.push_back (FullMatrix<double> (dofs_1d,dofs_1d) );
+
+    for (unsigned int j=0; j<dofs_1d; ++j)
+    {
+      for (unsigned int k=0; k<dofs_1d; ++k)
+      {        
+                                 // upper diagonal block
+       if ((j<=1) && (k<=1))
+       {
+         if (((c==0) && (j==0) && (k==0)) || 
+             ((c==1) && (j==1) && (k==1)))
+           dofs_cell[c](j,k) = 1.;
+         else 
+           dofs_cell[c](j,k) = 0.;
+         
+         if      (((c==0) && (j==1)) || ((c==1) && (j==0)))
+           dofs_subcell[c](j,k) = .5;
+         else if (((c==0) && (k==0)) || ((c==1) && (k==1)))
+           dofs_subcell[c](j,k) = 1.;
+         else
+           dofs_subcell[c](j,k) = 0.;
+       }
+                                // upper right block
+       else if ((j<=1) && (k>=2))
+       {
+         if (((c==0) && (j==1) && ((k % 2)==0)) ||
+             ((c==1) && (j==0) && ((k % 2)==0)))
+           dofs_subcell[c](j,k) = -1.;
+       }
+                               // lower diagonal block
+       else if ((j>=2) && (k>=2) && (j<=k))
+       {
+         double factor = 1.;
+         for (unsigned int i=1; i<=j;++i)
+           factor *= ((double) (k-i+1))/((double) i);
+         if (c==0)
+         {
+           dofs_subcell[c](j,k) = ((k+j) % 2 == 0) ? 
+             pow(.5,k)*factor: -pow(.5,k)*factor;
+           dofs_cell[c](j,k) = pow(2.,j)*factor;
+         }
+         else
+         {
+           dofs_subcell[c](j,k) = pow(.5,k)*factor;
+           dofs_cell[c](j,k) = ((k+j) % 2 == 0) ? 
+             pow(2.,j)*factor : -pow(2.,j)*factor;
+         }
+       }
+      }
+    }
+  }
+                                  // fill constraint matrices
+  if (dim==2 || dim==3)
+  {
+    this->interface_constraints.reinit ( (dim==2) ? 1 + 2*(degree-1) : 
+                         5 + 12*(degree-1) + 4*(degree-1)*(degree-1),
+                                        (dim==2) ? (degree+1) : 
+                                        (degree+1)*(degree+1) );
+    switch (dim)
+    {
+      case 2:
+       // vertex node
+       for (unsigned int i=0; i<dofs_1d; ++i)
+         interface_constraints(0,i) = dofs_subcell[0](1,i); 
+       // edge nodes
+       for (unsigned int c=0; c<GeometryInfo<1>::children_per_cell; ++c)
+         for (unsigned int i=0; i<dofs_1d; ++i)
+           for (unsigned int j=2; j<dofs_1d; ++j)
+             interface_constraints(1 + c*(degree-1) + j - 2,i) = 
+                                       dofs_subcell[c](j,i);
+       break;
+      case 3:
+       for (unsigned int i=0; i<dofs_1d * dofs_1d; i++)
+       {
+         // center vertex node   
+         interface_constraints(0,face_renumber[i]) = 
+           dofs_subcell[0](1,i % dofs_1d) * 
+           dofs_subcell[0](1,(i - (i % dofs_1d)) / dofs_1d);
+
+         // boundary vertex nodes
+         interface_constraints(1,face_renumber[i]) = 
+           dofs_subcell[0](1, i % dofs_1d) * 
+           dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
+         interface_constraints(2,face_renumber[i]) = 
+           dofs_subcell[1](1, i % dofs_1d) * 
+           dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
+         interface_constraints(3,face_renumber[i]) = 
+           dofs_subcell[1](0, i % dofs_1d) * 
+           dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
+         interface_constraints(4,face_renumber[i]) = 
+           dofs_subcell[0](0, i % dofs_1d) * 
+           dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
+         
+         // interior edges
+         for (unsigned int j=0; j<degree-1; j++)
+         {
+           interface_constraints(5 + j,face_renumber[i]) = 
+             dofs_subcell[0](1, i % dofs_1d) * 
+             dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
+           interface_constraints(5 + (degree-1) + j,face_renumber[i]) = 
+             dofs_subcell[1](2 + j, i % dofs_1d) * 
+             dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
+           interface_constraints(5 + 2*(degree-1) + j,face_renumber[i]) = 
+             dofs_subcell[0](1,i % dofs_1d) * 
+             dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
+           interface_constraints(5 + 3*(degree-1) + j,face_renumber[i]) = 
+             dofs_subcell[0](2 + j,i % dofs_1d) * 
+             dofs_subcell[1](0, (i - (i % dofs_1d)) / dofs_1d);
+         }
+
+         // boundary edges
+         for (unsigned int j=0; j<degree-1; j++)
+         {
+           // bottom edge 
+           interface_constraints(5 + 4*(degree-1) + j,face_renumber[i]) =
+             dofs_subcell[0](2 + j, i % dofs_1d) * 
+             dofs_subcell[0](0,     (i - (i % dofs_1d)) / dofs_1d);
+           interface_constraints(5 + 4*(degree-1) + (degree-1) + j,face_renumber[i]) =
+             dofs_subcell[1](2 + j, i % dofs_1d) * 
+             dofs_subcell[0](0,     (i - (i % dofs_1d)) / dofs_1d);
+           // right edge
+           interface_constraints(5 + 4*(degree-1) + 2*(degree-1) + j,face_renumber[i]) =
+             dofs_subcell[1](1,     i % dofs_1d) * 
+             dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
+           interface_constraints(5 + 4*(degree-1) + 3*(degree-1) + j,face_renumber[i]) =
+             dofs_subcell[1](1,     i % dofs_1d) * 
+             dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
+           // top edge
+           interface_constraints(5 + 4*(degree-1) + 4*(degree-1) + j,face_renumber[i]) =
+             dofs_subcell[0](2 + j, i % dofs_1d) * 
+             dofs_subcell[1](1,     (i - (i % dofs_1d)) / dofs_1d);
+           interface_constraints(5 + 4*(degree-1) + 5*(degree-1) + j,face_renumber[i]) =
+             dofs_subcell[1](2 + j, i % dofs_1d) * 
+             dofs_subcell[1](1,     (i - (i % dofs_1d)) / dofs_1d);
+           // left edge
+           interface_constraints(5 + 4*(degree-1) + 6*(degree-1) + j,face_renumber[i]) =
+             dofs_subcell[0](0,     i % dofs_1d) * 
+             dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
+           interface_constraints(5 + 4*(degree-1) + 7*(degree-1) + j,face_renumber[i]) =
+             dofs_subcell[0](0,     i % dofs_1d) * 
+             dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
+         }
+
+         // interior faces
+         for (unsigned int j=0; j<degree-1; j++)
+           for (unsigned int k=0; k<degree-1; k++)
+           {
+             // subcell 0
+             interface_constraints(5 + 12*(degree-1) + j + k*(degree-1),face_renumber[i]) =
+               dofs_subcell[0](2 + j, i % dofs_1d) * 
+               dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
+             // subcell 1
+             interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + (degree-1)*(degree-1),face_renumber[i]) =
+               dofs_subcell[1](2 + j, i % dofs_1d) * 
+               dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
+             // subcell 2
+             interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + 2*(degree-1)*(degree-1),face_renumber[i]) =
+               dofs_subcell[1](2 + j, i % dofs_1d) * 
+               dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
+             // subcell 3
+             interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + 3*(degree-1)*(degree-1),face_renumber[i]) =
+               dofs_subcell[0](2 + j, i % dofs_1d) * 
+               dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
+           }
+       }
+       break;
+    }
+  };
+
+                                  // fill prolongation and restriction 
+                                   // matrices
+  if (dim==1)
+  {
+    for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
+    {
+      this->prolongation[c].reinit (this->dofs_per_cell,this->dofs_per_cell);
+      this->prolongation[c].fill (dofs_subcell[c]);
+
+      this->restriction[c].reinit (this->dofs_per_cell,this->dofs_per_cell);
+      this->restriction[c].fill (dofs_cell[c]);
+    }
+  }
+  else if (dim==2 || dim==3)
+  {
+    for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
+    {
+      this->prolongation[c].reinit (this->dofs_per_cell,this->dofs_per_cell);
+      this->restriction[c].reinit (this->dofs_per_cell,this->dofs_per_cell);
+    }
+                                       // j loops over dofs in the subcell. 
+                                       // These are the rows in the 
+                                       // embedding matrix.
+    for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+    {
+                                       // i loops over the dofs in the master 
+                                       // cell. These are the columns in 
+                                       // the embedding matrix.
+      for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+      {
+       switch (dim)
+       {
+         case 2:
+           for (unsigned int c=0; c<GeometryInfo<2>::children_per_cell; ++c)
+           {
+             unsigned int c0 = ((c==1) || (c==2)) ? 1 : 0;
+             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);
+
+             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);
+           }
+           break;
+
+         case 3:
+           for (unsigned int c=0; c<GeometryInfo<3>::children_per_cell; ++c)
+           {
+             unsigned int c0 = ((c==1) || (c==2) || (c==5) || (c==6)) ? 1 : 0;
+             unsigned int c1 = ((c==4) || (c==5) || (c==6) || (c==7)) ? 1 : 0;
+              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);
+
+             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);
+           }
+           break;
+       }
+      }
+    }
+  }
+  else
+    Assert (false, ExcNotImplemented());
+                                  // finally fill in support points
+                                  // on cell and face
+  initialize_unit_support_points ();
+  initialize_unit_face_support_points ();
+};
+
+
+
+template <int dim>
+FiniteElement<dim> *
+FE_Q_Hierarchical<dim>::clone() const
+{
+  return new FE_Q_Hierarchical<dim>(degree);
+}
+
+
+
+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
+//----------------------------------------------------------------------
+
+
+
+template <int dim>
+void FE_Q_Hierarchical<dim>::initialize_unit_support_points ()
+{
+                                  // number of points: (degree+1)^dim
+  unsigned int n = degree+1;
+  for (unsigned int i=1; i<dim; ++i)
+    n *= degree+1;
+  
+  unit_support_points.resize(n);
+  
+  Point<dim> p;
+                                   // the method of numbering allows
+                                   // each dof to be associated with a
+                                   // support point. There is
+                                   // only one support point per
+                                   // vertex, line, quad, hex, etc.
+  unsigned int k=0;
+  for (unsigned int iz=0; iz <= ((dim>2) ? degree : 0) ; ++iz)
+    for (unsigned int iy=0; iy <= ((dim>1) ? degree : 0) ; ++iy)
+      for (unsigned int ix=0; ix<=degree; ++ix)
+      {
+       if (ix==0)
+         p(0) =  0.;
+       else if (ix==1)
+         p(0) =  1.;
+       else
+         p(0) = .5;
+       if (dim>1)
+       {
+         if (iy==0)
+           p(1) =  0.;
+         else if (iy==1)
+           p(1) =  1.;
+         else
+           p(1) = .5;
+       }
+       if (dim>2)
+       {
+         if (iz==0)
+           p(2) =  0.;
+         else if (iz==1)
+           p(2) =  1.;
+         else
+           p(2) = .5;
+       }
+       unit_support_points[renumber[k++]] = p;
+      };
+};
+
+
+#if deal_II_dimension == 1
+
+template <>
+void FE_Q_Hierarchical<1>::initialize_unit_face_support_points ()
+{
+                                  // no faces in 1d, so nothing to do
+};
+
+#endif
+
+
+template <int dim>
+void FE_Q_Hierarchical<dim>::initialize_unit_face_support_points ()
+{
+  const unsigned int codim = dim-1;
+  
+                                  // number of points: (degree+1)^codim
+  unsigned int n = degree+1;
+  for (unsigned int i=1; i<codim; ++i)
+    n *= degree+1;
+  
+  unit_face_support_points.resize(n);
+  
+  Point<codim> p;
+  
+  unsigned int k=0;
+  for (unsigned int iz=0; iz <= ((codim>2) ? degree : 0) ; ++iz)
+    for (unsigned int iy=0; iy <= ((codim>1) ? degree : 0) ; ++iy)
+      for (unsigned int ix=0; ix<=degree; ++ix)
+      {
+       if (ix==0)
+         p(0) =  0.;
+       else if (ix==1)
+         p(0) =  1.;
+       else
+         p(0) = .5;
+       if (codim>1)
+       {
+         if (iy==0)
+           p(1) =  0.;
+         else if (iy==1)
+           p(1) =  1.;
+         else
+           p(1) = .5;
+       }
+       if (codim>2)
+       {
+         if (iz==0)
+           p(2) =  0.;
+         else if (iz==1)
+           p(2) =  1.;
+         else
+           p(2) = .5;
+       }
+       unit_face_support_points[face_renumber[k++]] = p;
+      };
+};
+
+
+                                           // we use same dpo_vector as FE_Q
+template <int dim>
+std::vector<unsigned int>
+FE_Q_Hierarchical<dim>::get_dpo_vector(const unsigned int deg)
+{
+  std::vector<unsigned int> dpo(dim+1, static_cast<unsigned int>(1));
+  for (unsigned int i=1; i<dpo.size(); ++i)
+    dpo[i]=dpo[i-1]*(deg-1);
+  return dpo;
+}
+
+
+
+template <int dim>
+void
+FE_Q_Hierarchical<dim>::lexicographic_to_hierarchic_numbering (
+                            const FiniteElementData<dim> &fe_data,
+                             const unsigned int            degree,
+                            std::vector<unsigned int>    &renumber)
+{
+  const unsigned int n = degree+1;
+
+
+  if (degree == 0)
+  {
+    Assert ( (fe_data.dofs_per_vertex == 0) &&
+            ((fe_data.dofs_per_line == 0) || (dim == 1)) &&
+            ((fe_data.dofs_per_quad == 0) || (dim == 2)) &&
+            ((fe_data.dofs_per_hex == 0)  || (dim == 3)),
+            ExcInternalError() );
+    renumber[0] = 0;
+  };
+
+  if (degree > 0)
+  {
+    Assert (fe_data.dofs_per_vertex == 1, ExcInternalError());
+    for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+    {
+      unsigned int index = 0;
+                                          // Find indices of vertices.
+                                          // Unfortunately, somebody
+                                          // switched the upper corner
+                                          // points of a quad. The same
+                                          // person decided to find a very
+                                          // creative numbering of the
+                                          // vertices of a hexahedron.
+                                          // Therefore, this looks quite
+                                          // sophisticated.
+                                          //
+                                          // NB: This same person
+                                          // claims to have had good
+                                          // reasons then, but seems to
+                                          // have forgotten about
+                                          // them. At least, the
+                                          // numbering was discussed
+                                          // with the complaining
+                                          // person back then when all
+                                          // began :-)
+      switch (dim)
+      {
+        case 1:
+       {
+         const unsigned int values[GeometryInfo<1>::vertices_per_cell]
+           = { 0, 1 };
+         index = values[i];
+         break;
+       };
+            
+        case 2:
+       {
+         const unsigned int values[GeometryInfo<2>::vertices_per_cell]
+           = { 0, 1, n + 1, n };
+         index = values[i];
+         break;
+       };
+            
+        case 3:
+       {
+         const unsigned int values[GeometryInfo<3>::vertices_per_cell]
+           = { 0,             1,
+               n * n + 1,     n * n,
+               n,             n + 1,
+               n * n + n + 1, n * n + n};
+         index = values[i];
+         break;
+       };
+       
+        default:
+         Assert(false, ExcNotImplemented());
+      }
+      
+      renumber[index] = i;
+    }
+  };
+                                  // for degree 2 and higher: Lines,
+                                  // quads, hexes etc also carry
+                                  // degrees of freedom
+  if (degree > 1)
+  {
+    Assert (fe_data.dofs_per_line == degree-1, ExcInternalError());
+    Assert ((fe_data.dofs_per_quad == (degree-1)*(degree-1)) ||
+           (dim < 2), ExcInternalError());
+    Assert ((fe_data.dofs_per_hex == (degree-1)*(degree-1)*(degree-1)) ||
+           (dim < 3), ExcInternalError());
+    
+    for (unsigned int i=0; i<GeometryInfo<dim>::lines_per_cell; ++i)
+    {
+      unsigned int index = fe_data.first_line_index
+                      + i*fe_data.dofs_per_line;
+      unsigned int incr = 0;
+      unsigned int tensorstart = 0;
+                                          // This again looks quite
+                                          // strange because of the odd
+                                          // numbering scheme.
+      switch (i+100*dim)
+      {
+                                               // lines in x-direction
+        case 100:
+        case 200: case 202:
+        case 300: case 302: case 304: case 306:
+         incr = 1;
+         break;
+                                               // lines in y-direction
+        case 201: case 203:
+        case 308: case 309: case 310: case 311:
+         incr = n;
+         break;
+                                               // lines in z-direction
+        case 301: case 303: case 305: case 307:
+         incr = n * n;
+         break;
+        default:
+         Assert(false, ExcNotImplemented());
+      }
+      switch (i+100*dim)
+      {
+                                                // x=y=z=0
+        case 100:
+        case 200: case 203:
+        case 300: case 303: case 308:
+         tensorstart = 0;
+         break;
+                                                // x=1 y=z=0
+        case 201:
+        case 301: case 309:
+         tensorstart = 1;
+         break;
+                                                // y=1 x=z=0
+        case 202:
+        case 304: case 307:
+         tensorstart = n;
+         break;
+                                                // x=z=1 y=0
+        case 310:
+         tensorstart = n * n + 1;
+         break;
+                                                // z=1 x=y=0
+        case 302: case 311:
+         tensorstart = n * n;
+         break;
+                                                // x=y=1 z=0
+        case 305:
+         tensorstart = n + 1;
+         break;
+                                                // y=z=1 x=0
+        case 306:
+         tensorstart = n * n + n;
+         break;
+        default:
+         Assert(false, ExcNotImplemented());         
+      }
+         
+      for (unsigned int jx = 2; jx<=degree ;++jx)
+      {
+       unsigned int tensorindex = tensorstart + jx * incr;
+       renumber[tensorindex] = index++;
+      }
+    }
+    
+    for (int i=0; i<static_cast<signed int>(GeometryInfo<dim>::quads_per_cell); ++i)
+    {
+      unsigned int index = fe_data.first_quad_index+i*fe_data.dofs_per_quad;
+      unsigned int tensorstart = 0;
+      unsigned int incx = 0;
+      unsigned int incy = 0;
+      switch (i)
+      {
+                                                // z=0 (dim==2), y=0 (dim==3)
+        case 0:
+         tensorstart = 0; incx = 1;
+         if (dim==2)
+           incy = n;
+         else
+           incy = n * n;
+         break;
+                                               // y=1
+        case 1:
+         tensorstart = n; incx = 1; incy = n * n;
+         break;
+                                               // z=0
+        case 2:
+         tensorstart = 0; incx = 1; incy = n;
+         break;
+                                               // x=1
+        case 3:
+         tensorstart = 1; incx = n; incy = n * n;
+         break;
+                                               // z=1
+        case 4:
+         tensorstart = n * n; incx = 1; incy = n;
+         break;
+                                               // x=0
+        case 5:
+         tensorstart = 0; incx = n; incy = n * n;
+         break;
+        default:
+         Assert(false, ExcNotImplemented());         
+      }
+         
+      for (unsigned int jy = 2; jy<=degree; jy++)
+       for (unsigned int jx = 2; jx<=degree ;++jx)
+       {
+         unsigned int tensorindex = tensorstart
+                                  + jx * incx + jy * incy;
+         renumber[tensorindex] = index++;
+       }
+    }
+    
+    if (GeometryInfo<dim>::hexes_per_cell > 0)
+      for (int i=0; i<static_cast<signed int>(GeometryInfo<dim>::hexes_per_cell); ++i)
+      {
+       unsigned int index = fe_data.first_hex_index;
+           
+       for (unsigned int jz = 2; jz<=degree; jz++)
+         for (unsigned int jy = 2; jy<=degree; jy++)
+           for (unsigned int jx = 2; jx<=degree; jx++)
+           {
+             const unsigned int tensorindex = jx + jy * n + jz * n * n;
+             renumber[tensorindex]=index++;
+           }  
+      } 
+  }
+}
+
+
+
+template <int dim>
+void
+FE_Q_Hierarchical<dim>::face_lexicographic_to_hierarchic_numbering (
+                        const unsigned int         degree,
+                        std::vector<unsigned int> &numbering)
+{
+  FiniteElementData<dim-1> fe_data(FE_Q_Hierarchical<dim-1>::get_dpo_vector(degree),1);
+  FE_Q_Hierarchical<dim-1>::lexicographic_to_hierarchic_numbering (fe_data,
+                                                       degree,
+                                                       numbering); 
+}
+
+
+#if (deal_II_dimension == 1)
+
+template <>
+void
+FE_Q_Hierarchical<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int,
+                                                std::vector<unsigned int>&)
+{}
+
+#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)
+    compute_2nd (mapping, cell, 0, 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 unsigned int offset = 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)
+    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 unsigned int offset = ((face * GeometryInfo<dim>::subfaces_per_face + subface)
+                               * 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)
+    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
+FE_Q_Hierarchical<dim>::has_support_on_face (const unsigned int shape_index,
+                                            const unsigned int face_index) const
+{
+  Assert (shape_index < this->dofs_per_cell,
+         ExcIndexRange (shape_index, 0, this->dofs_per_cell));
+  Assert (face_index < GeometryInfo<dim>::faces_per_cell,
+         ExcIndexRange (face_index, 0, GeometryInfo<dim>::faces_per_cell));
+
+
+                                  // in 1d, things are simple. since
+                                  // there is only one degree of
+                                  // freedom per vertex in this
+                                  // class, the first is on vertex 0
+                                  // (==face 0 in some sense), the
+                                  // second on face 1:
+  if (dim==1)
+    return (((shape_index == 0) && (face_index == 0)) ||
+           ((shape_index == 1) && (face_index == 1)));
+  else
+                                    // more dimensions
+    {
+                                       // first, special-case interior
+                                       // shape functions, since they
+                                       // have no support no-where on
+                                       // the boundary
+      if (((dim==2) && (shape_index>=first_quad_index))
+          ||
+          ((dim==3) && (shape_index>=first_hex_index)))
+        return false;
+                                       
+                                       // let's see whether this is a
+                                       // vertex
+      if (shape_index < first_line_index) 
+        {
+                                           // for Q elements, there is
+                                           // one dof per vertex, so
+                                           // shape_index==vertex_number. check
+                                           // whether this vertex is
+                                           // on the given face. thus,
+                                           // for each face, give a
+                                           // list of vertices
+          const unsigned int vertex_no = shape_index;
+          Assert (vertex_no < GeometryInfo<dim>::vertices_per_cell,
+                  ExcInternalError());
+          switch (dim)
+            {
+              case 2:
+              {
+                static const unsigned int face_vertices[4][2] =
+                  { {0,1},{1,2},{2,3},{0,3} };
+                return ((face_vertices[face_index][0] == vertex_no)
+                        ||
+                        (face_vertices[face_index][1] == vertex_no));
+              };
+
+              case 3:
+              {
+                static const unsigned int face_vertices[6][4] =
+                  { {0,1,2,3},{4,5,6,7},{0,1,5,4},
+                    {1,5,6,2},{3,2,6,7},{0,4,7,3} };
+                return ((face_vertices[face_index][0] == vertex_no)||
+                        (face_vertices[face_index][1] == vertex_no)||
+                        (face_vertices[face_index][2] == vertex_no)||
+                        (face_vertices[face_index][3] == vertex_no));
+              };
+
+              default:
+                    Assert (false, ExcNotImplemented());
+            };
+        }
+      else if (shape_index < first_quad_index)
+                                         // ok, dof is on a line
+        {
+          const unsigned int line_index
+            = (shape_index - first_line_index) / dofs_per_line;
+          Assert (line_index < GeometryInfo<dim>::lines_per_cell,
+                  ExcInternalError());
+
+                                           // in 2d, the line is the
+                                           // face, so get the line
+                                           // index
+          if (dim == 2)
+            return (line_index == face_index);
+          else if (dim == 3)
+            {
+                                               // see whether the
+                                               // given line is on the
+                                               // given face. use
+                                               // table technique
+                                               // again
+              static const unsigned int face_lines[6][4] =
+                { {0,1,2,3},{4,5,6,7},{0,8,9,4},
+                  {1,9,5,10},{2,10,6,11},{3,8,7,11} };
+              return ((face_lines[face_index][0] == line_index)||
+                      (face_lines[face_index][1] == line_index)||
+                      (face_lines[face_index][2] == line_index)||
+                      (face_lines[face_index][3] == line_index));
+            }
+          else
+            Assert (false, ExcNotImplemented());
+        }
+      else if (shape_index < first_hex_index)
+                                         // dof is on a quad
+        {
+          const unsigned int quad_index 
+            = (shape_index - first_quad_index) / dofs_per_quad;
+          Assert (quad_index < GeometryInfo<dim>::quads_per_cell,
+                  ExcInternalError());
+          
+                                           // in 2d, cell bubble are
+                                           // zero on all faces. but
+                                           // we have treated this
+                                           // case above already
+          Assert (dim != 2, ExcInternalError());
+
+                                           // in 3d,
+                                           // quad_index=face_index
+          if (dim == 3)
+            return (quad_index == face_index);
+          else
+            Assert (false, ExcNotImplemented());
+        }
+      else
+                                         // dof on hex
+        {
+                                           // can only happen in 3d,
+                                           // but this case has
+                                           // already been covered
+                                           // above
+          Assert (false, ExcNotImplemented());
+          return false;
+        };
+    };
+
+                                   // we should not have gotten here
+  Assert (false, ExcInternalError());
+  return false;
+
+}
+
+
+
+template <int dim>
+std::vector<unsigned int> 
+FE_Q_Hierarchical<dim>::get_embedding_dofs (const unsigned int sub_degree) const
+{
+  Assert ((sub_degree>0) && (sub_degree<=degree),
+         ExcIndexRange(sub_degree, 1, degree));
+
+  if (dim==1)
+  {
+    std::vector<unsigned int> embedding_dofs (sub_degree+1);
+    for (unsigned int i=0; i<(sub_degree+1); ++i)
+      embedding_dofs[i] = i;
+
+    return embedding_dofs;
+  }
+  
+  if (sub_degree==1)
+  {
+    std::vector<unsigned int> embedding_dofs (GeometryInfo<dim>::vertices_per_cell);
+    for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+      embedding_dofs[i] = i;
+
+    return embedding_dofs;
+  }
+  else if (sub_degree==degree)
+  {
+    std::vector<unsigned int> embedding_dofs (this->dofs_per_cell);
+    for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+      embedding_dofs[i] = i;
+
+    return embedding_dofs;
+  }
+
+  if ((dim==2) || (dim==3))
+  {
+    std::vector<unsigned int> embedding_dofs ( (dim==2) ? 
+                                              (sub_degree+1) * (sub_degree+1) :
+                                              (sub_degree+1) * (sub_degree+1) * (sub_degree+1) );
+
+    for (unsigned int i=0; i<( (dim==2) ? 
+                              (sub_degree+1) * (sub_degree+1) :
+                              (sub_degree+1) * (sub_degree+1) * (sub_degree+1) ); ++i)
+    {
+      // vertex
+      if (i<GeometryInfo<dim>::vertices_per_cell)
+        embedding_dofs[i] = i;
+      // line
+      else if (i<(GeometryInfo<dim>::vertices_per_cell + 
+                 GeometryInfo<dim>::lines_per_cell * (sub_degree-1)))
+      {
+        const unsigned int j = (i - GeometryInfo<dim>::vertices_per_cell) %
+                              (sub_degree-1);
+        const unsigned int line = (i - GeometryInfo<dim>::vertices_per_cell - j) / (sub_degree-1);
+
+        embedding_dofs[i] = GeometryInfo<dim>::vertices_per_cell + 
+                           line * (degree-1) + j;
+      }
+      // quad
+      else if (i<(GeometryInfo<dim>::vertices_per_cell + 
+                 GeometryInfo<dim>::lines_per_cell * (sub_degree-1)) + 
+                 GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1))
+      {
+        const unsigned int j = (i - GeometryInfo<dim>::vertices_per_cell -
+                                GeometryInfo<dim>::lines_per_cell * (sub_degree-1)) % (sub_degree-1);
+        const unsigned int k = ( (i - GeometryInfo<dim>::vertices_per_cell -
+                                  GeometryInfo<dim>::lines_per_cell * (sub_degree-1) - j) / (sub_degree-1) ) % (sub_degree-1);
+        const unsigned int face = (i - GeometryInfo<dim>::vertices_per_cell - 
+                                   GeometryInfo<dim>::lines_per_cell * (sub_degree-1) - k * (sub_degree-1) - j) / ( (sub_degree-1) * (sub_degree-1) );
+
+        embedding_dofs[i] = GeometryInfo<dim>::vertices_per_cell + 
+                            GeometryInfo<dim>::lines_per_cell * (degree-1) +
+                           face * (degree-1) * (degree-1) +
+                            k * (degree-1) + j;
+      }
+      // hex
+      else if (i<(GeometryInfo<dim>::vertices_per_cell + 
+                 GeometryInfo<dim>::lines_per_cell * (sub_degree-1)) + 
+                 GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1) +
+                 GeometryInfo<dim>::hexes_per_cell * (sub_degree-1) * (sub_degree-1) * (sub_degree-1))
+      {
+       const unsigned int j = (i - GeometryInfo<dim>::vertices_per_cell -
+                                GeometryInfo<dim>::lines_per_cell * (sub_degree-1) -
+                                GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1) ) % (sub_degree-1);
+        const unsigned int k = ( (i - GeometryInfo<dim>::vertices_per_cell -
+                                  GeometryInfo<dim>::lines_per_cell * (sub_degree-1) -
+                                  GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1) - j) / (sub_degree-1) ) % (sub_degree-1);
+        const unsigned int l = (i - GeometryInfo<dim>::vertices_per_cell -
+                                GeometryInfo<dim>::lines_per_cell * (sub_degree-1) -
+                               GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1) - j - k * (sub_degree-1)) / ( (sub_degree-1) * (sub_degree-1) );
+        
+        embedding_dofs[i] = GeometryInfo<dim>::vertices_per_cell + 
+                            GeometryInfo<dim>::lines_per_cell * (degree-1) +
+                            GeometryInfo<dim>::quads_per_cell * (degree-1) * (degree-1) +
+                            l * (degree-1) * (degree-1) + k * (degree-1) + j;
+      }
+    }
+
+  return embedding_dofs;
+  }
+  else
+  {
+    Assert(false, ExcNotImplemented ());
+    return std::vector<unsigned int> ();
+  }
+}
+
+
+
+template <int dim>
+unsigned int
+FE_Q_Hierarchical<dim>::memory_consumption () const
+{
+  Assert (false, ExcNotImplemented ());
+  return 0;
+}
+
+
+
+template <int dim>
+unsigned int
+FE_Q_Hierarchical<dim>::get_degree () const
+{
+  return degree;
+};
+
+
+
+template class FE_Q_Hierarchical<deal_II_dimension>;

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.