From 37b4c54c496d4bcfec4315629f9d039a7a0def84 Mon Sep 17 00:00:00 2001 From: brian Date: Thu, 12 Dec 2002 16:47:00 +0000 Subject: [PATCH] initial version of hierarchical continuous fe. git-svn-id: https://svn.dealii.org/trunk@6811 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/fe/fe_q_hierarchical.h | 799 ++++++++++ .../deal.II/source/fe/fe_q_hierarchical.cc | 1338 +++++++++++++++++ 2 files changed, 2137 insertions(+) create mode 100644 deal.II/deal.II/include/fe/fe_q_hierarchical.h create mode 100644 deal.II/deal.II/source/fe/fe_q_hierarchical.cc 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 index 0000000000..c3cdbe7a8a --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_q_hierarchical.h @@ -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 +#include +#include +#include +#include + +template class TensorProductPolynomials; +template 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 +class FE_Q_Hierarchical : public FiniteElement +{ + 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 &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 &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 &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 &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 &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 &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 & 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 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 * clone() const; + + /** + * Prepare internal data + * structures and fill in values + * independent of the cell. + */ + virtual + typename Mapping::InternalDataBase * + get_data (const UpdateFlags, + const Mapping& mapping, + const Quadrature& quadrature) const ; + + /** + * Implementation of the same + * function in + * @ref{FiniteElement}. + */ + virtual void + fill_fe_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData& data) const; + + /** + * Implementation of the same + * function in + * @ref{FiniteElement}. + */ + virtual void + fill_fe_face_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData& data) const ; + + /** + * Implementation of the same + * function in + * @ref{FiniteElement}. + */ + virtual void + fill_fe_subface_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData& 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 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 &fe_data, + const unsigned int degree, + std::vector &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 &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 renumber; + + /** + * Inverse renumber + * vector. i.e. mapping from + * shape function numbering to + * lexicographic numbering. + */ + std::vector renumber_inverse; + + /** + * Mapping from lexicographic to + * shape function numbering on first face. + */ + std::vector 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 > dofs_cell; + std::vector > dofs_subcell; + + /** + * Pointer to the tensor + * product polynomials. + */ + const TensorProductPolynomials 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::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 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&); + +#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 index 0000000000..bc6a907517 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + +template +FE_Q_Hierarchical::FE_Q_Hierarchical (const unsigned int degree) + : + FiniteElement (FiniteElementData(get_dpo_vector(degree),1), + std::vector (FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, + false), + std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, + std::vector(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::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; idofs_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::children_per_cell; ++c) + { + dofs_cell.push_back (FullMatrix (dofs_1d,dofs_1d) ); + dofs_subcell.push_back (FullMatrix (dofs_1d,dofs_1d) ); + + for (unsigned int j=0; j=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::children_per_cell; ++c) + for (unsigned int i=0; i::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::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; jdofs_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; idofs_per_cell; ++i) + { + switch (dim) + { + case 2: + for (unsigned int c=0; c::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::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 +FiniteElement * +FE_Q_Hierarchical::clone() const +{ + return new FE_Q_Hierarchical(degree); +} + + + +template +double +FE_Q_Hierarchical::shape_value (const unsigned int i, + const Point &p) const +{ + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + return polynomial_space.compute_value(renumber_inverse[i], p); +} + + + +template +double +FE_Q_Hierarchical::shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (idofs_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 +Tensor<1,dim> +FE_Q_Hierarchical::shape_grad (const unsigned int i, + const Point &p) const +{ + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + return polynomial_space.compute_grad(renumber_inverse[i], p); +} + + + +template +Tensor<1,dim> +FE_Q_Hierarchical::shape_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (idofs_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 +Tensor<2,dim> +FE_Q_Hierarchical::shape_grad_grad (const unsigned int i, + const Point &p) const +{ + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + return polynomial_space.compute_grad_grad(renumber_inverse[i], p); +} + + + +template +Tensor<2,dim> +FE_Q_Hierarchical::shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (idofs_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 +void FE_Q_Hierarchical::initialize_unit_support_points () +{ + // number of points: (degree+1)^dim + unsigned int n = degree+1; + for (unsigned int i=1; i 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 +void FE_Q_Hierarchical::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 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 +std::vector +FE_Q_Hierarchical::get_dpo_vector(const unsigned int deg) +{ + std::vector dpo(dim+1, static_cast(1)); + for (unsigned int i=1; i +void +FE_Q_Hierarchical::lexicographic_to_hierarchic_numbering ( + const FiniteElementData &fe_data, + const unsigned int degree, + std::vector &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::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::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(GeometryInfo::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::hexes_per_cell > 0) + for (int i=0; i(GeometryInfo::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 +void +FE_Q_Hierarchical::face_lexicographic_to_hierarchic_numbering ( + const unsigned int degree, + std::vector &numbering) +{ + FiniteElementData fe_data(FE_Q_Hierarchical::get_dpo_vector(degree),1); + FE_Q_Hierarchical::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&) +{} + +#endif + + +template +UpdateFlags +FE_Q_Hierarchical::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 +UpdateFlags +FE_Q_Hierarchical::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 +typename Mapping::InternalDataBase * +FE_Q_Hierarchical::get_data (const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &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 values(0); + std::vector > grads(0); + std::vector > 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; idofs_per_cell; ++k) + data->shape_values[renumber[k]][i] = values[k]; + + if (flags & update_gradients) + for (unsigned int k=0; kdofs_per_cell; ++k) + data->shape_gradients[renumber[k]][i] = grads[k]; + } + return data; +} + + + + +//---------------------------------------------------------------------- +// Fill data of FEValues +//---------------------------------------------------------------------- + +template +void +FE_Q_Hierarchical::fill_fe_values ( + const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &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 (fedata); + + const UpdateFlags flags(fe_data.current_update_flags()); + + for (unsigned int k=0; kdofs_per_cell; ++k) + { + if (flags & update_values) + for (unsigned int i=0; i +void +FE_Q_Hierarchical::fill_fe_face_values ( + const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &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 (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; kdofs_per_cell; ++k) + { + if (flags & update_values) + for (unsigned int i=0; i +void +FE_Q_Hierarchical::fill_fe_subface_values ( + const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face, + const unsigned int subface, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &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 (fedata); + + // offset determines which data set + // to take (all data sets for all + // sub-faces are stored contiguously) + const unsigned int offset = ((face * GeometryInfo::subfaces_per_face + subface) + * quadrature.n_quadrature_points); + + const UpdateFlags flags(fe_data.update_once | fe_data.update_each); + + for (unsigned int k=0; kdofs_per_cell; ++k) + { + if (flags & update_values) + for (unsigned int i=0; i +unsigned int +FE_Q_Hierarchical::n_base_elements () const +{ + return 1; +}; + + + +template +const FiniteElement & +FE_Q_Hierarchical::base_element (const unsigned int index) const +{ + Assert (index==0, ExcIndexRange(index, 0, 1)); + return *this; +}; + + + +template +unsigned int +FE_Q_Hierarchical::element_multiplicity (const unsigned int index) const +{ + Assert (index==0, ExcIndexRange(index, 0, 1)); + return 1; +}; + + + +template +bool +FE_Q_Hierarchical::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::faces_per_cell, + ExcIndexRange (face_index, 0, GeometryInfo::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::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::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::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 +std::vector +FE_Q_Hierarchical::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 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 embedding_dofs (GeometryInfo::vertices_per_cell); + for (unsigned int i=0; i::vertices_per_cell; ++i) + embedding_dofs[i] = i; + + return embedding_dofs; + } + else if (sub_degree==degree) + { + std::vector embedding_dofs (this->dofs_per_cell); + for (unsigned int i=0; idofs_per_cell; ++i) + embedding_dofs[i] = i; + + return embedding_dofs; + } + + if ((dim==2) || (dim==3)) + { + std::vector 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::vertices_per_cell) + embedding_dofs[i] = i; + // line + else if (i<(GeometryInfo::vertices_per_cell + + GeometryInfo::lines_per_cell * (sub_degree-1))) + { + const unsigned int j = (i - GeometryInfo::vertices_per_cell) % + (sub_degree-1); + const unsigned int line = (i - GeometryInfo::vertices_per_cell - j) / (sub_degree-1); + + embedding_dofs[i] = GeometryInfo::vertices_per_cell + + line * (degree-1) + j; + } + // quad + else if (i<(GeometryInfo::vertices_per_cell + + GeometryInfo::lines_per_cell * (sub_degree-1)) + + GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1)) + { + const unsigned int j = (i - GeometryInfo::vertices_per_cell - + GeometryInfo::lines_per_cell * (sub_degree-1)) % (sub_degree-1); + const unsigned int k = ( (i - GeometryInfo::vertices_per_cell - + GeometryInfo::lines_per_cell * (sub_degree-1) - j) / (sub_degree-1) ) % (sub_degree-1); + const unsigned int face = (i - GeometryInfo::vertices_per_cell - + GeometryInfo::lines_per_cell * (sub_degree-1) - k * (sub_degree-1) - j) / ( (sub_degree-1) * (sub_degree-1) ); + + embedding_dofs[i] = GeometryInfo::vertices_per_cell + + GeometryInfo::lines_per_cell * (degree-1) + + face * (degree-1) * (degree-1) + + k * (degree-1) + j; + } + // hex + else if (i<(GeometryInfo::vertices_per_cell + + GeometryInfo::lines_per_cell * (sub_degree-1)) + + GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1) + + GeometryInfo::hexes_per_cell * (sub_degree-1) * (sub_degree-1) * (sub_degree-1)) + { + const unsigned int j = (i - GeometryInfo::vertices_per_cell - + GeometryInfo::lines_per_cell * (sub_degree-1) - + GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1) ) % (sub_degree-1); + const unsigned int k = ( (i - GeometryInfo::vertices_per_cell - + GeometryInfo::lines_per_cell * (sub_degree-1) - + GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1) - j) / (sub_degree-1) ) % (sub_degree-1); + const unsigned int l = (i - GeometryInfo::vertices_per_cell - + GeometryInfo::lines_per_cell * (sub_degree-1) - + GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1) - j - k * (sub_degree-1)) / ( (sub_degree-1) * (sub_degree-1) ); + + embedding_dofs[i] = GeometryInfo::vertices_per_cell + + GeometryInfo::lines_per_cell * (degree-1) + + GeometryInfo::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 (); + } +} + + + +template +unsigned int +FE_Q_Hierarchical::memory_consumption () const +{ + Assert (false, ExcNotImplemented ()); + return 0; +} + + + +template +unsigned int +FE_Q_Hierarchical::get_degree () const +{ + return degree; +}; + + + +template class FE_Q_Hierarchical; -- 2.39.5