From fc7043dc2865e9f94db9e0fb055fe1e46a3a3c6a Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 9 Jun 2003 15:56:35 +0000 Subject: [PATCH] Merge branch_raviart_thomas git-svn-id: https://svn.dealii.org/trunk@7770 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_base.h | 52 +- deal.II/deal.II/include/fe/fe_dgp.h | 11 + .../deal.II/include/fe/fe_dgp_nonparametric.h | 11 + deal.II/deal.II/include/fe/fe_dgq.h | 31 + deal.II/deal.II/include/fe/fe_nedelec.h | 13 +- deal.II/deal.II/include/fe/fe_q.h | 98 +- .../deal.II/include/fe/fe_q_hierarchical.h | 11 + .../deal.II/include/fe/fe_raviart_thomas.h | 636 ++++++ deal.II/deal.II/include/fe/fe_system.h | 56 +- deal.II/deal.II/include/fe/fe_tools.h | 2 +- deal.II/deal.II/source/fe/fe.cc | 19 +- deal.II/deal.II/source/fe/fe_dgp.cc | 25 + .../deal.II/source/fe/fe_dgp_nonparametric.cc | 25 + deal.II/deal.II/source/fe/fe_dgq.cc | 224 ++ deal.II/deal.II/source/fe/fe_nedelec.cc | 36 +- deal.II/deal.II/source/fe/fe_q.cc | 1730 ++++++++------ deal.II/deal.II/source/fe/fe_q_1d.cc | 90 - deal.II/deal.II/source/fe/fe_q_2d.cc | 242 +- deal.II/deal.II/source/fe/fe_q_3d.cc | 373 ---- .../deal.II/source/fe/fe_q_hierarchical.cc | 26 + .../deal.II/source/fe/fe_raviart_thomas.cc | 1983 +++++++++++++++++ deal.II/deal.II/source/fe/fe_system.cc | 124 ++ deal.II/deal.II/source/fe/fe_tools.cc | 73 +- 23 files changed, 4460 insertions(+), 1431 deletions(-) create mode 100644 deal.II/deal.II/include/fe/fe_raviart_thomas.h create mode 100644 deal.II/deal.II/source/fe/fe_raviart_thomas.cc diff --git a/deal.II/deal.II/include/fe/fe_base.h b/deal.II/deal.II/include/fe/fe_base.h index 9a00eee928..85f16c573d 100644 --- a/deal.II/deal.II/include/fe/fe_base.h +++ b/deal.II/deal.II/include/fe/fe_base.h @@ -24,6 +24,8 @@ #include #include +#include + template class FESystem; @@ -440,6 +442,26 @@ class FiniteElementBase : public Subscriptor, const std::vector &restriction_is_additive_flags, const std::vector > &nonzero_components); + /** + * Return a string that uniquely + * identifies a finite + * element. The general + * convention is that this is the + * class name, followed by the + * space dimension in angle + * brackets, and the polynomial + * degree and whatever else is + * necessary in parentheses. For + * example, @p{FE_Q<2>(3)} is the + * value returned for a cubic + * element in 2d. + * + * Systems of elements have their + * own naming convention, see the + * @ref{FESystem} class. + */ + virtual std::string get_name () const = 0; + /** * Return the value of the * @p{i}th shape function at the @@ -774,7 +796,31 @@ class FiniteElementBase : public Subscriptor, * just expresses. */ bool constraints_are_implemented () const; - + + /** + * Return the matrix + * interpolating from the given + * finite element to the present + * one. The size of the matrix is + * then @p{dofs_per_cell} times + * @p{source.dofs_per_cell}. + * + * Derived elements will have to + * implement this function. They + * may only provide interpolation + * matrices for certain source + * finite elements, for example + * those from the same family. If + * they don't implement + * interpolation from a given + * element, then they must throw + * an exception of type + * @ref{FiniteElementBase::ExcInterpolationNotImplemented}. + */ + virtual void + get_interpolation_matrix (const FiniteElementBase &source, + FullMatrix &matrix) const; + /** * Comparison operator. We also * check for equality of the @@ -1277,6 +1323,10 @@ class FiniteElementBase : public Subscriptor, << "The interface matrix has a size of " << arg1 << "x" << arg2 << ", which is not reasonable in the present dimension."); + /** + * Exception + */ + DeclException0 (ExcInterpolationNotImplemented); protected: /** diff --git a/deal.II/deal.II/include/fe/fe_dgp.h b/deal.II/deal.II/include/fe/fe_dgp.h index ad0ff15b0c..09c7be6112 100644 --- a/deal.II/deal.II/include/fe/fe_dgp.h +++ b/deal.II/deal.II/include/fe/fe_dgp.h @@ -43,6 +43,17 @@ class FE_DGP : public FiniteElement */ FE_DGP (const unsigned int k); + /** + * Return a string that uniquely + * identifies a finite + * element. This class returns + * @p{FE_DGP(degree)}, with + * @p{dim} and @p{degree} + * replaced by appropriate + * values. + */ + virtual std::string get_name () const; + /** * Return the value of the * @p{i}th shape function at the diff --git a/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h b/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h index 40d0ed1cdd..72e594c0f2 100644 --- a/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h +++ b/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h @@ -52,6 +52,17 @@ class FE_DGPNonparametric : public FiniteElement */ FE_DGPNonparametric (const unsigned int k); + /** + * Return a string that uniquely + * identifies a finite + * element. This class returns + * @p{FE_DGPNonparametric(degree)}, + * with @p{dim} and @p{degree} + * replaced by appropriate + * values. + */ + virtual std::string get_name () const; + /** * Return the value of the * @p{i}th shape function at the diff --git a/deal.II/deal.II/include/fe/fe_dgq.h b/deal.II/deal.II/include/fe/fe_dgq.h index 907ba8fd60..f5cd910154 100644 --- a/deal.II/deal.II/include/fe/fe_dgq.h +++ b/deal.II/deal.II/include/fe/fe_dgq.h @@ -42,6 +42,17 @@ class FE_DGQ : public FiniteElement */ FE_DGQ (const unsigned int k); + /** + * Return a string that uniquely + * identifies a finite + * element. This class returns + * @p{FE_DGQ(degree)}, with + * @p{dim} and @p{degree} + * replaced by appropriate + * values. + */ + virtual std::string get_name () const; + /** * Return the value of the * @p{i}th shape function at the @@ -152,6 +163,26 @@ class FE_DGQ : public FiniteElement * constructor. */ unsigned int get_degree () const; + + /** + * Return the matrix + * interpolating from the given + * finite element to the present + * one. The size of the matrix is + * then @p{dofs_per_cell} times + * @p{source.dofs_per_cell}. + * + * These matrices are only + * available if the source + * element is also a @p{FE_DGQ} + * element. Otherwise, an + * exception of type + * @ref{FiniteElementBase::ExcInterpolationNotImplemented} + * is thrown. + */ + virtual void + get_interpolation_matrix (const FiniteElementBase &source, + FullMatrix &matrix) const; /** * Number of base elements in a diff --git a/deal.II/deal.II/include/fe/fe_nedelec.h b/deal.II/deal.II/include/fe/fe_nedelec.h index 2f4ea09c50..8cb0977a14 100644 --- a/deal.II/deal.II/include/fe/fe_nedelec.h +++ b/deal.II/deal.II/include/fe/fe_nedelec.h @@ -174,7 +174,7 @@ template class MappingQ; * implemented there. * * - * @author Wolfgang Bangerth, Anna Schneebeli, 2002 + * @author Wolfgang Bangerth, Anna Schneebeli, 2002, 2003 */ template class FE_Nedelec : public FiniteElement @@ -186,6 +186,17 @@ class FE_Nedelec : public FiniteElement */ FE_Nedelec (const unsigned int p); + /** + * Return a string that uniquely + * identifies a finite + * element. This class returns + * @p{FE_Nedelec(degree)}, with + * @p{dim} and @p{degree} + * replaced by appropriate + * values. + */ + virtual std::string get_name () const; + /** * Return the value of the * @p{component}th vector diff --git a/deal.II/deal.II/include/fe/fe_q.h b/deal.II/deal.II/include/fe/fe_q.h index 218949d220..66f3627906 100644 --- a/deal.II/deal.II/include/fe/fe_q.h +++ b/deal.II/deal.II/include/fe/fe_q.h @@ -246,6 +246,17 @@ class FE_Q : public FiniteElement */ FE_Q (const unsigned int p); + /** + * Return a string that uniquely + * identifies a finite + * element. This class returns + * @p{FE_Q(degree)}, with + * @p{dim} and @p{degree} + * replaced by appropriate + * values. + */ + virtual std::string get_name () const; + /** * Return the value of the * @p{i}th shape function at the @@ -358,6 +369,26 @@ class FE_Q : public FiniteElement unsigned int get_degree () const; /** + * Return the matrix + * interpolating from the given + * finite element to the present + * one. The size of the matrix is + * then @p{dofs_per_cell} times + * @p{source.dofs_per_cell}. + * + * These matrices are only + * available if the source + * element is also a @p{FE_Q} + * element. Otherwise, an + * exception of type + * @ref{FiniteElementBase::ExcInterpolationNotImplemented} + * is thrown. + */ + virtual void + get_interpolation_matrix (const FiniteElementBase &source, + FullMatrix &matrix) const; + + /** * Number of base elements in a * mixed discretization. Since * this is a scalar element, @@ -426,42 +457,6 @@ class FE_Q : public FiniteElement */ struct Matrices { - /** - * Embedding matrices. For - * each element type (the - * first index) there are as - * many embedding matrices as - * there are children per - * cell. The first index - * starts with linear - * elements and goes up in - * polynomial degree. The - * array may grow in the - * future with the number of - * elements for which these - * matrices have been - * computed. If for some - * element, the matrices have - * not been computed then you - * may use the element - * nevertheless but can not - * access the respective - * fields. - */ - static const double * const - embedding[][GeometryInfo::children_per_cell]; - - /** - * Number of elements (first - * index) the above field - * has. Equals the highest - * polynomial degree for - * which the embedding - * matrices have been - * computed. - */ - static const unsigned int n_embedding_matrices; - /** * As the * @p{embedding_matrices} @@ -858,41 +853,30 @@ template <> std::vector FE_Q<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int); -// declaration of explicit specializations of member variables, if the -// compiler allows us to do that (the standard says we must) -#ifndef DEAL_II_MEMBER_VAR_SPECIALIZATION_BUG -template <> -const double * const -FE_Q<1>::Matrices::embedding[][GeometryInfo<1>::children_per_cell]; - template <> -const unsigned int FE_Q<1>::Matrices::n_embedding_matrices; +void FE_Q<1>::initialize_constraints (); template <> -const double * const FE_Q<1>::Matrices::constraint_matrices[]; +void FE_Q<2>::initialize_constraints (); template <> -const unsigned int FE_Q<1>::Matrices::n_constraint_matrices; +void FE_Q<3>::initialize_constraints (); -template <> -const double * const -FE_Q<2>::Matrices::embedding[][GeometryInfo<2>::children_per_cell]; +// declaration of explicit specializations of member variables, if the +// compiler allows us to do that (the standard says we must) +#ifndef DEAL_II_MEMBER_VAR_SPECIALIZATION_BUG template <> -const unsigned int FE_Q<2>::Matrices::n_embedding_matrices; +const double * const FE_Q<1>::Matrices::constraint_matrices[]; template <> -const double * const FE_Q<2>::Matrices::constraint_matrices[]; +const unsigned int FE_Q<1>::Matrices::n_constraint_matrices; template <> -const unsigned int FE_Q<2>::Matrices::n_constraint_matrices; - -template <> -const double * const -FE_Q<3>::Matrices::embedding[][GeometryInfo<3>::children_per_cell]; +const double * const FE_Q<2>::Matrices::constraint_matrices[]; template <> -const unsigned int FE_Q<3>::Matrices::n_embedding_matrices; +const unsigned int FE_Q<2>::Matrices::n_constraint_matrices; template <> const double * const FE_Q<3>::Matrices::constraint_matrices[]; diff --git a/deal.II/deal.II/include/fe/fe_q_hierarchical.h b/deal.II/deal.II/include/fe/fe_q_hierarchical.h index 816bb73f68..9823c87351 100644 --- a/deal.II/deal.II/include/fe/fe_q_hierarchical.h +++ b/deal.II/deal.II/include/fe/fe_q_hierarchical.h @@ -241,6 +241,17 @@ class FE_Q_Hierarchical : public FiniteElement */ FE_Q_Hierarchical (const unsigned int p); + /** + * Return a string that uniquely + * identifies a finite + * element. This class returns + * @p{FE_Q_Hierarchical(degree)}, + * with @p{dim} and @p{degree} + * replaced by appropriate + * values. + */ + virtual std::string get_name () const; + /** * Return the value of the * @p{i}th shape function at the diff --git a/deal.II/deal.II/include/fe/fe_raviart_thomas.h b/deal.II/deal.II/include/fe/fe_raviart_thomas.h new file mode 100644 index 0000000000..00011181cb --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_raviart_thomas.h @@ -0,0 +1,636 @@ +//--------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2002, 2003 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_raviart_thomas_h +#define __deal2__fe_raviart_thomas_h + +#include +#include +#include +#include + +template class MappingQ; + + + +/** + * Implementation of continuous Raviart-Thomas elements for the space + * H_div. Note, however, that continuity only concerns the normal + * component of the vector field. + * + * The constructor of this class takes the degree @p{p} of this finite + * element. The numbering of the degree of this element in the + * literature is somewhat funny: the degree is defined not as the + * polynomial degree of the finite element space, but as that of the + * normal component of the traces onto the boundary. Thus, the lowest + * order, zero, has linear shape functions, but on the faces, the + * traces of the normal component of these elements is constant on + * each face. + * + * + * @sect3{Interpolation to finer and coarser meshes} + * + * Each finite element class in deal.II provides matrices that are + * used to interpolate from coarser to finer meshes and the other way + * round. Interpolation from a mother cell to its children is usually + * trivial, since finite element spaces are normally nested and this + * kind of interpolation is therefore exact. On the other hand, when + * we interpolate from child cells to the mother cell, we usually have + * to throw away some information. + * + * For continuous elements, this transfer usually happens by + * interpolating the values on the child cells at the support points + * of the shape functions of the mother cell. However, for + * discontinuous elements, we often use a projection from the child + * cells to the mother cell. The projection approach is only possible + * for discontinuous elements, since it cannot be guaranteed that the + * values of the projected functions on one cell and its neighbor + * match. In this case, only an interpolation can be + * used. (Internally, whether the values of a shape function are + * interpolated or projected, or better: whether the matrices the + * finite element provides are to be treated with the properties of a + * projection or of an interpolation, is controlled by the + * @p{restriction_is_additive} flag. See there for more information.) + * + * Here, things are not so simple: since the element has some + * continuity requirements across faces, we can only resort to some + * kind of interpolation. On the other hand, for the lowest order + * elements, the values of generating functionals are the (constant) + * tangential values of the shape functions. We would therefore really + * like to take the mean value of the tangential values of the child + * faces, and make this the value of the mother face. Then, however, + * taking a mean value of two piecewise constant function is not an + * interpolation, but a restriction. Since this is not possible, we + * cannot use this. + * + * To make a long story somewhat shorter, when interpolating from + * refined edges to a coarse one, we do not take the mean value, but + * pick only one (the one from the first child edge). While this is + * not optimal, it is certainly a valid choice (using an interpolation + * point that is not in the middle of the cell, but shifted to one + * side), and it also preserves the order of the interpolation. + * + * + * @sect3{Numbering of the degrees of freedom (DoFs)} + * + * Nedelec elements have their degrees of freedom on edges, with shape + * functions being vector valued and pointing in tangential + * direction. We use the standard enumeration and direction of edges + * in deal.II, yielding the following shape functions in 2d: + * + * @begin{verbatim} + * 2 + * *---^---* + * | | + * 3> >1 + * | | + * *---^---* + * 0 + * @end{verbatim} + * + * For the 3d case, the ordering follows the same scheme: the lines + * are numbered as described in the documentation of the + * @ref{Triangulation} class, i.e. + * @begin{verbatim} + * *---6---* *---6---* + * /| | / /| + * 11 | 5 11 10 5 + * / 7 | / / | + * * | | *---2---* | + * | *---4---* | | * + * | / / | 1 / + * 3 8 9 3 | 9 + * |/ / | |/ + * *---0---* *---0---* + * @end{verbatim} + * and their directions are as follows: + * @begin{verbatim} + * *--->---* *--->---* + * /| | / /| + * ^ | ^ ^ ^ ^ + * / ^ | / / | + * * | | *--->---* | + * | *--->---* | | * + * | / / | ^ / + * ^ ^ ^ ^ | ^ + * |/ / | |/ + * *--->---* *--->---* + * @end{verbatim} + * + * The element does not make much sense in 1d, so it is not + * implemented there. + * + * + * @author Wolfgang Bangerth, 2003 + */ +template +class FE_RaviartThomas : public FiniteElement +{ + public: + /** + * Constructor for the Nedelec + * element of degree @p{p}. + */ + FE_RaviartThomas (const unsigned int p); + + /** + * Return a string that uniquely + * identifies a finite + * element. This class returns + * @p{FE_RaviartThomas(degree)}, with + * @p{dim} and @p{degree} + * replaced by appropriate + * values. + */ + virtual std::string get_name () 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. + */ + virtual double shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) 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. + */ + virtual Tensor<1,dim> shape_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) 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. + */ + 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; + + /** + * Return the matrix + * interpolating from the given + * finite element to the present + * one. The size of the matrix is + * then @p{dofs_per_cell} times + * @p{source.dofs_per_cell}. + * + * These matrices are only + * available if the source + * element is also a Raviart + * Thomas element. Otherwise, an + * exception of type + * @ref{FiniteElementBase::ExcInterpolationNotImplemented} + * is thrown. + */ + virtual void + get_interpolation_matrix (const FiniteElementBase &source, + FullMatrix &matrix) const; + + /** + * Number of base elements in a + * mixed discretization. Here, + * this is of course equal to + * one. + */ + virtual unsigned int n_base_elements () const; + + /** + * Access to base element + * objects. Since this element is + * atomic, @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 an + * atomic 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; + + /** + * This function returns + * @p{true}, if the shape + * function @p{shape_index} has + * non-zero values on the face + * @p{face_index}. For the lowest + * order Nedelec elements, this + * is actually the case for the + * one on which the shape + * function is defined and all + * neighboring ones. + * + * 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; + + /** + * Exception + */ + DeclException0 (ExcNotUsefulInThisDimension); + + 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: + /** + * Degree of the polynomials. + */ + const unsigned int degree; + + /** + * Spaces describing the + * anisotropic polynomial spaces + * for each vector component, + * i.e. there are @p{dim} + * elements of this field. The + * values for this member are + * created in + * @ref{create_polynomials}. + */ + const std::vector > polynomials; + + /** + * For each shape function, store + * to which vector component (on + * the unit cell, they are mixed + * on the real cell by the + * transformation) they belong, + * and which index they have + * within the anisotropic tensor + * product polynomial space + * describing this vector + * component. + * + * These values are computed by + * the @ref{compute_renumber} + * function. + */ + const std::vector > renumber; + + + /** + * Generate the polynomial spaces + * for the @ref{polynomials} + * member. + */ + static std::vector > + create_polynomials (const unsigned int degree); + + /** + * 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); + + /** + * Compute the vector used for + * the + * @p{restriction_is_additive} + * field passed to the base + * class's constructor. + */ + static std::vector + get_ria_vector (const unsigned int degree); + + /** + * Compute the values of the + * @p{renumber} field. + */ + static std::vector > + compute_renumber (const unsigned int); + + /** + * Initialize the hanging node + * constraints matrices. Called + * from the constructor. + */ + void initialize_constraints (); + + /** + * Initialize the embedding + * matrices. Called from the + * constructor. + */ + void initialize_embedding (); + + /** + * Initialize the restriction + * matrices. Called from the + * constructor. + */ + void initialize_restriction (); + + /** + * Initialize the + * @p{unit_support_points} field + * of the @ref{FiniteElementBase} + * class. Called from the + * constructor. + */ + void initialize_unit_support_points (); + + /** + * Initialize the + * @p{unit_face_support_points} field + * of the @ref{FiniteElementBase} + * class. Called from the + * constructor. + */ + void initialize_unit_face_support_points (); + + /** + * Given a set of flags indicating + * what quantities are requested + * from a @p{FEValues} object, + * return which of these can be + * precomputed once and for + * all. Often, the values of + * shape function at quadrature + * points can be precomputed, for + * example, in which case the + * return value of this function + * would be the logical and of + * the input @p{flags} and + * @p{update_values}. + * + * For the present kind of finite + * element, this is exactly the + * case. + */ + virtual UpdateFlags update_once (const UpdateFlags flags) const; + + /** + * This is the opposite to the + * above function: given a set of + * flags indicating what we want + * to know, return which of these + * need to be computed each time + * we visit a new cell. + * + * If for the computation of one + * quantity something else is + * also required (for example, we + * often need the covariant + * transformation when gradients + * need to be computed), include + * this in the result as well. + */ + virtual UpdateFlags update_each (const UpdateFlags flags) const; + + /** + * 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. Since + * the shape functions are + * vector-valued (with as + * many components as there + * are space dimensions), the + * value is a tensor. + * + * In this array, we store + * the values of the shape + * function in the quadrature + * points on the unit + * cell. The transformation + * to the real space cell is + * then simply done by + * multiplication with the + * Jacobian of the mapping. + */ + Table<2,Tensor<1,dim> > 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<2,dim> > shape_gradients; + }; + + /** + * Allow access from other + * dimensions. + */ + template friend class FE_RaviartThomas; +}; + + +/* -------------- declaration of explicit specializations ------------- */ + +template <> +void FE_RaviartThomas<1>::initialize_unit_face_support_points (); + +template <> +std::vector FE_RaviartThomas<1>::get_dpo_vector (const unsigned int); + +template <> +std::vector > +FE_RaviartThomas<1>::create_polynomials (const unsigned int); + +template <> +std::vector > +FE_RaviartThomas<2>::create_polynomials (const unsigned int); + +template <> +std::vector > +FE_RaviartThomas<3>::create_polynomials (const unsigned int); + +template <> +std::vector > +FE_RaviartThomas<1>::compute_renumber (const unsigned int); + +template <> +std::vector > +FE_RaviartThomas<2>::compute_renumber (const unsigned int); + +template <> +std::vector > +FE_RaviartThomas<3>::compute_renumber (const unsigned int); + +template <> +void +FE_RaviartThomas<1>::initialize_constraints (); + +template <> +void +FE_RaviartThomas<2>::initialize_constraints (); + +template <> +void +FE_RaviartThomas<3>::initialize_constraints (); + +template <> +void +FE_RaviartThomas<1>::initialize_embedding (); + +template <> +void +FE_RaviartThomas<1>::initialize_restriction (); + +template <> +void +FE_RaviartThomas<2>::initialize_restriction (); + +template <> +void +FE_RaviartThomas<3>::initialize_restriction (); + + +#endif diff --git a/deal.II/deal.II/include/fe/fe_system.h b/deal.II/deal.II/include/fe/fe_system.h index 349bf0efa6..6a9a557ac7 100644 --- a/deal.II/deal.II/include/fe/fe_system.h +++ b/deal.II/deal.II/include/fe/fe_system.h @@ -63,7 +63,7 @@ * coupled to @p{u} at the vertices and the line on the larger cell next to this * vertex, there is no interaction with @p{v} and @p{w} of this or the other cell. * - * @author Wolfgang Bangerth, Guido Kanschat, 1999, partial reimplementation Ralf Hartmann 2001. + * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2002, 2003, partial reimplementation Ralf Hartmann 2001. */ template class FESystem : public FiniteElement @@ -123,6 +123,24 @@ class FESystem : public FiniteElement */ virtual ~FESystem (); + /** + * Return a string that uniquely + * identifies a finite + * element. This element returns + * a string that is composed of + * the strings + * @p{name1}...@p{nameN} returned + * by the basis elements. From + * these, we create a sequence + * @p{FESystem[name1^m1-name2^m2-...-nameN^mN]}, + * where @p{mi} are the + * multiplicities of the basis + * elements. If a multiplicity is + * equal to one, then the + * superscript is omitted. + */ + virtual std::string get_name () const; + /** * Return the value of the * @p{i}th shape function at the @@ -273,11 +291,39 @@ class FESystem : public FiniteElement * the computation of these * values to the base elements. */ - virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i, - const Point &p, - const unsigned int component) const; + virtual + Tensor<2,dim> + shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + + /** + * Return the matrix + * interpolating from the given + * finite element to the present + * one. The size of the matrix is + * then @p{dofs_per_cell} times + * @p{source.dofs_per_cell}. + * + * These matrices are available + * if source and destination + * element are both @p{FESystem} + * elements, have the same number + * of base elements with same + * element multiplicity, and if + * these base elements also + * implement their + * @p{get_interpolation_matrix} + * functions. Otherwise, an + * exception of type + * @ref{FiniteElementBase::ExcInterpolationNotImplemented} + * is thrown. + */ + virtual void + get_interpolation_matrix (const FiniteElementBase &source, + FullMatrix &matrix) const; - /** + /** * Number of different base * elements of this object. * diff --git a/deal.II/deal.II/include/fe/fe_tools.h b/deal.II/deal.II/include/fe/fe_tools.h index e2d02f6f08..c10f96cd31 100644 --- a/deal.II/deal.II/include/fe/fe_tools.h +++ b/deal.II/deal.II/include/fe/fe_tools.h @@ -483,7 +483,7 @@ class FETools DeclException4 (ExcMatrixDimensionMismatch, int, int, int, int, << "This is a " << arg1 << "x" << arg2 << " matrix, " - << "but should be a " << arg1 << "x" << arg2 << " matrix."); + << "but should be a " << arg3 << "x" << arg4 << " matrix."); }; diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 2c24809b47..7b9b389887 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -328,11 +328,28 @@ FiniteElementBase::interface_constraints_size () const static_cast(-1)); } + + +template +void +FiniteElementBase:: +get_interpolation_matrix (const FiniteElementBase &, + FullMatrix &) const +{ + // by default, no interpolation + // implemented. so throw exception, + // as documentation says + AssertThrow (false, + typename FiniteElementBase:: + ExcInterpolationNotImplemented()); +} + template -bool FiniteElementBase::operator == (const FiniteElementBase &f) const +bool +FiniteElementBase::operator == (const FiniteElementBase &f) const { return ((static_cast&>(*this) == static_cast&>(f)) && diff --git a/deal.II/deal.II/source/fe/fe_dgp.cc b/deal.II/deal.II/source/fe/fe_dgp.cc index 878e972180..c0eb0d10b7 100644 --- a/deal.II/deal.II/source/fe/fe_dgp.cc +++ b/deal.II/deal.II/source/fe/fe_dgp.cc @@ -20,6 +20,11 @@ #include #include +#ifdef HAVE_STD_STRINGSTREAM +# include +#else +# include +#endif template @@ -73,6 +78,26 @@ FE_DGP::FE_DGP (const unsigned int degree) +template +std::string +FE_DGP::get_name () const +{ +#ifdef HAVE_STD_STRINGSTREAM + std::ostringstream namebuf; +#else + std::ostrstream namebuf; +#endif + + namebuf << "FE_DGP<" << dim << ">(" << degree << ")"; + +#ifndef HAVE_STD_STRINGSTREAM + namebuf << std::ends; +#endif + return namebuf.str(); +} + + + template FiniteElement * FE_DGP::clone() const diff --git a/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc b/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc index 361e71d0f6..c80225dc07 100644 --- a/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc +++ b/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc @@ -20,6 +20,11 @@ #include #include +#ifdef HAVE_STD_STRINGSTREAM +# include +#else +# include +#endif template @@ -76,6 +81,26 @@ FE_DGPNonparametric::FE_DGPNonparametric (const unsigned int degree) +template +std::string +FE_DGPNonparametric::get_name () const +{ +#ifdef HAVE_STD_STRINGSTREAM + std::ostringstream namebuf; +#else + std::ostrstream namebuf; +#endif + + namebuf << "FE_DGPNonparametric<" << dim << ">(" << degree << ")"; + +#ifndef HAVE_STD_STRINGSTREAM + namebuf << std::ends; +#endif + return namebuf.str(); +} + + + template FiniteElement * FE_DGPNonparametric::clone() const diff --git a/deal.II/deal.II/source/fe/fe_dgq.cc b/deal.II/deal.II/source/fe/fe_dgq.cc index 0e4cb2beed..fc5b6581e5 100644 --- a/deal.II/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/deal.II/source/fe/fe_dgq.cc @@ -22,6 +22,123 @@ #include #include +#ifdef HAVE_STD_STRINGSTREAM +# include +#else +# include +#endif + + + +// namespace for some functions that are used in this file. they are +// specific to numbering conventions used for the FE_DGQ element, and +// are thus not very interesting to the outside world +namespace +{ + // auxiliary type to allow for some + // kind of explicit template + // specialization of the following + // functions + template struct int2type {}; + + + // given an integer N, compute its + // integer square root (if it + // exists, otherwise give up) + unsigned int int_sqrt (const unsigned int N) + { + for (unsigned int i=0; i<=N; ++i) + if (i*i == N) + return i; + Assert (false, ExcInternalError()); + return static_cast(-1); + } + + + // given an integer N, compute its + // integer cube root (if it + // exists, otherwise give up) + unsigned int int_cuberoot (const unsigned int N) + { + for (unsigned int i=0; i<=N; ++i) + if (i*i*i == N) + return i; + Assert (false, ExcInternalError()); + return static_cast(-1); + } + + + // given N, generate i=0...N-1 + // equidistant points in the + // interior of the interval [0,1] + Point<1> + generate_unit_point (const unsigned int i, + const unsigned int N, + const int2type<1> ) + { + Assert (i (.5); + else + { + const double h = 1./(N-1); + return Point<1>(i*h); + } + } + + + // given N, generate i=0...N-1 + // equidistant points in the domain + // [0,1]^2 + Point<2> + generate_unit_point (const unsigned int i, + const unsigned int N, + const int2type<2> ) + { + Assert (i (.5, .5); + else + { + Assert (N>=4, ExcInternalError()); + const unsigned int N1d = int_sqrt(N); + const double h = 1./(N1d-1); + + return Point<2> (i%N1d * h, + i/N1d * h); + } + } + + + + + // given N, generate i=0...N-1 + // equidistant points in the domain + // [0,1]^3 + Point<3> + generate_unit_point (const unsigned int i, + const unsigned int N, + const int2type<3> ) + { + Assert (i (.5, .5, .5); + else + { + Assert (N>=8, ExcInternalError()); + + const unsigned int N1d = int_cuberoot(N); + const double h = 1./(N1d-1); + + return Point<3> (i%N1d * h, + (i/N1d)%N1d * h, + i/(N1d*N1d) * h); + } + } +} + + template @@ -209,6 +326,26 @@ FE_DGQ::FE_DGQ (const unsigned int degree) +template +std::string +FE_DGQ::get_name () const +{ +#ifdef HAVE_STD_STRINGSTREAM + std::ostringstream namebuf; +#else + std::ostrstream namebuf; +#endif + + namebuf << "FE_DGQ<" << dim << ">(" << degree << ")"; + +#ifndef HAVE_STD_STRINGSTREAM + namebuf << std::ends; +#endif + return namebuf.str(); +} + + + template FiniteElement * FE_DGQ::clone() const @@ -414,6 +551,93 @@ FE_DGQ::rotate_indices (std::vector &numbers, +template +void +FE_DGQ:: +get_interpolation_matrix (const FiniteElementBase &x_source_fe, + FullMatrix &interpolation_matrix) const +{ + // this is only implemented, if the + // source FE is also a + // DGQ element + AssertThrow ((x_source_fe.get_name().find ("FE_DGQ<") == 0) + || + (dynamic_cast*>(&x_source_fe) != 0), + typename FiniteElementBase:: + ExcInterpolationNotImplemented()); + + // ok, source is a Q element, so + // we will be able to do the work + const FE_DGQ &source_fe + = dynamic_cast&>(x_source_fe); + + Assert (interpolation_matrix.m() == this->dofs_per_cell, + ExcDimensionMismatch (interpolation_matrix.m(), + this->dofs_per_cell)); + Assert (interpolation_matrix.n() == source_fe.dofs_per_cell, + ExcDimensionMismatch (interpolation_matrix.m(), + source_fe.dofs_per_cell)); + + + // compute the interpolation + // matrices in much the same way as + // we do for the embedding matrices + // from mother to child. + FullMatrix cell_interpolation (this->dofs_per_cell, + this->dofs_per_cell); + FullMatrix source_interpolation (this->dofs_per_cell, + source_fe.dofs_per_cell); + FullMatrix tmp (this->dofs_per_cell, + source_fe.dofs_per_cell); + for (unsigned int j=0; jdofs_per_cell; ++j) + { + // generate a point on this + // cell and evaluate the + // shape functions there + const Point p = generate_unit_point (j, this->dofs_per_cell, + int2type()); + for (unsigned int i=0; idofs_per_cell; ++i) + cell_interpolation(j,i) + = polynomial_space.compute_value (i, p); + + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int j=0; jdofs_per_cell; ++i) + { + double sum = 0.; + for (unsigned int j=0; j #include +#ifdef HAVE_STD_STRINGSTREAM +# include +#else +# include +#endif + template FE_Nedelec::FE_Nedelec (const unsigned int degree) @@ -58,6 +64,26 @@ FE_Nedelec::FE_Nedelec (const unsigned int degree) +template +std::string +FE_Nedelec::get_name () const +{ +#ifdef HAVE_STD_STRINGSTREAM + std::ostringstream namebuf; +#else + std::ostrstream namebuf; +#endif + + namebuf << "FE_Nedelec<" << dim << ">(" << degree << ")"; + +#ifndef HAVE_STD_STRINGSTREAM + namebuf << std::ends; +#endif + return namebuf.str(); +} + + + template FiniteElement * FE_Nedelec::clone() const @@ -1249,8 +1275,10 @@ FE_Nedelec::fill_fe_subface_values (const Mapping &m if (flags & update_values) { Assert (fe_data.shape_values.n_cols() == - GeometryInfo::faces_per_cell * n_q_points, - ExcInternalError()); + GeometryInfo::subfaces_per_face * + GeometryInfo::faces_per_cell * + n_q_points, + ExcInternalError()); std::vector > shape_values (n_q_points); @@ -1279,7 +1307,9 @@ FE_Nedelec::fill_fe_subface_values (const Mapping &m if (flags & update_gradients) { Assert (fe_data.shape_gradients.n_cols() == - GeometryInfo::faces_per_cell * n_q_points, + GeometryInfo::faces_per_cell * + GeometryInfo::subfaces_per_face * + n_q_points, ExcInternalError()); std::vector > shape_grads1 (n_q_points); diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 84bbfaef08..2d2979464d 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -22,9 +22,27 @@ #include #include +#ifdef HAVE_STD_STRINGSTREAM +# include +#else +# include +#endif + -namespace +// namespace for some functions that are used in this file. they are +// specific to numbering conventions used for the FE_Q element, and +// are thus not very interesting to the outside world +namespace { + // auxiliary type to allow for some + // kind of explicit template + // specialization of the following + // functions + template struct int2type {}; + + // given a permutation array, + // compute and return the inverse + // permutation #ifdef DEAL_II_ANON_NAMESPACE_BUG static #endif @@ -36,6 +54,299 @@ namespace out[in[i]]=i; return out; } + + + // given an integer N, compute its + // integer square root (if it + // exists, otherwise give up) + unsigned int int_sqrt (const unsigned int N) + { + for (unsigned int i=0; i<=N; ++i) + if (i*i == N) + return i; + Assert (false, ExcInternalError()); + return static_cast(-1); + } + + + // given an integer N, compute its + // integer cube root (if it + // exists, otherwise give up) + unsigned int int_cuberoot (const unsigned int N) + { + for (unsigned int i=0; i<=N; ++i) + if (i*i*i == N) + return i; + Assert (false, ExcInternalError()); + return static_cast(-1); + } + + + // given N, generate i=0...N-1 + // equidistant points in the + // interior of the interval [0,1] + Point<1> + generate_unit_point (const unsigned int i, + const unsigned int N, + const int2type<1> ) + { + Assert (i(i*h); + } + + + // given N, generate i=0...N-1 + // equidistant points in the domain + // [0,1]^2 + Point<2> + generate_unit_point (const unsigned int i, + const unsigned int N, + const int2type<2> ) + { + Assert (i=4, ExcInternalError()); + + const unsigned int N1d = int_sqrt(N); + const double h = 1./(N1d-1); + + return Point<2> (i%N1d * h, + i/N1d * h); + } + + + + // given N, generate i=0...N-1 + // equidistant points in the domain + // [0,1]^3 + Point<3> + generate_unit_point (const unsigned int i, + const unsigned int N, + const int2type<3> ) + { + Assert (i=8, ExcInternalError()); + + const unsigned int N1d = int_cuberoot(N); + const double h = 1./(N1d-1); + + return Point<3> (i%N1d * h, + (i/N1d)%N1d * h, + i/(N1d*N1d) * h); + } + + + + // given N, generate i=0...N-1 + // equidistant points in the + // interior of the interval [0,1] + Point<1> + generate_face_unit_point (const unsigned int i, + const unsigned int N, + const int2type<1> ) + { + Assert (i((1+i)*h); + } + + + // given N, generate i=0...N-1 + // equidistant points in the domain + // [0,1]^2, but excluding the four + // vertices (since we don't have to + // consider shape functions on + // child cells that are located on + // existing vertices) + Point<2> + generate_face_unit_point (const unsigned int i, + const unsigned int N, + const int2type<2> ) + { + Assert (i= N1d-2 ? 1 : 0) + + + (i >= N1d*(N1d-1)-2 ? 1 : 0)); + return Point<2> ((true_i%N1d)*h, + (true_i/N1d)*h); + } + + + + // return whether shape function j, + // as given in the numbering + // specific to the computation of + // the constraint matrix, is active + // on the given subface + bool + constraint_function_is_active_on_child (const unsigned int j, + const unsigned int subface, + const FiniteElementData<2> &fe_data) + { + // note that in our weird + // numbering, the zeroth function + // is the one associated with the + // center node, then come the + // ones on subface 0, then those + // on subface 1. the initial one + // is active on both subfaces, + // all other ones only on one of + // the subfaces + return !(((j>=1) && (j<1+fe_data.dofs_per_line) && (subface == 1)) || + ((j>=1+fe_data.dofs_per_line) && (subface == 0))); + } + + + + bool + constraint_function_is_active_on_child (const unsigned int j, + const unsigned int subface, + const FiniteElementData<3> &fe_data) + { + // in 3d: in our weird numbering, + // the zeroth function is the one + // associated with the center + // node, then come the four edge + // midpoints, then the ones on + // the 12 edges then those on + // subfaces. some are active on + // more than one child + + if (j < 5) + // one one of the five vertices + { + switch (j) + { + case 0: return true; + case 1: return (subface == 0) || (subface == 1); + case 2: return (subface == 1) || (subface == 2); + case 3: return (subface == 2) || (subface == 3); + case 4: return (subface == 3) || (subface == 0); + } + } + else if (j < 5 + 12*fe_data.dofs_per_line) + // one one of the 12 lines + { + const unsigned int line = (j-5)/fe_data.dofs_per_line; + Assert (line<12, ExcInternalError()); + + switch (line) + { + case 0: return (subface == 0) || (subface == 1); + case 1: return (subface == 1) || (subface == 2); + case 2: return (subface == 2) || (subface == 3); + case 3: return (subface == 3) || (subface == 0); + case 4: return (subface == 0); + case 5: return (subface == 1); + case 6: return (subface == 1); + case 7: return (subface == 2); + case 8: return (subface == 3); + case 9: return (subface == 2); + case 10: return (subface == 0); + case 11: return (subface == 2); + } + } + else + // interior + { + const unsigned int quad = (j-5-12*fe_data.dofs_per_line)/fe_data.dofs_per_quad; + Assert (quad<4, ExcInternalError()); + return quad == subface; + } + + Assert (false, ExcInternalError()); + return static_cast(-1); + } + + + // given index j in the weird + // constraint numbering, compute + // its index in the polynomials + // space of a given subface + unsigned int + constraint_get_local_j (const unsigned int j, + const unsigned int subface, + const FiniteElementData<2> &fe_data) + { + // the zeroth shape function is a + // little special, since it has + // index N on subface 0 and index + // 0 on subface 1 + + return (subface == 0 ? + (j == 0 ? 1+fe_data.dofs_per_line : j) : + (j == 0 ? 0 : j-fe_data.dofs_per_line)); + } + + + unsigned int + constraint_get_local_j (const unsigned int /*j*/, + const unsigned int /*subface*/, + const FiniteElementData<3> &/*fe_data*/) + { + Assert (false, ExcNotImplemented()); +// const unsigned int N1d = 2+fe_data.dofs_per_line; + return static_cast(-1); + } + + + + // in the constraint numbering: + // return true if the support point + // of shape function j and + // evaluation point i coincide. to + // make things simpler, also pass + // the subface on which j is + // located + bool + constraint_is_support_point (const unsigned int i, + const unsigned int j, + const unsigned int subface, + const FiniteElementData<2> &fe_data) + { + return ((subface == 0) && (((j==0) && (i==fe_data.dofs_per_line)) + || + ((j!=0) && (i==j-1)))) + || + ((subface == 1) && (((j==0) && (i==fe_data.dofs_per_line)) + || + ((j!=0) && (i==j)))); + } + + + bool + constraint_is_support_point (const unsigned int /*i*/, + const unsigned int /*j*/, + const unsigned int /*subface*/, + const FiniteElementData<3> &/*fe_data*/) + { + Assert (false, ExcNotImplemented()); + return false; + } } @@ -71,6 +382,26 @@ FE_Q::FE_Q (const unsigned int degree) +template +std::string +FE_Q::get_name () const +{ +#ifdef HAVE_STD_STRINGSTREAM + std::ostringstream namebuf; +#else + std::ostrstream namebuf; +#endif + + namebuf << "FE_Q<" << dim << ">(" << degree << ")"; + +#ifndef HAVE_STD_STRINGSTREAM + namebuf << std::ends; +#endif + return namebuf.str(); +} + + + template FiniteElement * FE_Q::clone() const @@ -150,6 +481,92 @@ FE_Q::shape_grad_grad_component (const unsigned int i, } + +template +void +FE_Q:: +get_interpolation_matrix (const FiniteElementBase &x_source_fe, + FullMatrix &interpolation_matrix) const +{ + // this is only implemented, if the + // source FE is also a + // Q element + AssertThrow ((x_source_fe.get_name().find ("FE_Q<") == 0) + || + (dynamic_cast*>(&x_source_fe) != 0), + typename FiniteElementBase:: + ExcInterpolationNotImplemented()); + + // ok, source is a Q element, so + // we will be able to do the work + const FE_Q &source_fe + = dynamic_cast&>(x_source_fe); + + Assert (interpolation_matrix.m() == this->dofs_per_cell, + ExcDimensionMismatch (interpolation_matrix.m(), + this->dofs_per_cell)); + Assert (interpolation_matrix.n() == source_fe.dofs_per_cell, + ExcDimensionMismatch (interpolation_matrix.m(), + source_fe.dofs_per_cell)); + + + // compute the interpolation + // matrices in much the same way as + // we do for the embedding matrices + // from mother to child. + FullMatrix cell_interpolation (this->dofs_per_cell, + this->dofs_per_cell); + FullMatrix source_interpolation (this->dofs_per_cell, + source_fe.dofs_per_cell); + FullMatrix tmp (this->dofs_per_cell, + source_fe.dofs_per_cell); + for (unsigned int j=0; jdofs_per_cell; ++j) + { + // generate a point on this + // cell and evaluate the + // shape functions there + const Point p = generate_unit_point (j, this->dofs_per_cell, + int2type()); + for (unsigned int i=0; idofs_per_cell; ++i) + cell_interpolation(renumber[j],renumber[i]) + = polynomial_space.compute_value (i, p); + + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int j=0; jdofs_per_cell; ++i) + { + double sum = 0.; + for (unsigned int j=0; j::lexicographic_to_hierarchic_numbering (const FiniteElementData & const unsigned int n = degree+1; - if (degree == 0) { Assert ((fe_data.dofs_per_vertex == 0) && @@ -288,34 +704,34 @@ FE_Q::lexicographic_to_hierarchic_numbering (const FiniteElementData & switch (dim) { case 1: - { - const unsigned int values[GeometryInfo<1>::vertices_per_cell] - = { 0, degree }; - index = values[i]; - break; - }; + { + const unsigned int values[GeometryInfo<1>::vertices_per_cell] + = { 0, degree }; + index = values[i]; + break; + }; case 2: - { - const unsigned int values[GeometryInfo<2>::vertices_per_cell] - = { 0, degree, n*degree+degree, n*degree }; - index = values[i]; - break; - }; + { + const unsigned int values[GeometryInfo<2>::vertices_per_cell] + = { 0, degree, n*degree+degree, n*degree }; + index = values[i]; + break; + }; case 3: - { - const unsigned int values[GeometryInfo<3>::vertices_per_cell] - = { 0, degree, - n*n*degree + degree, n*n*degree, - n*degree, n*degree+degree, - n*n*degree + n*degree+degree, n*n*degree + n*degree}; - index = values[i]; - break; - }; + { + const unsigned int values[GeometryInfo<3>::vertices_per_cell] + = { 0, degree, + n*n*degree + degree, n*n*degree, + n*degree, n*degree+degree, + n*n*degree + n*degree+degree, n*n*degree + n*degree}; + index = values[i]; + break; + }; default: - Assert(false, ExcNotImplemented()); + Assert(false, ExcNotImplemented()); } Assert (index::lexicographic_to_hierarchic_numbering (const FiniteElementData & case 100: case 200: case 202: case 300: case 302: case 304: case 306: - incr = 1; - break; - // lines in y-direction + 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 + incr = n; + break; + // lines in z-direction case 301: case 303: case 305: case 307: - incr = n*n; - break; + incr = n*n; + break; default: - Assert(false, ExcNotImplemented()); + Assert(false, ExcNotImplemented()); } switch (i+100*dim) { @@ -368,36 +784,36 @@ FE_Q::lexicographic_to_hierarchic_numbering (const FiniteElementData & case 100: case 200: case 203: case 300: case 303: case 308: - tensorstart = 0; - break; - // x=1 y=z=0 + tensorstart = 0; + break; + // x=1 y=z=0 case 201: case 301: case 309: - tensorstart = degree; - break; - // y=1 x=z=0 + tensorstart = degree; + break; + // y=1 x=z=0 case 202: case 304: case 307: - tensorstart = n*degree; - break; - // x=z=1 y=0 + tensorstart = n*degree; + break; + // x=z=1 y=0 case 310: - tensorstart = n*n*degree+degree; - break; - // z=1 x=y=0 + tensorstart = n*n*degree+degree; + break; + // z=1 x=y=0 case 302: case 311: - tensorstart = n*n*degree; - break; - // x=y=1 z=0 + tensorstart = n*n*degree; + break; + // x=y=1 z=0 case 305: - tensorstart = n*degree+degree; - break; - // y=z=1 x=0 + tensorstart = n*degree+degree; + break; + // y=z=1 x=0 case 306: - tensorstart = n*n*n-n; - break; + tensorstart = n*n*n-n; + break; default: - Assert(false, ExcNotImplemented()); + Assert(false, ExcNotImplemented()); } for (unsigned int jx = 1; jx::lexicographic_to_hierarchic_numbering (const FiniteElementData & switch (i) { case 0: - tensorstart = 0; incx = 1; - if (dim==2) - incy = n; - else - incy = n*n; - break; + tensorstart = 0; incx = 1; + if (dim==2) + incy = n; + else + incy = n*n; + break; case 1: - tensorstart = n*degree; incx = 1; incy = n*n; - break; + tensorstart = n*degree; incx = 1; incy = n*n; + break; case 2: - tensorstart = 0; incx = 1; incy = n; - break; + tensorstart = 0; incx = 1; incy = n; + break; case 3: - tensorstart = degree; incx = n; incy = n*n; - break; + tensorstart = degree; incx = n; incy = n*n; + break; case 4: - tensorstart = n*n*degree; incx = 1; incy = n; - break; + tensorstart = n*n*degree; incx = 1; incy = n; + break; case 5: - tensorstart = 0; incx = n; incy = n*n; - break; + tensorstart = 0; incx = n; incy = n*n; + break; default: - Assert(false, ExcNotImplemented()); + Assert(false, ExcNotImplemented()); } for (unsigned int jy = 1; jy::face_lexicographic_to_hierarchic_numbering (const unsigned int degree } -#if (deal_II_dimension == 1) +#if deal_II_dimension == 1 template <> std::vector @@ -493,50 +909,488 @@ FE_Q<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int) #endif +#if deal_II_dimension == 1 +template <> +void +FE_Q<1>::initialize_constraints () +{ + // no constraints in 1d +} -template +#endif + + +#if deal_II_dimension == 2 + +template <> void -FE_Q::initialize_constraints () -{ - // copy constraint matrices if they - // are defined. otherwise leave them - // at invalid size - if ((dim > 1) && (degree < Matrices::n_constraint_matrices+1)) +FE_Q<2>::initialize_constraints () +{ + const unsigned int dim = 2; + + // restricted to each face, the + // traces of the shape functions is + // an element of P_{k} (in 2d), or + // Q_{k} (in 3d), where k is the + // degree of the element + // + // from this, we interpolate + // between mother and cell + // face. for the general case, this + // may be a little complicated if + // we don't use Lagrange + // interpolation polynomials, since + // then we can't just use point + // interpolation. what we do + // instead is to evaluate at a + // number of points and then invert + // the interpolation matrix. here, + // for the FE_Q elements, we + // actually do have Lagrange + // polynomials, but we still follow + // the general scheme since this + // code here is the master copy for + // what we use in other elements as + // well. however, there are places + // where we make use of the fact + // that we have Lagrange + // interpolation polynomials. + + // mathematically speaking, the + // interpolation process works in + // the following way: on each + // subface, we want that finite + // element solututions from both + // sides coincide. i.e. if a and b + // are expansion coefficients for + // the shape functions from both + // sides, we seek a relation + // between x and y such that + // sum_i a_i phi^c_i(x) + // == sum_j b_j phi_j(x) + // for all points x on the + // interface. here, phi^c_i are the + // shape functions on the small + // cell on one side of the face, + // and phi_j those on the big cell + // on the other side. To get this + // relation, it suffices to look at + // a sufficient number of points + // for which this has to hold. if + // there are n functions, then we + // need n evaluation points, and we + // choose them equidistantly. + // + // what one then gets is a matrix + // system + // a A == b B + // where + // A_ij = phi^c_i(x_j) + // B_ij = phi_i(x_j) + // and the relation we are looking for + // is + // a = (A^T)^-1 B^T b + // + // below, we build up these + // matrices, but rather than + // transposing them after the + // fact, we do so while building + // them. A will be + // subface_interpolation, B will be + // face_interpolation. note that we + // build up these matrices for all + // faces at once, rather than + // considering them separately. the + // reason is that we finally will + // want to have them in this order + // anyway, as this is the format we + // need inside deal.II + TensorProductPolynomials + face_polynomials (Polynomials::LagrangeEquidistant:: + generate_complete_basis (degree)); + Assert (face_polynomials.n() == this->dofs_per_face, ExcInternalError()); + + const unsigned int n_small_functions = this->interface_constraints_size()[0]; + + FullMatrix face_interpolation (n_small_functions, this->dofs_per_face); + FullMatrix subface_interpolation (n_small_functions, n_small_functions); + + const std::vector + face_renumber_inverse (invert_numbering(face_renumber)); + + for (unsigned int i=0; iinterface_constraints. - TableBase<2,double>::reinit (this->interface_constraints_size()); + // generate a quadrature point + // xi. it is actually not so + // important where this point + // lies, as long as we make + // sure that they are not + // equal. however, we will want + // them to be the (equidistant) + // Lagrange points, since then + // the subface_interpolation + // matrix has a most positive + // property: it is a + // permutation of the identity + // matrix. so create an + // equidistant mesh of points + // in the interior of the face + // (in 2d). for 3d, things are + // somewhat more convoluted as + // usual, since the new (child) + // shape functions are not only + // located in the interior of + // the face, but also on the + // edges, with the exception of + // the four vertices of the + // face. the function we call + // takes care of all this + const Point p_face = generate_face_unit_point (i, n_small_functions, + int2type()); + + // evaluate the big face + // shape function at this + // point. note that the + // numbering of our shape + // functions is different + // from that of the + // polynomial, which orders + // them in the order of + // interpolation points. + // + // face_renumber_inverse will + // get us over this little + // conversion + for (unsigned int j=0; jdofs_per_face; ++j) + { + face_interpolation(i,j) + = face_polynomials.compute_value(face_renumber_inverse[j], p_face); + // if the value is small up + // to round-off, then + // simply set it to zero to + // avoid unwanted fill-in + // of the constraint + // matrices (which would + // then increase the number + // of other DoFs a + // constrained DoF would + // couple to) + if (std::fabs(face_interpolation(i,j)) < 1e-14) + face_interpolation(i,j) = 0; + } + + // then evaluate all the + // small shape functions at + // this point. + for (unsigned int j=0; j::child_cell_from_point (p_face); + + // then check whether small + // shape function number j + // is nonzero on this + // face. as usual with our + // numbering of shape + // functions in constraint + // matrices, this is messy, + // so have a function that + // does this for us + // + // if not active, then the + // entry in the matrix will + // remain zero, and we + // simply go on with the + // next entry + if (! constraint_function_is_active_on_child (j, subface, *this)) + continue; + + // otherwise: compute the + // coordinates of this + // evaluation point on + // the small face + const Point p_subface + = GeometryInfo::cell_to_child_coordinates (p_face, subface); + + // then get the index of + // small shape function j + // on this subface. again, + // divert to a function + // that is specialized for + // this + const unsigned int local_j + = constraint_get_local_j (j, subface, *this); + + // so evaluate this shape + // function there. now, + // since we have been + // careful with our choice + // of evaluation points, + // this is not actually + // necessary: the values of + // the small shape + // functions at these + // points should be either + // zero, and we can + // precompute which they + // are. However, we double + // check just to be sure we + // didn't do something + // stupid... + // + // (we could just set the + // evaluated value, but + // we'd end up with a lot + // of almost-zero entries, + // which will then carry + // over to the final + // result. this clutters up + // the constraint matrices, + // which we want to keep as + // small as possible.) + if (constraint_is_support_point (i, j, subface, *this)) + subface_interpolation(i, j) = 1.; + else + subface_interpolation(i, j) = 0.; + Assert (std::fabs (subface_interpolation(i, j) - + face_polynomials.compute_value(local_j, p_subface)) + < 1e-12, + ExcInternalError()); + } + } + + // what we now want to do is to + // compute + // (subface_intp)^-1 face_intp + // which should give us the + // desired hanging node constraints. + // rather than actually doing this, + // we note that we have constructed + // subface_interpolation to be a + // permutation of the unit matrix. + // rather than doing a gauss jordan + // inversion, we note that the + // inverse is actually given by the + // transpose of the matrix. This has + // the additional benefit of being + // more stable and in particular of + // not adding almost-zeros + this->interface_constraints + .TableBase<2,double>::reinit (this->interface_constraints_size()); + subface_interpolation.Tmmult (this->interface_constraints, + face_interpolation); + + // in 3d we still have the + // constraint matrices, so make the + // check + if (dim == 3) + if (degree < Matrices::n_constraint_matrices+1) + { + FullMatrix x; + x.TableBase<2,double>::reinit (this->interface_constraints_size()); + x.fill (Matrices::constraint_matrices[degree-1]); + + for (unsigned int i=0; iinterface_constraints(i,j)) + < + 1e-14, + ExcInternalError()); + } +} + +#endif + +#if deal_II_dimension == 3 + +template <> +void +FE_Q<3>::initialize_constraints () +{ + // the algorithm for 2d is written + // in a way so that it can be + // extended to 3d as well. however, + // the weird numbering convention + // makes this really really hard, + // so we abandoned this project at + // one point. the plan is to change + // the numbering convention for the + // constraint matrices, and then + // the approach for 2d will be + // readily extendable to 3d as + // well, but until this happens we + // rather prefer to go back to the + // precomputed matrices in 3d + if (degree < Matrices::n_constraint_matrices+1) + { + this->interface_constraints + .TableBase<2,double>::reinit (this->interface_constraints_size()); this->interface_constraints.fill (Matrices::constraint_matrices[degree-1]); - }; + } } +#endif template void FE_Q::initialize_embedding () { - // copy over embedding matrices if - // they are defined - if ((degree < Matrices::n_embedding_matrices+1) && - (Matrices::embedding[degree-1][0] != 0)) - for (unsigned int c=0; c::children_per_cell; ++c) - { - this->prolongation[c].reinit (this->dofs_per_cell, - this->dofs_per_cell); - this->prolongation[c].fill (Matrices::embedding[degree-1][c]); - - // and make sure that the row - // sum is 1 - for (unsigned int row=0; rowdofs_per_cell; ++row) - { - double sum = 0; - for (unsigned int col=0; coldofs_per_cell; ++col) - sum += this->prolongation[c](row,col); - Assert (std::fabs(sum-1.) < 1e-14, - ExcInternalError()); - }; - }; + // compute the interpolation + // matrices in much the same way as + // we do for the constraints. it's + // actually simpler here, since we + // don't have this weird + // renumbering stuff going on + FullMatrix cell_interpolation (this->dofs_per_cell, + this->dofs_per_cell); + FullMatrix subcell_interpolation (this->dofs_per_cell, + this->dofs_per_cell); + for (unsigned int child=0; child::children_per_cell; ++child) + this->prolongation[child].reinit (this->dofs_per_cell, + this->dofs_per_cell); + for (unsigned int child=0; child::children_per_cell; ++child) + { + for (unsigned int j=0; jdofs_per_cell; ++j) + { + // generate a point on + // the child cell and + // evaluate the shape + // functions there + const Point p_subcell = generate_unit_point (j, this->dofs_per_cell, + int2type()); + const Point p_cell = + GeometryInfo::child_to_cell_coordinates (p_subcell, child); + + for (unsigned int i=0; idofs_per_cell; ++i) + { + const double + cell_value = polynomial_space.compute_value (i, p_cell), + subcell_value = polynomial_space.compute_value (i, p_subcell); + + // cut off values that + // are too small. note + // that we have here + // Lagrange + // interpolation + // functions, so they + // should be zero at + // almost all points, + // and one at the + // others, at least on + // the subcells. so set + // them to their exact + // values + // + // the actual cut-off + // value is somewhat + // fuzzy, but it works + // for + // 1e-14*degree*dim, + // which is kind of + // reasonable given + // that we compute the + // values of the + // polynomials via an + // degree-step + // recursion and then + // multiply the + // 1d-values. this + // gives us a linear + // growth in + // degree*dim, times a + // small constant. + if (std::fabs(cell_value) < 2e-14*degree*dim) + cell_interpolation(renumber[j], renumber[i]) = 0.; + else + cell_interpolation(renumber[j], renumber[i]) = cell_value; + + if (std::fabs(subcell_value) < 2e-14*degree*dim) + subcell_interpolation(renumber[j], renumber[i]) = 0.; + else + if (std::fabs(subcell_value-1) < 2e-14*degree*dim) + subcell_interpolation(renumber[j], renumber[i]) = 1.; + else + // we have put our + // evaluation + // points onto the + // interpolation + // points, so we + // should either + // get zeros or + // ones! + Assert (false, ExcInternalError()); + } + } + + // then compute the embedding + // matrix for this child and + // this coordinate + // direction. by the same trick + // as with the constraint + // matrices, don't compute the + // inverse of + // subcell_interpolation, but + // use the fact that we have + // put our interpolation points + // onto the interpolation + // points of the Lagrange + // polynomials used here. then, + // the subcell_interpolation + // matrix is just a permutation + // of the identity matrix and + // its inverse is also its + // transpose + subcell_interpolation.Tmmult (this->prolongation[child], + cell_interpolation); + + // cut off very small values + // here + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int j=0; jdofs_per_cell; ++j) + if (std::fabs(this->prolongation[child](i,j)) < 2e-14*degree*dim) + this->prolongation[child](i,j) = 0.; + + // and make sure that the row + // sum is 1. this must be so + // since for this element, the + // shape functions add up to on + for (unsigned int row=0; rowdofs_per_cell; ++row) + { + double sum = 0; + for (unsigned int col=0; coldofs_per_cell; ++col) + sum += this->prolongation[child](row,col); + Assert (std::fabs(sum-1.) < 2e-14*degree*dim, + ExcInternalError()); + }; + } } @@ -545,559 +1399,131 @@ template void FE_Q::initialize_restriction () { - - // then fill restriction - // matrices. they are hardcoded for - // the first few elements. in - // contrast to the other matrices, - // these are not stored in the - // files fe_q_[123]d.cc, since they - // contain only a rather small - // number of zeros, and storing - // them element-wise is more - // expensive than just setting the - // nonzero elements as done here + // for these Lagrange interpolation + // polynomials, construction of the + // restriction matrices is + // relatively simple. the reason is + // that the interpolation points on + // the mother cell are always also + // interpolation points for some + // shape function on one or the + // other child, because we have + // chosen equidistant Lagrange + // interpolation points for the + // polynomials + // + // so the only thing we have to + // find out is: for each shape + // function on the mother cell, + // which is the child cell + // (possibly more than one) on + // which it is located, and which + // is the corresponding shape + // function there. rather than + // doing it for the shape functions + // on the mother cell, we take the + // interpolation points there are + // also search which shape function + // corresponds to it (too lazy to + // do this mapping by hand) + // + // note that the interpolation + // point of a shape function can be + // on the boundary between + // subcells. in that case, + // restriction from children to + // mother may produce two or more + // entries for a dof on the mother + // cell. however, this doesn't + // hurt: since the element is + // continuous, the contribution + // from each child should yield the + // same result, and since the + // element is non-additive we just + // overwrite one value (compute one + // one child) by the same value + // (compute on a later child), so + // we don't have to care about this for (unsigned int c=0; c::children_per_cell; ++c) this->restriction[c].reinit (this->dofs_per_cell, this->dofs_per_cell); - switch (dim) + for (unsigned int i=0; idofs_per_cell; ++i) { - case 1: // 1d - { - switch (degree) - { - case 1: - this->restriction[0](0,0) = 1; - this->restriction[1](1,1) = 1; - break; - case 2: - this->restriction[0](0,0) = 1; - this->restriction[0](2,1) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](2,0) = 1; - break; - case 3: - this->restriction[0](0,0) = 1; - this->restriction[0](2,3) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](3,2) = 1; - break; - case 4: - this->restriction[0](0,0) = 1; - this->restriction[0](2,3) = 1; - this->restriction[0](3,1) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](3,0) = 1; - this->restriction[1](4,3) = 1; - break; - - default: - { - // in case we don't - // have the matrices - // (yet), reset them to - // zero size. this does - // not prevent the use - // of this FE, but will - // prevent the use of - // these matrices - for (unsigned int i=0; - i::children_per_cell; - ++i) - this->restriction[i].reinit(0,0); - }; - } - break; - }; - - case 2: // 2d - { - switch (degree) - { - case 1: - this->restriction[0](0,0) = 1; - this->restriction[1](1,1) = 1; - this->restriction[2](2,2) = 1; - this->restriction[3](3,3) = 1; - break; - case 2: - this->restriction[0](0,0) = 1; - this->restriction[0](4,1) = 1; - this->restriction[0](7,3) = 1; - this->restriction[0](8,2) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](4,0) = 1; - this->restriction[1](5,2) = 1; - this->restriction[1](8,3) = 1; - this->restriction[2](2,2) = 1; - this->restriction[2](5,1) = 1; - this->restriction[2](6,3) = 1; - this->restriction[2](8,0) = 1; - this->restriction[3](3,3) = 1; - this->restriction[3](6,2) = 1; - this->restriction[3](7,0) = 1; - this->restriction[3](8,1) = 1; - break; - case 3: - this->restriction[0](0,0) = 1; - this->restriction[0](4,5) = 1; - this->restriction[0](10,11) = 1; - this->restriction[0](12,15) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](5,4) = 1; - this->restriction[1](6,7) = 1; - this->restriction[1](13,14) = 1; - this->restriction[2](2,2) = 1; - this->restriction[2](7,6) = 1; - this->restriction[2](9,8) = 1; - this->restriction[2](15,12) = 1; - this->restriction[3](3,3) = 1; - this->restriction[3](8,9) = 1; - this->restriction[3](11,10) = 1; - this->restriction[3](14,13) = 1; - break; - case 4: - this->restriction[0](0,0) = 1; - this->restriction[0](4,5) = 1; - this->restriction[0](5,1) = 1; - this->restriction[0](13,14) = 1; - this->restriction[0](14,3) = 1; - this->restriction[0](16,20) = 1; - this->restriction[0](17,8) = 1; - this->restriction[0](19,11) = 1; - this->restriction[0](20,2) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](5,0) = 1; - this->restriction[1](6,5) = 1; - this->restriction[1](7,8) = 1; - this->restriction[1](8,2) = 1; - this->restriction[1](17,14) = 1; - this->restriction[1](18,20) = 1; - this->restriction[1](20,3) = 1; - this->restriction[1](21,11) = 1; - this->restriction[2](2,2) = 1; - this->restriction[2](8,1) = 1; - this->restriction[2](9,8) = 1; - this->restriction[2](11,3) = 1; - this->restriction[2](12,11) = 1; - this->restriction[2](20,0) = 1; - this->restriction[2](21,5) = 1; - this->restriction[2](23,14) = 1; - this->restriction[2](24,20) = 1; - this->restriction[3](3,3) = 1; - this->restriction[3](10,11) = 1; - this->restriction[3](11,2) = 1; - this->restriction[3](14,0) = 1; - this->restriction[3](15,14) = 1; - this->restriction[3](19,5) = 1; - this->restriction[3](20,1) = 1; - this->restriction[3](22,20) = 1; - this->restriction[3](23,8) = 1; - break; - - default: - { - // in case we don't - // have the matrices - // (yet), reset them to - // zero size. this does - // not prevent the use - // of this FE, but will - // prevent the use of - // these matrices - for (unsigned int i=0; - i::children_per_cell; - ++i) - this->restriction[i].reinit(0,0); - }; - } - break; - }; - - case 3: // 3d - { - switch (degree) - { - case 1: - this->restriction[0](0,0) = 1; - this->restriction[1](1,1) = 1; - this->restriction[2](2,2) = 1; - this->restriction[3](3,3) = 1; - this->restriction[4](4,4) = 1; - this->restriction[5](5,5) = 1; - this->restriction[6](6,6) = 1; - this->restriction[7](7,7) = 1; - break; - case 2: - this->restriction[0](0,0) = 1; - this->restriction[0](8,1) = 1; - this->restriction[0](11,3) = 1; - this->restriction[0](16,4) = 1; - this->restriction[0](20,2) = 1; - this->restriction[0](22,5) = 1; - this->restriction[0](25,7) = 1; - this->restriction[0](26,6) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](8,0) = 1; - this->restriction[1](9,2) = 1; - this->restriction[1](17,5) = 1; - this->restriction[1](20,3) = 1; - this->restriction[1](22,4) = 1; - this->restriction[1](23,6) = 1; - this->restriction[1](26,7) = 1; - this->restriction[2](2,2) = 1; - this->restriction[2](9,1) = 1; - this->restriction[2](10,3) = 1; - this->restriction[2](18,6) = 1; - this->restriction[2](20,0) = 1; - this->restriction[2](23,5) = 1; - this->restriction[2](24,7) = 1; - this->restriction[2](26,4) = 1; - this->restriction[3](3,3) = 1; - this->restriction[3](10,2) = 1; - this->restriction[3](11,0) = 1; - this->restriction[3](19,7) = 1; - this->restriction[3](20,1) = 1; - this->restriction[3](24,6) = 1; - this->restriction[3](25,4) = 1; - this->restriction[3](26,5) = 1; - this->restriction[4](4,4) = 1; - this->restriction[4](12,5) = 1; - this->restriction[4](15,7) = 1; - this->restriction[4](16,0) = 1; - this->restriction[4](21,6) = 1; - this->restriction[4](22,1) = 1; - this->restriction[4](25,3) = 1; - this->restriction[4](26,2) = 1; - this->restriction[5](5,5) = 1; - this->restriction[5](12,4) = 1; - this->restriction[5](13,6) = 1; - this->restriction[5](17,1) = 1; - this->restriction[5](21,7) = 1; - this->restriction[5](22,0) = 1; - this->restriction[5](23,2) = 1; - this->restriction[5](26,3) = 1; - this->restriction[6](6,6) = 1; - this->restriction[6](13,5) = 1; - this->restriction[6](14,7) = 1; - this->restriction[6](18,2) = 1; - this->restriction[6](21,4) = 1; - this->restriction[6](23,1) = 1; - this->restriction[6](24,3) = 1; - this->restriction[6](26,0) = 1; - this->restriction[7](7,7) = 1; - this->restriction[7](14,6) = 1; - this->restriction[7](15,4) = 1; - this->restriction[7](19,3) = 1; - this->restriction[7](21,5) = 1; - this->restriction[7](24,2) = 1; - this->restriction[7](25,0) = 1; - this->restriction[7](26,1) = 1; - break; - case 3: - this->restriction[0](0,0) = 1; - this->restriction[0](8,9) = 1; - this->restriction[0](14,15) = 1; - this->restriction[0](24,25) = 1; - this->restriction[0](32,35) = 1; - this->restriction[0](40,43) = 1; - this->restriction[0](52,55) = 1; - this->restriction[0](56,63) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](9,8) = 1; - this->restriction[1](10,11) = 1; - this->restriction[1](26,27) = 1; - this->restriction[1](33,34) = 1; - this->restriction[1](41,42) = 1; - this->restriction[1](44,47) = 1; - this->restriction[1](57,62) = 1; - this->restriction[2](2,2) = 1; - this->restriction[2](11,10) = 1; - this->restriction[2](13,12) = 1; - this->restriction[2](28,29) = 1; - this->restriction[2](35,32) = 1; - this->restriction[2](46,45) = 1; - this->restriction[2](49,50) = 1; - this->restriction[2](61,58) = 1; - this->restriction[3](3,3) = 1; - this->restriction[3](12,13) = 1; - this->restriction[3](15,14) = 1; - this->restriction[3](30,31) = 1; - this->restriction[3](34,33) = 1; - this->restriction[3](48,51) = 1; - this->restriction[3](54,53) = 1; - this->restriction[3](60,59) = 1; - this->restriction[4](4,4) = 1; - this->restriction[4](16,17) = 1; - this->restriction[4](22,23) = 1; - this->restriction[4](25,24) = 1; - this->restriction[4](36,39) = 1; - this->restriction[4](42,41) = 1; - this->restriction[4](53,54) = 1; - this->restriction[4](58,61) = 1; - this->restriction[5](5,5) = 1; - this->restriction[5](17,16) = 1; - this->restriction[5](18,19) = 1; - this->restriction[5](27,26) = 1; - this->restriction[5](37,38) = 1; - this->restriction[5](43,40) = 1; - this->restriction[5](45,46) = 1; - this->restriction[5](59,60) = 1; - this->restriction[6](6,6) = 1; - this->restriction[6](19,18) = 1; - this->restriction[6](21,20) = 1; - this->restriction[6](29,28) = 1; - this->restriction[6](39,36) = 1; - this->restriction[6](47,44) = 1; - this->restriction[6](51,48) = 1; - this->restriction[6](63,56) = 1; - this->restriction[7](7,7) = 1; - this->restriction[7](20,21) = 1; - this->restriction[7](23,22) = 1; - this->restriction[7](31,30) = 1; - this->restriction[7](38,37) = 1; - this->restriction[7](50,49) = 1; - this->restriction[7](55,52) = 1; - this->restriction[7](62,57) = 1; - break; - case 4: - this->restriction[0](0,0) = 1; - this->restriction[0](8,9) = 1; - this->restriction[0](9,1) = 1; - this->restriction[0](17,18) = 1; - this->restriction[0](18,3) = 1; - this->restriction[0](32,33) = 1; - this->restriction[0](33,4) = 1; - this->restriction[0](44,48) = 1; - this->restriction[0](45,12) = 1; - this->restriction[0](47,15) = 1; - this->restriction[0](48,2) = 1; - this->restriction[0](62,66) = 1; - this->restriction[0](63,36) = 1; - this->restriction[0](65,21) = 1; - this->restriction[0](66,5) = 1; - this->restriction[0](89,93) = 1; - this->restriction[0](90,30) = 1; - this->restriction[0](92,42) = 1; - this->restriction[0](93,7) = 1; - this->restriction[0](98,111) = 1; - this->restriction[0](99,75) = 1; - this->restriction[0](101,57) = 1; - this->restriction[0](102,24) = 1; - this->restriction[0](107,84) = 1; - this->restriction[0](108,39) = 1; - this->restriction[0](110,27) = 1; - this->restriction[0](111,6) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](9,0) = 1; - this->restriction[1](10,9) = 1; - this->restriction[1](11,12) = 1; - this->restriction[1](12,2) = 1; - this->restriction[1](35,36) = 1; - this->restriction[1](36,5) = 1; - this->restriction[1](45,18) = 1; - this->restriction[1](46,48) = 1; - this->restriction[1](48,3) = 1; - this->restriction[1](49,15) = 1; - this->restriction[1](63,33) = 1; - this->restriction[1](64,66) = 1; - this->restriction[1](66,4) = 1; - this->restriction[1](67,21) = 1; - this->restriction[1](71,75) = 1; - this->restriction[1](72,24) = 1; - this->restriction[1](74,39) = 1; - this->restriction[1](75,6) = 1; - this->restriction[1](99,93) = 1; - this->restriction[1](100,111) = 1; - this->restriction[1](102,30) = 1; - this->restriction[1](103,57) = 1; - this->restriction[1](108,42) = 1; - this->restriction[1](109,84) = 1; - this->restriction[1](111,7) = 1; - this->restriction[1](112,27) = 1; - this->restriction[2](2,2) = 1; - this->restriction[2](12,1) = 1; - this->restriction[2](13,12) = 1; - this->restriction[2](15,3) = 1; - this->restriction[2](16,15) = 1; - this->restriction[2](38,39) = 1; - this->restriction[2](39,6) = 1; - this->restriction[2](48,0) = 1; - this->restriction[2](49,9) = 1; - this->restriction[2](51,18) = 1; - this->restriction[2](52,48) = 1; - this->restriction[2](74,36) = 1; - this->restriction[2](75,5) = 1; - this->restriction[2](77,75) = 1; - this->restriction[2](78,24) = 1; - this->restriction[2](81,42) = 1; - this->restriction[2](82,84) = 1; - this->restriction[2](84,7) = 1; - this->restriction[2](85,27) = 1; - this->restriction[2](108,33) = 1; - this->restriction[2](109,66) = 1; - this->restriction[2](111,4) = 1; - this->restriction[2](112,21) = 1; - this->restriction[2](117,93) = 1; - this->restriction[2](118,111) = 1; - this->restriction[2](120,30) = 1; - this->restriction[2](121,57) = 1; - this->restriction[3](3,3) = 1; - this->restriction[3](14,15) = 1; - this->restriction[3](15,2) = 1; - this->restriction[3](18,0) = 1; - this->restriction[3](19,18) = 1; - this->restriction[3](41,42) = 1; - this->restriction[3](42,7) = 1; - this->restriction[3](47,9) = 1; - this->restriction[3](48,1) = 1; - this->restriction[3](50,48) = 1; - this->restriction[3](51,12) = 1; - this->restriction[3](80,84) = 1; - this->restriction[3](81,39) = 1; - this->restriction[3](83,27) = 1; - this->restriction[3](84,6) = 1; - this->restriction[3](92,33) = 1; - this->restriction[3](93,4) = 1; - this->restriction[3](95,93) = 1; - this->restriction[3](96,30) = 1; - this->restriction[3](107,66) = 1; - this->restriction[3](108,36) = 1; - this->restriction[3](110,21) = 1; - this->restriction[3](111,5) = 1; - this->restriction[3](116,111) = 1; - this->restriction[3](117,75) = 1; - this->restriction[3](119,57) = 1; - this->restriction[3](120,24) = 1; - this->restriction[4](4,4) = 1; - this->restriction[4](20,21) = 1; - this->restriction[4](21,5) = 1; - this->restriction[4](29,30) = 1; - this->restriction[4](30,7) = 1; - this->restriction[4](33,0) = 1; - this->restriction[4](34,33) = 1; - this->restriction[4](53,57) = 1; - this->restriction[4](54,24) = 1; - this->restriction[4](56,27) = 1; - this->restriction[4](57,6) = 1; - this->restriction[4](65,9) = 1; - this->restriction[4](66,1) = 1; - this->restriction[4](68,66) = 1; - this->restriction[4](69,36) = 1; - this->restriction[4](90,18) = 1; - this->restriction[4](91,93) = 1; - this->restriction[4](93,3) = 1; - this->restriction[4](94,42) = 1; - this->restriction[4](101,48) = 1; - this->restriction[4](102,12) = 1; - this->restriction[4](104,111) = 1; - this->restriction[4](105,75) = 1; - this->restriction[4](110,15) = 1; - this->restriction[4](111,2) = 1; - this->restriction[4](113,84) = 1; - this->restriction[4](114,39) = 1; - this->restriction[5](5,5) = 1; - this->restriction[5](21,4) = 1; - this->restriction[5](22,21) = 1; - this->restriction[5](23,24) = 1; - this->restriction[5](24,6) = 1; - this->restriction[5](36,1) = 1; - this->restriction[5](37,36) = 1; - this->restriction[5](54,30) = 1; - this->restriction[5](55,57) = 1; - this->restriction[5](57,7) = 1; - this->restriction[5](58,27) = 1; - this->restriction[5](66,0) = 1; - this->restriction[5](67,9) = 1; - this->restriction[5](69,33) = 1; - this->restriction[5](70,66) = 1; - this->restriction[5](72,12) = 1; - this->restriction[5](73,75) = 1; - this->restriction[5](75,2) = 1; - this->restriction[5](76,39) = 1; - this->restriction[5](102,18) = 1; - this->restriction[5](103,48) = 1; - this->restriction[5](105,93) = 1; - this->restriction[5](106,111) = 1; - this->restriction[5](111,3) = 1; - this->restriction[5](112,15) = 1; - this->restriction[5](114,42) = 1; - this->restriction[5](115,84) = 1; - this->restriction[6](6,6) = 1; - this->restriction[6](24,5) = 1; - this->restriction[6](25,24) = 1; - this->restriction[6](27,7) = 1; - this->restriction[6](28,27) = 1; - this->restriction[6](39,2) = 1; - this->restriction[6](40,39) = 1; - this->restriction[6](57,4) = 1; - this->restriction[6](58,21) = 1; - this->restriction[6](60,30) = 1; - this->restriction[6](61,57) = 1; - this->restriction[6](75,1) = 1; - this->restriction[6](76,36) = 1; - this->restriction[6](78,12) = 1; - this->restriction[6](79,75) = 1; - this->restriction[6](84,3) = 1; - this->restriction[6](85,15) = 1; - this->restriction[6](87,42) = 1; - this->restriction[6](88,84) = 1; - this->restriction[6](111,0) = 1; - this->restriction[6](112,9) = 1; - this->restriction[6](114,33) = 1; - this->restriction[6](115,66) = 1; - this->restriction[6](120,18) = 1; - this->restriction[6](121,48) = 1; - this->restriction[6](123,93) = 1; - this->restriction[6](124,111) = 1; - this->restriction[7](7,7) = 1; - this->restriction[7](26,27) = 1; - this->restriction[7](27,6) = 1; - this->restriction[7](30,4) = 1; - this->restriction[7](31,30) = 1; - this->restriction[7](42,3) = 1; - this->restriction[7](43,42) = 1; - this->restriction[7](56,21) = 1; - this->restriction[7](57,5) = 1; - this->restriction[7](59,57) = 1; - this->restriction[7](60,24) = 1; - this->restriction[7](83,15) = 1; - this->restriction[7](84,2) = 1; - this->restriction[7](86,84) = 1; - this->restriction[7](87,39) = 1; - this->restriction[7](93,0) = 1; - this->restriction[7](94,33) = 1; - this->restriction[7](96,18) = 1; - this->restriction[7](97,93) = 1; - this->restriction[7](110,9) = 1; - this->restriction[7](111,1) = 1; - this->restriction[7](113,66) = 1; - this->restriction[7](114,36) = 1; - this->restriction[7](119,48) = 1; - this->restriction[7](120,12) = 1; - this->restriction[7](122,111) = 1; - this->restriction[7](123,75) = 1; - break; - default: + const Point p_cell = generate_unit_point (i, this->dofs_per_cell, + int2type()); + unsigned int mother_dof = 0; + for (; mother_dofdofs_per_cell; ++mother_dof) + { + const double val + = polynomial_space.compute_value(renumber_inverse[mother_dof], + p_cell); + if (std::fabs (val-1.) < 2e-14*degree*dim) + // ok, this is the right + // dof + break; + else + // make sure that all + // other shape functions + // are zero there + Assert (std::fabs(val) < 2e-14*degree*dim, + ExcInternalError()); + } + // check also the shape + // functions after tat + for (unsigned int j=mother_dof+1; jdofs_per_cell; ++j) + Assert (std::fabs (polynomial_space.compute_value(renumber_inverse[j], + p_cell)) + < 2e-14*degree*dim, + ExcInternalError()); + + // then find the children on + // which the interpolation + // point is located + for (unsigned int child=0; child::children_per_cell; + ++child) + { + // first initialize this + // column of the matrix + for (unsigned int j=0; jdofs_per_cell; ++j) + this->restriction[child](mother_dof, j) = 0.; + + // then check whether this + // interpolation point is + // inside this child cell + const Point p_subcell + = GeometryInfo::cell_to_child_coordinates (p_cell, child); + if (GeometryInfo::is_inside_unit_cell (p_subcell)) { - // in case we don't - // have the matrices - // (yet), reset them to - // zero size. this does - // not prevent the use - // of this FE, but will - // prevent the use of - // these matrices - for (unsigned int i=0; - i::children_per_cell; - ++i) - this->restriction[i].reinit(0,0); - }; - } - break; - }; - - default: - Assert (false, ExcNotImplemented()); + // find the one child + // shape function + // corresponding to + // this point. do it in + // the same way as + // above + unsigned int child_dof = 0; + for (; child_dofdofs_per_cell; ++child_dof) + { + const double val + = polynomial_space.compute_value(renumber_inverse[child_dof], + p_subcell); + if (std::fabs (val-1.) < 2e-14*degree*dim) + break; + else + Assert (std::fabs(val) < 2e-14*degree*dim, + ExcInternalError()); + } + for (unsigned int j=child_dof+1; jdofs_per_cell; ++j) + Assert (std::fabs (polynomial_space.compute_value(renumber_inverse[j], + p_subcell)) + < 2e-14*degree*dim, + ExcInternalError()); + + // so now that we have + // it, set the + // corresponding value + // in the matrix + this->restriction[child](mother_dof, child_dof) = 1.; + } + } } } diff --git a/deal.II/deal.II/source/fe/fe_q_1d.cc b/deal.II/deal.II/source/fe/fe_q_1d.cc index 247f791255..fac6527640 100644 --- a/deal.II/deal.II/source/fe/fe_q_1d.cc +++ b/deal.II/deal.II/source/fe/fe_q_1d.cc @@ -19,96 +19,6 @@ #include -// Transfer matrices for finite elements - -namespace FE_Q_1d -{ - static const double q1_into_q1_refined_0[] = - { - 1., 0., - 13.5/27., 13.5/27., - }; - - static const double q1_into_q1_refined_1[] = - { - 13.5/27., 13.5/27., - 0., 1., - }; - - static const double q2_into_q2_refined_0[] = - { - 1., 0., 0., - 0., 0., 1., - 10.125/27., -3.375/27., 20.25/27., - }; - - static const double q2_into_q2_refined_1[] = - { - 0., 0., 1., - 0., 1., 0., - -3.375/27., 10.125/27., 20.25/27., - }; - - static const double q3_into_q3_refined_0[] = - { - 1., 0., 0., 0., - -1.6875/27., -1.6875/27., 15.1875/27., 15.1875/27., - 8.4375/27., 1.6875/27., 25.3125/27., -8.4375/27., - 0., 0., 1., 0. - }; - - static const double q3_into_q3_refined_1[] = - { - -1.6875/27., -1.6875/27., 15.1875/27., 15.1875/27., - 0., 1., 0., 0., - 0., 0., 0., 1., - 1.6875/27., 8.4375/27., -8.4375/27., 25.3125/27., - }; - - static const double q4_into_q4_refined_0[] = - { - 1., 0., 0., 0., 0., - 0., 0., 0., 1., 0., - 7.3828125/27., -1.0546875/27., 29.53125/27., -14.765625/27., 5.90625/27., - 0., 0., 1., 0., 0., - -1.0546875/27., 0.6328125/27., 12.65625/27., 18.984375/27., -4.21875/27., - }; - - static const double q4_into_q4_refined_1[] = - { - 0., 0., 0., 1., 0., - 0., 1., 0., 0., 0., - 0.6328125/27., -1.0546875/27., -4.21875/27., 18.984375/27., 12.65625/27., - 0., 0., 0., 0., 1., - -1.0546875/27., 7.3828125/27., 5.90625/27., -14.765625/27., 29.53125/27., - }; - -} // namespace FE_Q_1d - - - -// embedding matrices - - -template <> -const double * const -FE_Q<1>::Matrices::embedding[][GeometryInfo<1>::children_per_cell] = -{ - {FE_Q_1d::q1_into_q1_refined_0, FE_Q_1d::q1_into_q1_refined_1}, - {FE_Q_1d::q2_into_q2_refined_0, FE_Q_1d::q2_into_q2_refined_1}, - {FE_Q_1d::q3_into_q3_refined_0, FE_Q_1d::q3_into_q3_refined_1}, - {FE_Q_1d::q4_into_q4_refined_0, FE_Q_1d::q4_into_q4_refined_1}, -}; - - - -template <> -const unsigned int -FE_Q<1>::Matrices::n_embedding_matrices - = sizeof(FE_Q<1>::Matrices::embedding) / - sizeof(FE_Q<1>::Matrices::embedding[0]); - - // No constraints in 1d template <> diff --git a/deal.II/deal.II/source/fe/fe_q_2d.cc b/deal.II/deal.II/source/fe/fe_q_2d.cc index ea7764cc2d..83cdec5544 100644 --- a/deal.II/deal.II/source/fe/fe_q_2d.cc +++ b/deal.II/deal.II/source/fe/fe_q_2d.cc @@ -18,251 +18,17 @@ #include -// Transfer matrices for finite elements -namespace FE_Q_2d -{ - static const double q1_into_q1_refined_0[] = - { - 1., 0., 0., 0., - 13.5/27., 13.5/27., 0., 0., - 6.75/27., 6.75/27., 6.75/27., 6.75/27., - 13.5/27., 0., 0., 13.5/27., - }; - - static const double q1_into_q1_refined_1[] = - { - 13.5/27., 13.5/27., 0., 0., - 0., 1., 0., 0., - 0., 13.5/27., 13.5/27., 0., - 6.75/27., 6.75/27., 6.75/27., 6.75/27., - }; - - static const double q1_into_q1_refined_2[] = - { - 6.75/27., 6.75/27., 6.75/27., 6.75/27., - 0., 13.5/27., 13.5/27., 0., - 0., 0., 1., 0., - 0., 0., 13.5/27., 13.5/27., - }; - - static const double q1_into_q1_refined_3[] = - { - 13.5/27., 0., 0., 13.5/27., - 6.75/27., 6.75/27., 6.75/27., 6.75/27., - 0., 0., 13.5/27., 13.5/27., - 0., 0., 0., 1., - }; - - static const double q2_into_q2_refined_0[] = - { - 1., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 1., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 1., - 0., 0., 0., 0., 0., 0., 0., 1., 0., - 10.125/27., -3.375/27., 0., 0., 20.25/27., 0., 0., 0., 0., - 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 0., 20.25/27., - 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 20.25/27., - 10.125/27., 0., 0., -3.375/27., 0., 0., 0., 20.25/27., 0., - 3.796875/27., -1.265625/27., 0.421875/27., -1.265625/27., 7.59375/27., -2.53125/27., -2.53125/27., 7.59375/27., 15.1875/27., - }; - - static const double q2_into_q2_refined_1[] = - { - 0., 0., 0., 0., 1., 0., 0., 0., 0., - 0., 1., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 1., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 1., - -3.375/27., 10.125/27., 0., 0., 20.25/27., 0., 0., 0., 0., - 0., 10.125/27., -3.375/27., 0., 0., 20.25/27., 0., 0., 0., - 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 20.25/27., - 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 0., 20.25/27., - -1.265625/27., 3.796875/27., -1.265625/27., 0.421875/27., 7.59375/27., 7.59375/27., -2.53125/27., -2.53125/27., 15.1875/27., - }; - - static const double q2_into_q2_refined_2[] = - { - 0., 0., 0., 0., 0., 0., 0., 0., 1., - 0., 0., 0., 0., 0., 1., 0., 0., 0., - 0., 0., 1., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 1., 0., 0., - 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 20.25/27., - 0., -3.375/27., 10.125/27., 0., 0., 20.25/27., 0., 0., 0., - 0., 0., 10.125/27., -3.375/27., 0., 0., 20.25/27., 0., 0., - 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 0., 20.25/27., - 0.421875/27., -1.265625/27., 3.796875/27., -1.265625/27., -2.53125/27., 7.59375/27., 7.59375/27., -2.53125/27., 15.1875/27., - }; - - static const double q2_into_q2_refined_3[] = - { - 0., 0., 0., 0., 0., 0., 0., 1., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 1., - 0., 0., 0., 0., 0., 0., 1., 0., 0., - 0., 0., 0., 1., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 20.25/27., - 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 0., 20.25/27., - 0., 0., -3.375/27., 10.125/27., 0., 0., 20.25/27., 0., 0., - -3.375/27., 0., 0., 10.125/27., 0., 0., 0., 20.25/27., 0., - -1.265625/27., 0.421875/27., -1.265625/27., 3.796875/27., -2.53125/27., -2.53125/27., 7.59375/27., 7.59375/27., 15.1875/27., - }; - - static const double q3_into_q3_refined_0[] = - { - 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - -1.6875/27., -1.6875/27., 0., 0., 15.1875/27., 15.1875/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0.10546875/27., 0.10546875/27., 0.10546875/27., 0.10546875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 8.54296875/27., 8.54296875/27., 8.54296875/27., 8.54296875/27., - -1.6875/27., 0., 0., -1.6875/27., 0., 0., 0., 0., 0., 0., 15.1875/27., 15.1875/27., 0., 0., 0., 0., - 8.4375/27., 1.6875/27., 0., 0., 25.3125/27., -8.4375/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - -0.52734375/27., -0.52734375/27., -0.10546875/27., -0.10546875/27., 4.74609375/27., 4.74609375/27., -1.58203125/27., 0.52734375/27., 0.94921875/27., 0.94921875/27., -1.58203125/27., 0.52734375/27., 14.23828125/27., 14.23828125/27., -4.74609375/27., -4.74609375/27., - 0., 0., 0., 0., 0., 0., -1.6875/27., 0., 0., 0., -1.6875/27., 0., 15.1875/27., 15.1875/27., 0., 0., - -0.52734375/27., -0.10546875/27., -0.10546875/27., -0.52734375/27., -1.58203125/27., 0.52734375/27., 0.94921875/27., 0.94921875/27., -1.58203125/27., 0.52734375/27., 4.74609375/27., 4.74609375/27., 14.23828125/27., -4.74609375/27., 14.23828125/27., -4.74609375/27., - 0., 0., 0., 0., -1.6875/27., 0., 0., 0., -1.6875/27., 0., 0., 0., 15.1875/27., 0., 15.1875/27., 0., - 8.4375/27., 0., 0., 1.6875/27., 0., 0., 0., 0., 0., 0., 25.3125/27., -8.4375/27., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., - 2.63671875/27., 0.52734375/27., 0.10546875/27., 0.52734375/27., 7.91015625/27., -2.63671875/27., 1.58203125/27., -0.52734375/27., 1.58203125/27., -0.52734375/27., 7.91015625/27., -2.63671875/27., 23.73046875/27., -7.91015625/27., -7.91015625/27., 2.63671875/27., - 0., 0., 0., 0., 8.4375/27., 0., 0., 0., 1.6875/27., 0., 0., 0., 25.3125/27., 0., -8.4375/27., 0., - 0., 0., 0., 0., 0., 0., 1.6875/27., 0., 0., 0., 8.4375/27., 0., 25.3125/27., -8.4375/27., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., - }; - - static const double q3_into_q3_refined_1[] = - { - -1.6875/27., -1.6875/27., 0., 0., 15.1875/27., 15.1875/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., -1.6875/27., -1.6875/27., 0., 0., 0., 15.1875/27., 15.1875/27., 0., 0., 0., 0., 0., 0., 0., 0., - 0.10546875/27., 0.10546875/27., 0.10546875/27., 0.10546875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 8.54296875/27., 8.54296875/27., 8.54296875/27., 8.54296875/27., - 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 1.6875/27., 8.4375/27., 0., 0., -8.4375/27., 25.3125/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 8.4375/27., 1.6875/27., 0., 0., 0., 25.3125/27., -8.4375/27., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., -1.6875/27., 0., 0., 0., -1.6875/27., 0., 0., 0., 15.1875/27., 0., 15.1875/27., - -0.10546875/27., -0.52734375/27., -0.52734375/27., -0.10546875/27., 0.52734375/27., -1.58203125/27., 4.74609375/27., 4.74609375/27., 0.52734375/27., -1.58203125/27., 0.94921875/27., 0.94921875/27., -4.74609375/27., 14.23828125/27., -4.74609375/27., 14.23828125/27., - -0.52734375/27., -0.52734375/27., -0.10546875/27., -0.10546875/27., 4.74609375/27., 4.74609375/27., -1.58203125/27., 0.52734375/27., 0.94921875/27., 0.94921875/27., -1.58203125/27., 0.52734375/27., 14.23828125/27., 14.23828125/27., -4.74609375/27., -4.74609375/27., - 0., 0., 0., 0., 0., 0., -1.6875/27., 0., 0., 0., -1.6875/27., 0., 15.1875/27., 15.1875/27., 0., 0., - 0., 0., 0., 0., 0., 8.4375/27., 0., 0., 0., 1.6875/27., 0., 0., 0., 25.3125/27., 0., -8.4375/27., - 0.52734375/27., 2.63671875/27., 0.52734375/27., 0.10546875/27., -2.63671875/27., 7.91015625/27., 7.91015625/27., -2.63671875/27., -0.52734375/27., 1.58203125/27., 1.58203125/27., -0.52734375/27., -7.91015625/27., 23.73046875/27., 2.63671875/27., -7.91015625/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., - 0., 0., 0., 0., 0., 0., 8.4375/27., 0., 0., 0., 1.6875/27., 0., -8.4375/27., 25.3125/27., 0., 0., - }; - - static const double q3_into_q3_refined_2[] = - { - 0.10546875/27., 0.10546875/27., 0.10546875/27., 0.10546875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 8.54296875/27., 8.54296875/27., 8.54296875/27., 8.54296875/27., - 0., -1.6875/27., -1.6875/27., 0., 0., 0., 15.1875/27., 15.1875/27., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., -1.6875/27., -1.6875/27., 0., 0., 0., 0., 15.1875/27., 15.1875/27., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., -1.6875/27., 0., 0., 0., -1.6875/27., 0., 0., 0., 15.1875/27., 0., 15.1875/27., - -0.10546875/27., -0.52734375/27., -0.52734375/27., -0.10546875/27., 0.52734375/27., -1.58203125/27., 4.74609375/27., 4.74609375/27., 0.52734375/27., -1.58203125/27., 0.94921875/27., 0.94921875/27., -4.74609375/27., 14.23828125/27., -4.74609375/27., 14.23828125/27., - 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 1.6875/27., 8.4375/27., 0., 0., 0., -8.4375/27., 25.3125/27., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., - 0., 0., 8.4375/27., 1.6875/27., 0., 0., 0., 0., -8.4375/27., 25.3125/27., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., -1.6875/27., 0., 0., 0., -1.6875/27., 0., 0., 15.1875/27., 15.1875/27., - -0.10546875/27., -0.10546875/27., -0.52734375/27., -0.52734375/27., 0.94921875/27., 0.94921875/27., 0.52734375/27., -1.58203125/27., 4.74609375/27., 4.74609375/27., 0.52734375/27., -1.58203125/27., -4.74609375/27., -4.74609375/27., 14.23828125/27., 14.23828125/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., - 0., 0., 0., 0., 0., 0., 0., 8.4375/27., 0., 0., 0., 1.6875/27., 0., 0., -8.4375/27., 25.3125/27., - 0., 0., 0., 0., 0., 1.6875/27., 0., 0., 0., 8.4375/27., 0., 0., 0., -8.4375/27., 0., 25.3125/27., - 0.10546875/27., 0.52734375/27., 2.63671875/27., 0.52734375/27., -0.52734375/27., 1.58203125/27., -2.63671875/27., 7.91015625/27., -2.63671875/27., 7.91015625/27., -0.52734375/27., 1.58203125/27., 2.63671875/27., -7.91015625/27., -7.91015625/27., 23.73046875/27., - }; - - static const double q3_into_q3_refined_3[] = - { - -1.6875/27., 0., 0., -1.6875/27., 0., 0., 0., 0., 0., 0., 15.1875/27., 15.1875/27., 0., 0., 0., 0., - 0.10546875/27., 0.10546875/27., 0.10546875/27., 0.10546875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 8.54296875/27., 8.54296875/27., 8.54296875/27., 8.54296875/27., - 0., 0., -1.6875/27., -1.6875/27., 0., 0., 0., 0., 15.1875/27., 15.1875/27., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - -0.52734375/27., -0.10546875/27., -0.10546875/27., -0.52734375/27., -1.58203125/27., 0.52734375/27., 0.94921875/27., 0.94921875/27., -1.58203125/27., 0.52734375/27., 4.74609375/27., 4.74609375/27., 14.23828125/27., -4.74609375/27., 14.23828125/27., -4.74609375/27., - 0., 0., 0., 0., -1.6875/27., 0., 0., 0., -1.6875/27., 0., 0., 0., 15.1875/27., 0., 15.1875/27., 0., - 0., 0., 0., 0., 0., 0., 0., -1.6875/27., 0., 0., 0., -1.6875/27., 0., 0., 15.1875/27., 15.1875/27., - -0.10546875/27., -0.10546875/27., -0.52734375/27., -0.52734375/27., 0.94921875/27., 0.94921875/27., 0.52734375/27., -1.58203125/27., 4.74609375/27., 4.74609375/27., 0.52734375/27., -1.58203125/27., -4.74609375/27., -4.74609375/27., 14.23828125/27., 14.23828125/27., - 0., 0., 1.6875/27., 8.4375/27., 0., 0., 0., 0., 25.3125/27., -8.4375/27., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., - 1.6875/27., 0., 0., 8.4375/27., 0., 0., 0., 0., 0., 0., -8.4375/27., 25.3125/27., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 1.6875/27., 0., 0., 0., 8.4375/27., 0., 0., 25.3125/27., -8.4375/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., - 0.52734375/27., 0.10546875/27., 0.52734375/27., 2.63671875/27., 1.58203125/27., -0.52734375/27., -0.52734375/27., 1.58203125/27., 7.91015625/27., -2.63671875/27., -2.63671875/27., 7.91015625/27., -7.91015625/27., 2.63671875/27., 23.73046875/27., -7.91015625/27., - 0., 0., 0., 0., 1.6875/27., 0., 0., 0., 8.4375/27., 0., 0., 0., -8.4375/27., 0., 25.3125/27., 0., - }; - -} // namespace FE_Q_2d - - -// embedding matrices - -template <> -const double * const -FE_Q<2>::Matrices::embedding[][GeometryInfo<2>::children_per_cell] = -{ - { FE_Q_2d::q1_into_q1_refined_0, FE_Q_2d::q1_into_q1_refined_1, - FE_Q_2d::q1_into_q1_refined_2, FE_Q_2d::q1_into_q1_refined_3 }, - { FE_Q_2d::q2_into_q2_refined_0, FE_Q_2d::q2_into_q2_refined_1, - FE_Q_2d::q2_into_q2_refined_2, FE_Q_2d::q2_into_q2_refined_3 }, - { FE_Q_2d::q3_into_q3_refined_0, FE_Q_2d::q3_into_q3_refined_1, - FE_Q_2d::q3_into_q3_refined_2, FE_Q_2d::q3_into_q3_refined_3 } -}; - - -template <> -const unsigned int -FE_Q<2>::Matrices::n_embedding_matrices - = sizeof(FE_Q<2>::Matrices::embedding) / - sizeof(FE_Q<2>::Matrices::embedding[0]); - - -// Constraint matrices taken from Wolfgangs old version -namespace FE_Q_2d -{ - static const double constraint_q1[] = - { - .5, .5 - }; - - static const double constraint_q2[] = - { - 0., 0., 1., - .375, -.125, .75, - -.125, .375, .75 - }; - - static const double constraint_q3[] = - { - -.0625, -.0625, .5625, .5625, - .3125, .0625, .9375, -.3125, - 0., 0., 1., 0., - 0., 0., 0., 1., - .0625, .3125, -.3125, 0.9375 - }; - - static const double constraint_q4[] = - { - 0., 0., 0., 1., 0., - 0.2734375, -0.0390625, 1.09375, -0.546875, 0.21875, - 0., 0., 1., 0., 0., - -0.0390625, 0.0234375, 0.46875, 0.703125, -0.15625, - 0.0234375, -0.0390625, -0.15625, 0.703125, 0.46875, - 0., 0., 0., 0., 1., - -0.0390625, 0.2734375, 0.21875, -0.546875, 1.09375 - }; -} - +// constraint matrices in 2d are now implemented by computing them on +// the fly for all polynomial degrees template <> const double * const -FE_Q<2>::Matrices::constraint_matrices[] = -{ - FE_Q_2d::constraint_q1, - FE_Q_2d::constraint_q2, - FE_Q_2d::constraint_q3, - FE_Q_2d::constraint_q4, -}; +FE_Q<2>::Matrices::constraint_matrices[] = {}; template <> const unsigned int -FE_Q<2>::Matrices::n_constraint_matrices - = sizeof(FE_Q<2>::Matrices::constraint_matrices) / - sizeof(FE_Q<2>::Matrices::constraint_matrices[0]); +FE_Q<2>::Matrices::n_constraint_matrices = 0; diff --git a/deal.II/deal.II/source/fe/fe_q_3d.cc b/deal.II/deal.II/source/fe/fe_q_3d.cc index c82c465701..6b0de5a7df 100644 --- a/deal.II/deal.II/source/fe/fe_q_3d.cc +++ b/deal.II/deal.II/source/fe/fe_q_3d.cc @@ -19,379 +19,6 @@ #include -namespace FE_Q_3d -{ - static const double q1_into_q1_refined_0[] = - { - 1., 0., 0., 0., 0., 0., 0., 0., - 13.5/27., 13.5/27., 0., 0., 0., 0., 0., 0., - 6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 0., 0., 0., - 13.5/27., 0., 0., 13.5/27., 0., 0., 0., 0., - 13.5/27., 0., 0., 0., 13.5/27., 0., 0., 0., - 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., - 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., - 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., - }; - - static const double q1_into_q1_refined_1[] = - { - 13.5/27., 13.5/27., 0., 0., 0., 0., 0., 0., - 0., 1., 0., 0., 0., 0., 0., 0., - 0., 13.5/27., 13.5/27., 0., 0., 0., 0., 0., - 6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 0., 0., 0., - 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., - 0., 13.5/27., 0., 0., 0., 13.5/27., 0., 0., - 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., - 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., - }; - - static const double q1_into_q1_refined_2[] = - { - 6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 0., 0., 0., - 0., 13.5/27., 13.5/27., 0., 0., 0., 0., 0., - 0., 0., 1., 0., 0., 0., 0., 0., - 0., 0., 13.5/27., 13.5/27., 0., 0., 0., 0., - 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., - 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., - 0., 0., 13.5/27., 0., 0., 0., 13.5/27., 0., - 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., - }; - - static const double q1_into_q1_refined_3[] = - { - 13.5/27., 0., 0., 13.5/27., 0., 0., 0., 0., - 6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 0., 0., 0., - 0., 0., 13.5/27., 13.5/27., 0., 0., 0., 0., - 0., 0., 0., 1., 0., 0., 0., 0., - 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., - 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., - 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., - 0., 0., 0., 13.5/27., 0., 0., 0., 13.5/27., - }; - - static const double q1_into_q1_refined_4[] = - { - 13.5/27., 0., 0., 0., 13.5/27., 0., 0., 0., - 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., - 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., - 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., - 0., 0., 0., 0., 1., 0., 0., 0., - 0., 0., 0., 0., 13.5/27., 13.5/27., 0., 0., - 0., 0., 0., 0., 6.75/27., 6.75/27., 6.75/27., 6.75/27., - 0., 0., 0., 0., 13.5/27., 0., 0., 13.5/27., - }; - - static const double q1_into_q1_refined_5[] = - { - 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., - 0., 13.5/27., 0., 0., 0., 13.5/27., 0., 0., - 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., - 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., - 0., 0., 0., 0., 13.5/27., 13.5/27., 0., 0., - 0., 0., 0., 0., 0., 1., 0., 0., - 0., 0., 0., 0., 0., 13.5/27., 13.5/27., 0., - 0., 0., 0., 0., 6.75/27., 6.75/27., 6.75/27., 6.75/27., - }; - - static const double q1_into_q1_refined_6[] = - { - 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., - 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., - 0., 0., 13.5/27., 0., 0., 0., 13.5/27., 0., - 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., - 0., 0., 0., 0., 6.75/27., 6.75/27., 6.75/27., 6.75/27., - 0., 0., 0., 0., 0., 13.5/27., 13.5/27., 0., - 0., 0., 0., 0., 0., 0., 1., 0., - 0., 0., 0., 0., 0., 0., 13.5/27., 13.5/27., - }; - - static const double q1_into_q1_refined_7[] = - { - 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., - 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., - 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., - 0., 0., 0., 13.5/27., 0., 0., 0., 13.5/27., - 0., 0., 0., 0., 13.5/27., 0., 0., 13.5/27., - 0., 0., 0., 0., 6.75/27., 6.75/27., 6.75/27., 6.75/27., - 0., 0., 0., 0., 0., 0., 13.5/27., 13.5/27., - 0., 0., 0., 0., 0., 0., 0., 1., - }; - - static const double q2_into_q2_refined_0[] = - { - 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., - 10.125/27., -3.375/27., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., - 10.125/27., 0., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., -3.375/27., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 0., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., 0., -3.375/27., 0., 0., 0., 0., 0., 20.25/27., 0., - 10.125/27., 0., 0., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., 0., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., -3.375/27., 0., 0., 0., 0., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., 0., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., - 3.796875/27., -1.265625/27., 0.421875/27., -1.265625/27., 0., 0., 0., 0., 7.59375/27., -2.53125/27., -2.53125/27., 7.59375/27., 0., 0., 0., 0., 0., 0., 0., 0., 15.1875/27., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 3.796875/27., -1.265625/27., 0.421875/27., -1.265625/27., 0., 0., 7.59375/27., -2.53125/27., -2.53125/27., 7.59375/27., 15.1875/27., - 3.796875/27., -1.265625/27., 0., 0., -1.265625/27., 0.421875/27., 0., 0., 7.59375/27., 0., 0., 0., -2.53125/27., 0., 0., 0., 7.59375/27., -2.53125/27., 0., 0., 0., 0., 15.1875/27., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 3.796875/27., 0., -1.265625/27., 0., -1.265625/27., 0., 0.421875/27., 0., 0., 0., 0., 0., 7.59375/27., -2.53125/27., 7.59375/27., 0., -2.53125/27., 0., 15.1875/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1.265625/27., 0., 3.796875/27., 0., 0.421875/27., 0., -1.265625/27., 0., 0., 0., 0., 7.59375/27., -2.53125/27., 0., -2.53125/27., 0., 7.59375/27., 15.1875/27., - 3.796875/27., 0., 0., -1.265625/27., -1.265625/27., 0., 0., 0.421875/27., 0., 0., 0., 7.59375/27., 0., 0., 0., -2.53125/27., 7.59375/27., 0., 0., -2.53125/27., 0., 0., 0., 0., 0., 15.1875/27., 0., - 1.423828125/27., -0.474609375/27., 0.158203125/27., -0.474609375/27., -0.474609375/27., 0.158203125/27., -0.052734375/27., 0.158203125/27., 2.84765625/27., -0.94921875/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., 5.6953125/27., -1.8984375/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., 5.6953125/27., 11.390625/27., - }; - - static const double q2_into_q2_refined_1[] = - { - 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., - -3.375/27., 10.125/27., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 10.125/27., -3.375/27., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 10.125/27., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., -3.375/27., 0., 0., 0., 0., 20.25/27., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 0., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., 0., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., - 0., 10.125/27., 0., 0., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., 0., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., -3.375/27., 0., 0., 0., 0., 20.25/27., - -1.265625/27., 3.796875/27., -1.265625/27., 0.421875/27., 0., 0., 0., 0., 7.59375/27., 7.59375/27., -2.53125/27., -2.53125/27., 0., 0., 0., 0., 0., 0., 0., 0., 15.1875/27., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1.265625/27., 3.796875/27., -1.265625/27., 0.421875/27., 0., 0., 7.59375/27., 7.59375/27., -2.53125/27., -2.53125/27., 15.1875/27., - -1.265625/27., 3.796875/27., 0., 0., 0.421875/27., -1.265625/27., 0., 0., 7.59375/27., 0., 0., 0., -2.53125/27., 0., 0., 0., -2.53125/27., 7.59375/27., 0., 0., 0., 0., 15.1875/27., 0., 0., 0., 0., - 0., 3.796875/27., -1.265625/27., 0., 0., -1.265625/27., 0.421875/27., 0., 0., 7.59375/27., 0., 0., 0., -2.53125/27., 0., 0., 0., 7.59375/27., -2.53125/27., 0., 0., 0., 0., 15.1875/27., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 3.796875/27., 0., -1.265625/27., 0., -1.265625/27., 0., 0.421875/27., 0., 0., 0., 0., 7.59375/27., -2.53125/27., 0., 7.59375/27., 0., -2.53125/27., 15.1875/27., - 0., 0., 0., 0., 0., 0., 0., 0., 3.796875/27., 0., -1.265625/27., 0., -1.265625/27., 0., 0.421875/27., 0., 0., 0., 0., 0., 7.59375/27., -2.53125/27., 7.59375/27., 0., -2.53125/27., 0., 15.1875/27., - -0.474609375/27., 1.423828125/27., -0.474609375/27., 0.158203125/27., 0.158203125/27., -0.474609375/27., 0.158203125/27., -0.052734375/27., 2.84765625/27., 2.84765625/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 0.31640625/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., 5.6953125/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., 11.390625/27., - }; - - static const double q2_into_q2_refined_2[] = - { - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., - 0., -3.375/27., 10.125/27., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 10.125/27., -3.375/27., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 10.125/27., 0., 0., 0., 0., 20.25/27., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., -3.375/27., 0., 0., 0., 0., 20.25/27., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 0., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., -3.375/27., 0., 0., 0., 0., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., 0., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., - 0., 0., 10.125/27., 0., 0., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., 0., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., - 0.421875/27., -1.265625/27., 3.796875/27., -1.265625/27., 0., 0., 0., 0., -2.53125/27., 7.59375/27., 7.59375/27., -2.53125/27., 0., 0., 0., 0., 0., 0., 0., 0., 15.1875/27., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.421875/27., -1.265625/27., 3.796875/27., -1.265625/27., 0., 0., -2.53125/27., 7.59375/27., 7.59375/27., -2.53125/27., 15.1875/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 3.796875/27., 0., -1.265625/27., 0., -1.265625/27., 0., 0.421875/27., 0., 0., 0., 0., 7.59375/27., -2.53125/27., 0., 7.59375/27., 0., -2.53125/27., 15.1875/27., - 0., -1.265625/27., 3.796875/27., 0., 0., 0.421875/27., -1.265625/27., 0., 0., 7.59375/27., 0., 0., 0., -2.53125/27., 0., 0., 0., -2.53125/27., 7.59375/27., 0., 0., 0., 0., 15.1875/27., 0., 0., 0., - 0., 0., 3.796875/27., -1.265625/27., 0., 0., -1.265625/27., 0.421875/27., 0., 0., 7.59375/27., 0., 0., 0., -2.53125/27., 0., 0., 0., 7.59375/27., -2.53125/27., 0., 0., 0., 0., 15.1875/27., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., -1.265625/27., 0., 3.796875/27., 0., 0.421875/27., 0., -1.265625/27., 0., 0., 0., 0., 0., 7.59375/27., -2.53125/27., -2.53125/27., 0., 7.59375/27., 0., 15.1875/27., - 0.158203125/27., -0.474609375/27., 1.423828125/27., -0.474609375/27., -0.052734375/27., 0.158203125/27., -0.474609375/27., 0.158203125/27., -0.94921875/27., 2.84765625/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., -0.94921875/27., 0.31640625/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., -1.8984375/27., 11.390625/27., - }; - - static const double q2_into_q2_refined_3[] = - { - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., - 0., 0., -3.375/27., 10.125/27., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - -3.375/27., 0., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 0., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 10.125/27., 0., 0., 0., 0., 20.25/27., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 0., 10.125/27., 0., 0., 0., 0., 0., 20.25/27., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., 0., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., -3.375/27., 0., 0., 0., 0., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., 0., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., - 0., 0., 0., 10.125/27., 0., 0., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., - -1.265625/27., 0.421875/27., -1.265625/27., 3.796875/27., 0., 0., 0., 0., -2.53125/27., -2.53125/27., 7.59375/27., 7.59375/27., 0., 0., 0., 0., 0., 0., 0., 0., 15.1875/27., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1.265625/27., 0.421875/27., -1.265625/27., 3.796875/27., 0., 0., -2.53125/27., -2.53125/27., 7.59375/27., 7.59375/27., 15.1875/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1.265625/27., 0., 3.796875/27., 0., 0.421875/27., 0., -1.265625/27., 0., 0., 0., 0., 7.59375/27., -2.53125/27., 0., -2.53125/27., 0., 7.59375/27., 15.1875/27., - 0., 0., 0., 0., 0., 0., 0., 0., -1.265625/27., 0., 3.796875/27., 0., 0.421875/27., 0., -1.265625/27., 0., 0., 0., 0., 0., 7.59375/27., -2.53125/27., -2.53125/27., 0., 7.59375/27., 0., 15.1875/27., - 0., 0., -1.265625/27., 3.796875/27., 0., 0., 0.421875/27., -1.265625/27., 0., 0., 7.59375/27., 0., 0., 0., -2.53125/27., 0., 0., 0., -2.53125/27., 7.59375/27., 0., 0., 0., 0., 15.1875/27., 0., 0., - -1.265625/27., 0., 0., 3.796875/27., 0.421875/27., 0., 0., -1.265625/27., 0., 0., 0., 7.59375/27., 0., 0., 0., -2.53125/27., -2.53125/27., 0., 0., 7.59375/27., 0., 0., 0., 0., 0., 15.1875/27., 0., - -0.474609375/27., 0.158203125/27., -0.474609375/27., 1.423828125/27., 0.158203125/27., -0.052734375/27., 0.158203125/27., -0.474609375/27., -0.94921875/27., -0.94921875/27., 2.84765625/27., 2.84765625/27., 0.31640625/27., 0.31640625/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., 11.390625/27., - }; - - static const double q2_into_q2_refined_4[] = - { - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., - 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., -3.375/27., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 0., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., 0., -3.375/27., 0., 0., 0., 0., 0., 20.25/27., 0., - 0., 0., 0., 0., 10.125/27., -3.375/27., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 10.125/27., 0., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - -3.375/27., 0., 0., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 0., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 10.125/27., 0., 0., 0., 0., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 0., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 3.796875/27., -1.265625/27., 0.421875/27., -1.265625/27., 0., 0., 7.59375/27., -2.53125/27., -2.53125/27., 7.59375/27., 15.1875/27., - 0., 0., 0., 0., 3.796875/27., -1.265625/27., 0.421875/27., -1.265625/27., 0., 0., 0., 0., 7.59375/27., -2.53125/27., -2.53125/27., 7.59375/27., 0., 0., 0., 0., 0., 15.1875/27., 0., 0., 0., 0., 0., - -1.265625/27., 0.421875/27., 0., 0., 3.796875/27., -1.265625/27., 0., 0., -2.53125/27., 0., 0., 0., 7.59375/27., 0., 0., 0., 7.59375/27., -2.53125/27., 0., 0., 0., 0., 15.1875/27., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., -1.265625/27., 0., 0.421875/27., 0., 3.796875/27., 0., -1.265625/27., 0., 0., 0., 0., 0., -2.53125/27., 7.59375/27., 7.59375/27., 0., -2.53125/27., 0., 15.1875/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.421875/27., 0., -1.265625/27., 0., -1.265625/27., 0., 3.796875/27., 0., 0., 0., 0., -2.53125/27., 7.59375/27., 0., -2.53125/27., 0., 7.59375/27., 15.1875/27., - -1.265625/27., 0., 0., 0.421875/27., 3.796875/27., 0., 0., -1.265625/27., 0., 0., 0., -2.53125/27., 0., 0., 0., 7.59375/27., 7.59375/27., 0., 0., -2.53125/27., 0., 0., 0., 0., 0., 15.1875/27., 0., - -0.474609375/27., 0.158203125/27., -0.052734375/27., 0.158203125/27., 1.423828125/27., -0.474609375/27., 0.158203125/27., -0.474609375/27., -0.94921875/27., 0.31640625/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., -0.94921875/27., 2.84765625/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., 5.6953125/27., 11.390625/27., - }; - - static const double q2_into_q2_refined_5[] = - { - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 10.125/27., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., -3.375/27., 0., 0., 0., 0., 20.25/27., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 0., 20.25/27., - 0., 0., 0., 0., -3.375/27., 10.125/27., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 10.125/27., -3.375/27., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 0., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., - 0., -3.375/27., 0., 0., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 0., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 10.125/27., 0., 0., 0., 0., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1.265625/27., 3.796875/27., -1.265625/27., 0.421875/27., 0., 0., 7.59375/27., 7.59375/27., -2.53125/27., -2.53125/27., 15.1875/27., - 0., 0., 0., 0., -1.265625/27., 3.796875/27., -1.265625/27., 0.421875/27., 0., 0., 0., 0., 7.59375/27., 7.59375/27., -2.53125/27., -2.53125/27., 0., 0., 0., 0., 0., 15.1875/27., 0., 0., 0., 0., 0., - 0.421875/27., -1.265625/27., 0., 0., -1.265625/27., 3.796875/27., 0., 0., -2.53125/27., 0., 0., 0., 7.59375/27., 0., 0., 0., -2.53125/27., 7.59375/27., 0., 0., 0., 0., 15.1875/27., 0., 0., 0., 0., - 0., -1.265625/27., 0.421875/27., 0., 0., 3.796875/27., -1.265625/27., 0., 0., -2.53125/27., 0., 0., 0., 7.59375/27., 0., 0., 0., 7.59375/27., -2.53125/27., 0., 0., 0., 0., 15.1875/27., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1.265625/27., 0., 0.421875/27., 0., 3.796875/27., 0., -1.265625/27., 0., 0., 0., 0., -2.53125/27., 7.59375/27., 0., 7.59375/27., 0., -2.53125/27., 15.1875/27., - 0., 0., 0., 0., 0., 0., 0., 0., -1.265625/27., 0., 0.421875/27., 0., 3.796875/27., 0., -1.265625/27., 0., 0., 0., 0., 0., -2.53125/27., 7.59375/27., 7.59375/27., 0., -2.53125/27., 0., 15.1875/27., - 0.158203125/27., -0.474609375/27., 0.158203125/27., -0.052734375/27., -0.474609375/27., 1.423828125/27., -0.474609375/27., 0.158203125/27., -0.94921875/27., -0.94921875/27., 0.31640625/27., 0.31640625/27., 2.84765625/27., 2.84765625/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., 11.390625/27., - }; - - static const double q2_into_q2_refined_6[] = - { - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 10.125/27., 0., 0., 0., 0., 20.25/27., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., -3.375/27., 0., 0., 0., 0., 20.25/27., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 0., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 10.125/27., 0., -3.375/27., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., -3.375/27., 10.125/27., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 10.125/27., -3.375/27., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 10.125/27., 0., 0., 0., 0., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 0., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., - 0., 0., -3.375/27., 0., 0., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 0., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.421875/27., -1.265625/27., 3.796875/27., -1.265625/27., 0., 0., -2.53125/27., 7.59375/27., 7.59375/27., -2.53125/27., 15.1875/27., - 0., 0., 0., 0., 0.421875/27., -1.265625/27., 3.796875/27., -1.265625/27., 0., 0., 0., 0., -2.53125/27., 7.59375/27., 7.59375/27., -2.53125/27., 0., 0., 0., 0., 0., 15.1875/27., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1.265625/27., 0., 0.421875/27., 0., 3.796875/27., 0., -1.265625/27., 0., 0., 0., 0., -2.53125/27., 7.59375/27., 0., 7.59375/27., 0., -2.53125/27., 15.1875/27., - 0., 0.421875/27., -1.265625/27., 0., 0., -1.265625/27., 3.796875/27., 0., 0., -2.53125/27., 0., 0., 0., 7.59375/27., 0., 0., 0., -2.53125/27., 7.59375/27., 0., 0., 0., 0., 15.1875/27., 0., 0., 0., - 0., 0., -1.265625/27., 0.421875/27., 0., 0., 3.796875/27., -1.265625/27., 0., 0., -2.53125/27., 0., 0., 0., 7.59375/27., 0., 0., 0., 7.59375/27., -2.53125/27., 0., 0., 0., 0., 15.1875/27., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0.421875/27., 0., -1.265625/27., 0., -1.265625/27., 0., 3.796875/27., 0., 0., 0., 0., 0., -2.53125/27., 7.59375/27., -2.53125/27., 0., 7.59375/27., 0., 15.1875/27., - -0.052734375/27., 0.158203125/27., -0.474609375/27., 0.158203125/27., 0.158203125/27., -0.474609375/27., 1.423828125/27., -0.474609375/27., 0.31640625/27., -0.94921875/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., -1.8984375/27., 5.6953125/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., -1.8984375/27., 11.390625/27., - }; - - static const double q2_into_q2_refined_7[] = - { - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 0., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 10.125/27., 0., 0., 0., 0., 20.25/27., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 0., 10.125/27., 0., 0., 0., 0., 0., 20.25/27., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., -3.375/27., 10.125/27., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., -3.375/27., 0., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 0., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 10.125/27., 0., 0., 0., 0., 20.25/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -3.375/27., 0., 0., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., - 0., 0., 0., -3.375/27., 0., 0., 0., 10.125/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 20.25/27., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1.265625/27., 0.421875/27., -1.265625/27., 3.796875/27., 0., 0., -2.53125/27., -2.53125/27., 7.59375/27., 7.59375/27., 15.1875/27., - 0., 0., 0., 0., -1.265625/27., 0.421875/27., -1.265625/27., 3.796875/27., 0., 0., 0., 0., -2.53125/27., -2.53125/27., 7.59375/27., 7.59375/27., 0., 0., 0., 0., 0., 15.1875/27., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.421875/27., 0., -1.265625/27., 0., -1.265625/27., 0., 3.796875/27., 0., 0., 0., 0., -2.53125/27., 7.59375/27., 0., -2.53125/27., 0., 7.59375/27., 15.1875/27., - 0., 0., 0., 0., 0., 0., 0., 0., 0.421875/27., 0., -1.265625/27., 0., -1.265625/27., 0., 3.796875/27., 0., 0., 0., 0., 0., -2.53125/27., 7.59375/27., -2.53125/27., 0., 7.59375/27., 0., 15.1875/27., - 0., 0., 0.421875/27., -1.265625/27., 0., 0., -1.265625/27., 3.796875/27., 0., 0., -2.53125/27., 0., 0., 0., 7.59375/27., 0., 0., 0., -2.53125/27., 7.59375/27., 0., 0., 0., 0., 15.1875/27., 0., 0., - 0.421875/27., 0., 0., -1.265625/27., -1.265625/27., 0., 0., 3.796875/27., 0., 0., 0., -2.53125/27., 0., 0., 0., 7.59375/27., -2.53125/27., 0., 0., 7.59375/27., 0., 0., 0., 0., 0., 15.1875/27., 0., - 0.158203125/27., -0.052734375/27., 0.158203125/27., -0.474609375/27., -0.474609375/27., 0.158203125/27., -0.474609375/27., 1.423828125/27., 0.31640625/27., 0.31640625/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 2.84765625/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., -1.8984375/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., 11.390625/27., - }; -} // namespace FE_Q_3d - - -// embedding matrices - -template <> -const double * const -FE_Q<3>::Matrices::embedding[][GeometryInfo<3>::children_per_cell] = -{ - { FE_Q_3d::q1_into_q1_refined_0, FE_Q_3d::q1_into_q1_refined_1, - FE_Q_3d::q1_into_q1_refined_2, FE_Q_3d::q1_into_q1_refined_3, - FE_Q_3d::q1_into_q1_refined_4, FE_Q_3d::q1_into_q1_refined_5, - FE_Q_3d::q1_into_q1_refined_6, FE_Q_3d::q1_into_q1_refined_7 }, - { FE_Q_3d::q2_into_q2_refined_0, FE_Q_3d::q2_into_q2_refined_1, - FE_Q_3d::q2_into_q2_refined_2, FE_Q_3d::q2_into_q2_refined_3, - FE_Q_3d::q2_into_q2_refined_4, FE_Q_3d::q2_into_q2_refined_5, - FE_Q_3d::q2_into_q2_refined_6, FE_Q_3d::q2_into_q2_refined_7 } -}; - - -template <> -const unsigned int -FE_Q<3>::Matrices::n_embedding_matrices - = sizeof(FE_Q<3>::Matrices::embedding) / - sizeof(FE_Q<3>::Matrices::embedding[0]); - - - // Constraint matrices taken from Wolfgangs old version namespace FE_Q_3d diff --git a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc index 9c45ac62dd..5c924144cc 100644 --- a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc +++ b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc @@ -23,6 +23,12 @@ #include #include +#ifdef HAVE_STD_STRINGSTREAM +# include +#else +# include +#endif + namespace { @@ -92,6 +98,26 @@ FE_Q_Hierarchical::FE_Q_Hierarchical (const unsigned int degree) +template +std::string +FE_Q_Hierarchical::get_name () const +{ +#ifdef HAVE_STD_STRINGSTREAM + std::ostringstream namebuf; +#else + std::ostrstream namebuf; +#endif + + namebuf << "FE_Q_Hierarchical<" << dim << ">(" << degree << ")"; + +#ifndef HAVE_STD_STRINGSTREAM + namebuf << std::ends; +#endif + return namebuf.str(); +} + + + template FiniteElement * FE_Q_Hierarchical::clone() const diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc new file mode 100644 index 0000000000..d14c6fba94 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc @@ -0,0 +1,1983 @@ +//---------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2003 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 + +#ifdef HAVE_STD_STRINGSTREAM +# include +#else +# include +#endif + + +// namespace for some functions that are used in this file. they are +// specific to numbering conventions used for the FE_RT element, and +// are thus not very interesting to the outside world +namespace +{ + // auxiliary type to allow for some + // kind of explicit template + // specialization of the following + // functions + template struct int2type {}; + + + // generate the j-th out of a total + // of N points on the unit square + // in 2d. N needs not be a square + // number, but must be the product + // of two integers + // + // there is one complication: we + // want to generate interpolation + // points on the unit square for + // the shape functions for this + // element, but for that we need to + // make sure that these + // interpolation points make the + // resulting matrix rows linearly + // independent. this is a problem + // since we have anisotropic + // polynomials, so for example for + // the lowest order elements, we + // have as polynomials in for the + // x-component of the shape + // functions only "x" and "1-x", + // i.e. no y-dependence. if we + // select as interpolation points + // the points (.5,0) and (.5,1), + // we're hosed! + // + // thus, the third parameter gives + // the coordinate direction in + // which the polynomial degree is + // highest. we use this to select + // interpolation points primarily + // in this direction then + Point<2> generate_unit_point (const unsigned int j, + const unsigned int N, + const unsigned int d, + const int2type<2> &) + { + Assert (d<2, ExcInternalError()); + + // factorize N int N1*N2. note + // that we always have N1<=N2, + // since the square root is + // rounded down + const unsigned int N1 = static_cast(std::sqrt(1.*N)); + const unsigned int N2 = N/N1; + Assert (N1*N2 == N, ExcInternalError()); + + const unsigned int Nx = (d==0 ? N2 : N1), + Ny = (d==1 ? N2 : N1); + + return Point<2> (Nx == 1 ? .5 : 1.*(j%Nx)/(Nx-1), + Ny == 1 ? .5 : 1.*(j/Nx)/(Ny-1)); + } + + + // generate the j-th out of a total + // of N points on the unit cube + // in 3d. N needs not be a cube + // number, but must be the product + // of three integers + // + // the same applies as above for + // the meaning of the parameter "d" + Point<3> generate_unit_point (const unsigned int /*j*/, + const unsigned int N, + const unsigned int d, + const int2type<3> &) + { + Assert (d<3, ExcInternalError()); + + const unsigned int N1 = static_cast(std::pow(1.*N, 1./3.)); + const unsigned int N2 = static_cast(std::sqrt(1.*N/N1)); + const unsigned int N3 = N/(N1*N2); + Assert (N1*N2*N3 == N, ExcInternalError()); + + Assert (false, ExcNotImplemented()); + + return Point<3> (); + } + +} + + + +template +FE_RaviartThomas::FE_RaviartThomas (const unsigned int degree) + : + FiniteElement (FiniteElementData(get_dpo_vector(degree), + dim), + get_ria_vector (degree), + std::vector >(FiniteElementData(get_dpo_vector(degree),dim).dofs_per_cell, + std::vector(dim,true))), + degree(degree), + polynomials (create_polynomials(degree)), + renumber (compute_renumber(degree)) +{ + Assert (dim >= 2, ExcNotUsefulInThisDimension()); + + // check formula (III.3.22) in the + // book by Brezzi & Fortin about + // the number of degrees of freedom + // per cell + Assert (((dim==2) && + (this->dofs_per_cell == 2*(degree+1)*(degree+2))) + || + ((dim==3) && + (this->dofs_per_cell == 3*(degree+1)*(degree+1)*(degree+2))), + ExcInternalError()); + Assert (renumber.size() == this->dofs_per_cell, + ExcInternalError()); + + // initialize the various matrices + initialize_constraints (); + initialize_embedding (); + initialize_restriction (); + + // finally fill in support points + // on cell and face + initialize_unit_support_points (); + initialize_unit_face_support_points (); + + // then make + // system_to_component_table + // invalid, since this has no + // meaning for the present element + std::vector > tmp1, tmp2; + this->system_to_component_table.swap (tmp1); + this->face_system_to_component_table.swap (tmp2); +} + + + +template +std::string +FE_RaviartThomas::get_name () const +{ +#ifdef HAVE_STD_STRINGSTREAM + std::ostringstream namebuf; +#else + std::ostrstream namebuf; +#endif + + namebuf << "FE_RaviartThomas<" << dim << ">(" << degree << ")"; + +#ifndef HAVE_STD_STRINGSTREAM + namebuf << std::ends; +#endif + return namebuf.str(); +} + + + +template +FiniteElement * +FE_RaviartThomas::clone() const +{ + return new FE_RaviartThomas(degree); +} + + +template +double +FE_RaviartThomas::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 < dim, ExcIndexRange (component, 0, dim)); + + // check whether this shape + // function has a contribution in + // this component at all, and if so + // delegate to the respective + // polynomial + if (component == renumber[i].first) + return polynomials[component].compute_value(renumber[i].second, p); + else + return 0; +} + + + +template +Tensor<1,dim> +FE_RaviartThomas::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 < dim, ExcIndexRange (component, 0, dim)); + + // check whether this shape + // function has a contribution in + // this component at all, and if so + // delegate to the respective + // polynomial + if (component == renumber[i].first) + return polynomials[component].compute_grad(renumber[i].second, p); + else + return Tensor<1,dim>(); +} + + + +template +Tensor<2,dim> +FE_RaviartThomas::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 < dim, ExcIndexRange (component, 0, dim)); + + // check whether this shape + // function has a contribution in + // this component at all, and if so + // delegate to the respective + // polynomial + if (component == renumber[i].first) + return polynomials[component].compute_grad_grad(renumber[i].second, p); + else + return Tensor<2,dim>(); +} + + + +#if deal_II_dimension == 1 + +template <> +void +FE_RaviartThomas<1>:: +get_interpolation_matrix (const FiniteElementBase<1> &, + FullMatrix &) const +{ + Assert (false, ExcNotUsefulInThisDimension()); +} + +#endif + + +template +void +FE_RaviartThomas:: +get_interpolation_matrix (const FiniteElementBase &x_source_fe, + FullMatrix &interpolation_matrix) const +{ + // this is only implemented, if the + // source FE is also a + // Raviart-Thomas element, + // otherwise throw an exception, as + // the documentation says + AssertThrow ((x_source_fe.get_name().find ("FE_RaviartThomas<") == 0) + || + (dynamic_cast*>(&x_source_fe) != 0), + typename FiniteElementBase:: + ExcInterpolationNotImplemented()); + + // ok, source is a RT element, so + // we will be able to do the work + const FE_RaviartThomas &source_fe + = dynamic_cast&>(x_source_fe); + + Assert (interpolation_matrix.m() == this->dofs_per_cell, + ExcDimensionMismatch (interpolation_matrix.m(), + this->dofs_per_cell)); + Assert (interpolation_matrix.n() == source_fe.dofs_per_cell, + ExcDimensionMismatch (interpolation_matrix.m(), + source_fe.dofs_per_cell)); + + + // compute the interpolation + // matrices in much the same way as + // we do for the embedding matrices + // from mother to child. + const unsigned int dofs_per_coordinate = this->dofs_per_cell/dim; + Assert (dofs_per_coordinate*dim == this->dofs_per_cell, + ExcInternalError()); + for (unsigned int d=0; d cell_interpolation (dofs_per_coordinate, + dofs_per_coordinate); + FullMatrix source_interpolation (dofs_per_coordinate, + source_dofs_per_coordinate); + FullMatrix tmp (dofs_per_coordinate, + source_dofs_per_coordinate); + for (unsigned int d=0; d p = generate_unit_point (j, dofs_per_coordinate, + d, int2type()); + for (unsigned int i=0; idofs_per_cell; ++i) + if (renumber[i].first == d) + for (unsigned int j=0; j 1e-15) + interpolation_matrix(i,j) = tmp(renumber[i].second, + source_fe.renumber[j].second); + } + + // if this were a Lagrange + // interpolation element, we could + // make sure that the row sum of + // each of the matrices is 1 at + // this point. note that this won't + // work here, since we are working + // with hierarchical elements for + // which the shape functions don't + // sum up to 1 + // + // however, we can make sure that + // only components couple that have + // the same vector component + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int j=0; j +void +FE_RaviartThomas<1>::initialize_constraints () +{ + Assert (false, ExcNotUsefulInThisDimension()); +} + +#endif + +#if deal_II_dimension == 2 + +template <> +void +FE_RaviartThomas<2>::initialize_constraints () +{ + const unsigned int dim = 2; + + this->interface_constraints. + TableBase<2,double>::reinit (this->interface_constraints_size()); + + // this case is too easy, so + // special case it + if (degree == 0) + { + this->interface_constraints(0,0) = this->interface_constraints(1,0) = .5; + return; + } + + // for higher orders of the + // Raviart-Thomas element: + + // restricted to each face, the + // normal component of the shape + // functions is an element of P_{k} + // (in 2d), or Q_{k} (in 3d), where + // k is the degree of the element + // + // from this, we interpolate + // between mother and cell + // face. this is slightly + // complicated by the fact that we + // don't use Lagrange interpolation + // polynomials, but rather + // hierarchical polynomials, so we + // can't just use point + // interpolation. what we do + // instead is to evaluate at a + // number of points and then invert + // the interpolation matrix + + // mathematically speaking, this + // works in the following way: on + // each subface, we want that + // finite element solututions from + // both sides coincide. i.e. if a + // and b are expansion coefficients + // for the shape functions from + // both sides, we seek a relation + // between x and y such that + // sum_i a_i phi^c_i(x) + // == sum_j b_j phi_j(x) + // for all points x on the + // interface. here, phi^c_i are the + // shape functions on the small + // cell on one side of the face, + // and phi_j those on the big cell + // on the other side. To get this + // relation, it suffices to look at + // a sufficient number of points + // for which this has to hold. if + // there are n functions, then we + // need n evaluation points, and we + // choose them equidistantly. + // + // what one then gets is a matrix + // system + // a A == b B + // where + // A_ij = phi^c_i(x_j) + // B_ij = phi_i(x_j) + // and the relation we are looking for + // is + // a = (A^T)^-1 B^T b + // + // below, we build up these + // matrices, but rather than + // transposing them after the + // fact, we do so while building + // them. A will be + // subface_interpolation, B will be + // face_interpolation. note that we + // build up these matrices for all + // faces at once, rather than + // considering them separately. the + // reason is that we finally will + // want to have them in this order + // anyway, as this is the format we + // need inside deal.II + const std::vector > + face_polynomials (Polynomials::Hierarchical:: + generate_complete_basis (degree)); + Assert (face_polynomials.size() == this->dofs_per_face, ExcInternalError()); + + FullMatrix face_interpolation (2*this->dofs_per_face, this->dofs_per_face); + FullMatrix subface_interpolation (2*this->dofs_per_face, 2*this->dofs_per_face); + + // generate the matrix for the + // evaluation points on the big + // face, and the corresponding + // points in the coordinate system + // of the small face. order the + // shape functions in the same way + // we want to have them in the + // final matrix. extend shape + // functions on the small faces by + // zero to the other face on which + // they are not defined (we do this + // by simply not considering these + // entries in the matrix) + // + // note the agreeable fact that for + // this element, all the shape + // functions we presently care for + // are face-based (i.e. not vertex + // shape functions); thus, for this + // element, we can skip the + // annoying index shifting for the + // constraints matrix due to its + // weird format + for (unsigned int subface=0; subface::subfaces_per_face; ++subface) + for (unsigned int i=0; idofs_per_face; ++i) + { + const double p_face (1.*i/degree/2 + (subface == 0 ? 0. : .5)); + const double p_subface (1.*i/degree); + + for (unsigned int j=0; jdofs_per_face; ++j) + { + face_interpolation(subface*this->dofs_per_face+i, + j) + = face_polynomials[j].value(p_face); + subface_interpolation(subface*this->dofs_per_face+i, + subface*this->dofs_per_face+j) + = face_polynomials[j].value(p_subface); + } + } + + subface_interpolation.gauss_jordan (); + subface_interpolation.mmult (this->interface_constraints, + face_interpolation); + + // there is one additional thing to + // be considered: since the shape + // functions on the real cell + // contain the Jacobian (actually, + // the determinant of the inverse), + // there is an additional factor of + // 2 when going from the big to the + // small cell: + this->interface_constraints *= 1./2; + + // finally: constraints become + // really messy if the matrix in + // question has some entries that + // are almost zero, but not + // quite. this will happen in the + // above procedure due to + // round-off. let us simply delete + // these entries + for (unsigned int i=0; iinterface_constraints.m(); ++i) + for (unsigned int j=0; jinterface_constraints.n(); ++j) + if (fabs(this->interface_constraints(i,j)) < 1e-14) + this->interface_constraints(i,j) = 0.; +} + +#endif + +#if deal_II_dimension == 3 + +template <> +void +FE_RaviartThomas<3>::initialize_constraints () +{ + Assert (false, ExcNotImplemented()); +} + +#endif + + +#if deal_II_dimension == 1 + +template <> +void +FE_RaviartThomas<1>::initialize_embedding () +{ + Assert (false, ExcNotUsefulInThisDimension()); +} + +#endif + + +template +void +FE_RaviartThomas::initialize_embedding () +{ + // compute the interpolation + // matrices in much the same way as + // we do for the constraints. it's + // actually simpler here, since we + // don't have this weird + // renumbering stuff going on + // + // it is, however, slightly + // complicated by the fact that we + // have vector-valued elements + // here, so we do all the stuff for + // the degrees of freedom + // corresponding to each coordinate + // direction separately + const unsigned int dofs_per_coordinate = this->dofs_per_cell/dim; + Assert (dofs_per_coordinate*dim == this->dofs_per_cell, + ExcInternalError()); + for (unsigned int d=0; d cell_interpolation (dofs_per_coordinate, + dofs_per_coordinate); + FullMatrix subcell_interpolation (dofs_per_coordinate, + dofs_per_coordinate); + FullMatrix tmp (dofs_per_coordinate, + dofs_per_coordinate); + for (unsigned int child=0; child::children_per_cell; ++child) + this->prolongation[child].reinit (this->dofs_per_cell, + this->dofs_per_cell); + for (unsigned int child=0; child::children_per_cell; ++child) + for (unsigned int d=0; d p_subcell = generate_unit_point (j, dofs_per_coordinate, + d, int2type()); + const Point p_cell = + GeometryInfo::child_to_cell_coordinates (p_subcell, child); + + for (unsigned int i=0; idofs_per_cell; ++i) + if (renumber[i].first == d) + for (unsigned int j=0; jdofs_per_cell; ++j) + if (renumber[j].first == d) + if (std::fabs(tmp(renumber[i].second,renumber[j].second)) > 1e-15) + this->prolongation[child](i,j) = tmp(renumber[i].second, + renumber[j].second); + } + + // if this were a Lagrange + // interpolation element, we could + // make sure that the row sum of + // each of the matrices is 1 at + // this point. note that this won't + // work here, since we are working + // with hierarchical elements for + // which the shape functions don't + // sum up to 1 + // + // however, we can make sure that + // only components couple that have + // the same vector component + for (unsigned int child=0; child::children_per_cell; ++child) + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int j=0; jdofs_per_cell; ++j) + Assert ((this->prolongation[child](i,j) == 0.) || + (renumber[i].first == renumber[j].first), + ExcInternalError()); + + + // there is one additional thing to + // be considered: since the shape + // functions on the real cell + // contain the Jacobian (actually, + // the determinant of the inverse), + // there is an additional factor of + // 2 when going from the big to the + // small cell: + for (unsigned int child=0; child::children_per_cell; ++child) + this->prolongation[child] *= 1./2; +} + + +#if deal_II_dimension == 1 + +template <> +void +FE_RaviartThomas<1>::initialize_restriction () +{} + +#endif + + +#if deal_II_dimension == 2 + +template <> +void +FE_RaviartThomas<2>::initialize_restriction () +{ + const unsigned int dim = 2; + switch (degree) + { + case 0: + { + // this is a strange element, + // since it is both additive + // and then it is also + // not. ideally, we would + // like to have the value of + // the shape function on the + // coarse line to be the mean + // value of that on the two + // child ones. thus, one + // should make it + // additive. however, + // additivity only works if + // an element does not have + // any continuity + // requirements, since + // otherwise degrees of + // freedom are shared between + // adjacent elements, and + // when we make the element + // additive, that would mean + // that we end up adding up + // contributions not only + // from the child cells of + // this cell, but also from + // the child cells of the + // neighbor, and since we + // cannot know whether there + // even exists a neighbor we + // cannot simply make the + // element additive. + // + // so, until someone comes + // along with a better + // alternative, we do the + // following: make the + // element non-additive, and + // simply pick the value of + // one of the child lines for + // the value of the mother + // line (note that we have to + // multiply by two, since the + // shape functions scale with + // the inverse Jacobian). we + // thus throw away the + // information of one of the + // child lines, but there + // seems to be no other way + // than that... + // + // note: to make things + // consistent, and + // restriction independent of + // the order in which we + // travel across the cells of + // the coarse grid, we have + // to make sure that we take + // the same small line when + // visiting its two + // neighbors, to get the + // value for the mother + // line. we take the first + // line always, in the + // canonical direction of + // lines + for (unsigned int c=0; c::children_per_cell; ++c) + this->restriction[c].reinit (this->dofs_per_cell, + this->dofs_per_cell); + + this->restriction[0](0,0) = 2.; + this->restriction[1](1,1) = 2.; + this->restriction[3](2,2) = 2.; + this->restriction[0](3,3) = 2.; + + break; + }; + + + case 1: + { + for (unsigned int c=0; c::children_per_cell; ++c) + this->restriction[c].reinit (this->dofs_per_cell, + this->dofs_per_cell); + + // first set the corner + // nodes. note that they are + // non-additive + this->restriction[0](0,0) = 2.; + this->restriction[0](6,6) = 2.; + + this->restriction[1](1,1) = 2.; + this->restriction[1](2,2) = 2.; + + this->restriction[2](3,3) = 2.; + this->restriction[2](5,5) = 2.; + + this->restriction[3](4,4) = 2.; + this->restriction[3](7,7) = 2.; + + // then also set the bubble + // nodes. they _are_ + // additive. to understand + // what's going on, recall + // that the bubble shape + // functions have value -1 + // (!) at the center point, + // by construction of the + // polynomials, and that the + // corner nodes have values + // 1/2 there since they are + // just the linears, and not + // some interpolating + // polynomial + // + // (actually, the + // additive/non-additive + // business shouldn't make + // that much of a difference: + // node 4 on cell 0 and node + // 0 on cell 3 must have the + // same value, since normal + // components are + // continuous. so we could + // pick either and make these + // shape functions + // non-additive as well. we + // choose to take the mean + // value, which should be the + // same as either value, and + // make the shape function + // additive) + this->restriction[0](10,0) = 1.; + this->restriction[0](10,4) = -1.; + this->restriction[3](10,0) = -1.; + this->restriction[3](10,4) = 1.; + + this->restriction[1](11,1) = 1.; + this->restriction[1](11,5) = -1.; + this->restriction[2](11,1) = -1.; + this->restriction[2](11,5) = 1.; + + this->restriction[0](8,6) = 1.; + this->restriction[0](8,2) = -1.; + this->restriction[1](8,6) = -1.; + this->restriction[1](8,2) = 1.; + + this->restriction[3](8,7) = 1.; + this->restriction[3](8,3) = -1.; + this->restriction[2](8,7) = -1.; + this->restriction[2](8,3) = 1.; + + break; + }; + + // in case we don't have the + // matrices (yet), leave them + // empty. this does not + // prevent the use of this FE, + // but will prevent the use of + // these matrices + }; +} + +#endif + +#if deal_II_dimension == 3 + +template <> +void +FE_RaviartThomas<3>::initialize_restriction () +{ + Assert (false, ExcNotImplemented()); +} + +#endif + + +template +void +FE_RaviartThomas::initialize_unit_support_points () +{ + this->unit_support_points.resize (this->dofs_per_cell); + switch (dim) + { + case 2: + { + Assert (degree+1 == this->dofs_per_face, ExcInternalError()); + + // associate support points + // with mid-face points if a + // shape function has a + // non-zero normal component + // there, otherwise with the + // cell center. the reason + // for this non-unique + // support point is that we + // use hierarchical shape + // functions, rather than + // Lagrange functions, for + // which we get into the same + // trouble as in the + // FE_Q_Hierarchical element; + // see the respective + // function there + + // start with the face shape + // functions + for (unsigned int i=0; idofs_per_face; ++i) + this->unit_support_points[0*this->dofs_per_face+i] = Point(.5, .0); + for (unsigned int i=0; idofs_per_face; ++i) + this->unit_support_points[1*this->dofs_per_face+i] = Point(1., .5); + for (unsigned int i=0; idofs_per_face; ++i) + this->unit_support_points[2*this->dofs_per_face+i] = Point(.5, 1.); + for (unsigned int i=0; idofs_per_face; ++i) + this->unit_support_points[3*this->dofs_per_face+i] = Point(.0, .5); + + // associate the rest with + // the cell center + for (unsigned int i=4*this->dofs_per_face; idofs_per_cell; ++i) + this->unit_support_points[i] = Point(.5, .5); + + break; + } + + case 3: + { + // same as in 2d + Assert ((degree+1)*(degree+1) == this->dofs_per_face, ExcInternalError()); + + // start with the face shape + // functions + for (unsigned int i=0; idofs_per_face; ++i) + this->unit_support_points[0*this->dofs_per_face+i] = Point(.5, .0, .5); + for (unsigned int i=0; idofs_per_face; ++i) + this->unit_support_points[1*this->dofs_per_face+i] = Point(.5, 1., .5); + for (unsigned int i=0; idofs_per_face; ++i) + this->unit_support_points[2*this->dofs_per_face+i] = Point(.5, .5, 0.); + for (unsigned int i=0; idofs_per_face; ++i) + this->unit_support_points[3*this->dofs_per_face+i] = Point(1., .5, .5); + for (unsigned int i=0; idofs_per_face; ++i) + this->unit_support_points[4*this->dofs_per_face+i] = Point(.5, .5, 1.); + for (unsigned int i=0; idofs_per_face; ++i) + this->unit_support_points[5*this->dofs_per_face+i] = Point(.0, .5, .5); + + // associate the rest with + // the cell center + for (unsigned int i=6*this->dofs_per_face; idofs_per_cell; ++i) + this->unit_support_points[i] = Point(.5, .5, .5); + + break; + } + + default: + Assert (false, ExcNotImplemented()); + }; +} + + +#if deal_II_dimension == 1 + +template <> +void FE_RaviartThomas<1>::initialize_unit_face_support_points () +{ + // no faces in 1d, so nothing to do +} + +#endif + + +template +void FE_RaviartThomas::initialize_unit_face_support_points () +{ + this->unit_face_support_points.resize (this->dofs_per_face); + + // like with cell + // unit_support_points: + // associate all of the in + // the face mid-point, since + // there is no other useful + // way + for (unsigned int i=0; idofs_per_face; ++i) + this->unit_face_support_points[i] = (dim == 2 ? + Point(.5) : + Point(.5,.5)); +} + + +#if deal_II_dimension == 1 + +template <> +std::vector +FE_RaviartThomas<1>::get_dpo_vector (const unsigned int) +{ + Assert (false, ExcNotUsefulInThisDimension()); + return std::vector(); +} + +#endif + + +template +std::vector +FE_RaviartThomas::get_dpo_vector (const unsigned int degree) +{ + // the element is face-based (not + // to be confused with George + // W. Bush's Faith Based + // Initiative...), and we have + // (degree+1)^(dim-1) DoFs per face + unsigned int dofs_per_face = 1; + for (unsigned int d=0; d dpo(dim+1); + dpo[dim-1] = dofs_per_face; + dpo[dim] = interior_dofs; + + return dpo; +} + + + +#if deal_II_dimension == 1 + +template <> +std::vector +FE_RaviartThomas<1>::get_ria_vector (const unsigned int) +{ + Assert (false, ExcNotUsefulInThisDimension()); + return std::vector(); +} + +#endif + + +template +std::vector +FE_RaviartThomas::get_ria_vector (const unsigned int degree) +{ + unsigned int dofs_per_cell, dofs_per_face; + switch (dim) + { + case 2: + dofs_per_face = degree+1; + dofs_per_cell = 2*(degree+1)*(degree+2); + break; + case 3: + dofs_per_face = (degree+1)*(degree+1); + dofs_per_cell = 3*(degree+1)*(degree+1)*(degree+2); + break; + default: + Assert (false, ExcNotImplemented()); + } + Assert (FiniteElementData(get_dpo_vector(degree),dim).dofs_per_cell == + dofs_per_cell, + ExcInternalError()); + Assert (FiniteElementData(get_dpo_vector(degree),dim).dofs_per_face == + dofs_per_face, + ExcInternalError()); + + // all face dofs need to be + // non-additive, since they have + // continuity requirements. + // however, the interior dofs are + // made additive + std::vector ret_val(dofs_per_cell,false); + for (unsigned int i=GeometryInfo::faces_per_cell*dofs_per_face; + i < dofs_per_cell; ++i) + ret_val[i] = true; + + return ret_val; +} + + +#if deal_II_dimension == 1 + +template <> +std::vector > +FE_RaviartThomas<1>::create_polynomials (const unsigned int) +{ + Assert (false, ExcNotUsefulInThisDimension()); + return std::vector > (); +} + +#endif + + +#if deal_II_dimension == 2 + +template <> +std::vector > +FE_RaviartThomas<2>::create_polynomials (const unsigned int degree) +{ + const unsigned int dim = 2; + + // use the fact that the RT(k) + // spaces are spanned by the + // functions + // P_{k+1,k} \times P_{k,k+1}, + // see the book by Brezzi and + // Fortin + const std::vector > pols[2] + = { Polynomials::Hierarchical::generate_complete_basis (degree+1), + Polynomials::Hierarchical::generate_complete_basis (degree)}; + + // create spaces (k+1,k) and (k,k+1) + std::vector > > + pols_vector_1(dim), pols_vector_2(dim); + pols_vector_1[0] = pols[0]; + pols_vector_1[1] = pols[1]; + + pols_vector_2[0] = pols[1]; + pols_vector_2[1] = pols[0]; + + const AnisotropicPolynomials anisotropic[dim] + = { AnisotropicPolynomials (pols_vector_1), + AnisotropicPolynomials (pols_vector_2) }; + + // work around a stupid bug in + // gcc2.95 where the compiler + // complains about reaching the end + // of a non-void function when we + // simply return the following + // object unnamed, rather than + // first creating a named object + // and then returning it... + const std::vector > + ret_val (&anisotropic[0], &anisotropic[dim]); + return ret_val; +} + +#endif + + +#if deal_II_dimension == 3 + +template <> +std::vector > +FE_RaviartThomas<3>::create_polynomials (const unsigned int degree) +{ + const unsigned int dim = 3; + + // use the fact that the RT(k) + // spaces are spanned by the + // functions + // P_{k+1,k,k} \times P_{k,k+1,k} + // \times P_{k,k,k+1}, + // see the book by Brezzi and + // Fortin + const std::vector > pols[2] + = { Polynomials::Hierarchical::generate_complete_basis (degree+1), + Polynomials::Hierarchical::generate_complete_basis (degree)}; + + // create spaces (k+1,k,k), + // (k,k+1,k) and (k,k,k+1) + std::vector > > + pols_vector_1(dim), pols_vector_2(dim), pols_vector_3(dim); + pols_vector_1[0] = pols[0]; + pols_vector_1[1] = pols[1]; + pols_vector_1[2] = pols[1]; + + pols_vector_2[0] = pols[1]; + pols_vector_2[1] = pols[0]; + pols_vector_2[2] = pols[1]; + + pols_vector_3[0] = pols[1]; + pols_vector_3[1] = pols[1]; + pols_vector_3[2] = pols[0]; + + const AnisotropicPolynomials anisotropic[dim] + = { AnisotropicPolynomials (pols_vector_1), + AnisotropicPolynomials (pols_vector_2), + AnisotropicPolynomials (pols_vector_3) }; + + // work around a stupid bug in + // gcc2.95 where the compiler + // complains about reaching the end + // of a non-void function when we + // simply return the following + // object unnamed, rather than + // first creating a named object + // and then returning it... + const std::vector > + ret_val (&anisotropic[0], &anisotropic[dim]); + return ret_val; +} + +#endif + + + +#if deal_II_dimension == 1 + +template <> +std::vector > +FE_RaviartThomas<1>::compute_renumber (const unsigned int) +{ + Assert (false, ExcNotUsefulInThisDimension()); + return std::vector > (); +} + +#endif + + +#if deal_II_dimension == 2 + +template <> +std::vector > +FE_RaviartThomas<2>::compute_renumber (const unsigned int degree) +{ + const unsigned int dim = 2; + + std::vector > ret_val; + + // to explain the following: the + // first (degree+1) shape functions + // are on face 0, and point in + // y-direction, so are for the + // second vector component. then + // there are (degree+1) shape + // functions on face 1, which is + // for the x vector component, and + // so on. since the order of face + // degrees of freedom is arbitrary, + // we simply use the same order as + // that provided by the 1d + // polynomial class on which this + // element is based. after + // 4*(degree+1), the remaining + // shape functions are all bubbles, + // so we can number them in any way + // we want. we do so by first + // numbering the x-vectors, then + // the y-vectors + // + // now, we have to find a mapping + // from the above ordering to: + // first which vector component + // they belong to (easy), and + // second the index within this + // component as provided by the + // AnisotropicPolynomials class + // + // this is mostly a counting + // argument, tedious and error + // prone, and so boring to explain + // that we rather not try to do so + // here (it's simple, but boring, + // as said), aside from a few + // comments below + + // face 0 + for (unsigned int i=0; i test[dim] = { std::vector(ret_val.size()/dim, false), + std::vector(ret_val.size()/dim, false) }; + for (unsigned int i=0; i +std::vector > +FE_RaviartThomas<3>::compute_renumber (const unsigned int /*degree*/) +{ + Assert (false, ExcNotImplemented()); + return std::vector > (); +} + +#endif + + + + +template +UpdateFlags +FE_RaviartThomas::update_once (const UpdateFlags) const +{ + // even the values have to be + // computed on the real cell, so + // nothing can be done in advance + return update_default; +} + + + +template +UpdateFlags +FE_RaviartThomas::update_each (const UpdateFlags flags) const +{ + UpdateFlags out = update_default; + + if (flags & update_values) + out |= update_values | update_covariant_transformation; + 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_RaviartThomas::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; + + // initialize fields only if really + // necessary. otherwise, don't + // allocate memory + if (flags & update_values) + data->shape_values.reinit (this->dofs_per_cell, n_q_points); + + if (flags & update_gradients) + 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 values + // and gradients are only those on + // the unit cell, and need to be + // transformed when visiting an + // actual cell + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int q=0; qshape_values[i][q][c] + = shape_value_component(i,quadrature.point(q),c); + + if (flags & update_gradients) + for (unsigned int c=0; cshape_gradients[i][q][c] + = shape_grad_component(i,quadrature.point(q),c); + } + + return data; +} + + + + +//---------------------------------------------------------------------- +// Fill data of FEValues +//---------------------------------------------------------------------- + +template +void +FE_RaviartThomas::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); + + // get the flags indicating the + // fields that have to be filled + const UpdateFlags flags(fe_data.current_update_flags()); + + const unsigned int n_q_points = quadrature.n_quadrature_points; + + // fill shape function + // values. these are vector-valued, + // so we have to transform + // them. since the output format + // (in data.shape_values) is a + // sequence of doubles (one for + // each non-zero shape function + // value, and for each quadrature + // point, rather than a sequence of + // small vectors, we have to use a + // number of conversions + if (flags & update_values) + { + std::vector > shape_values (n_q_points); + + Assert (data.shape_values.n_rows() == this->dofs_per_cell * dim, + ExcInternalError()); + Assert (data.shape_values.n_cols() == n_q_points, + ExcInternalError()); + + for (unsigned int k=0; kdofs_per_cell; ++k) + { + // first transform shape + // values... + Assert (fe_data.shape_values[k].size() == n_q_points, + ExcInternalError()); + mapping.transform_covariant(&*shape_values.begin(), + &*shape_values.end(), + fe_data.shape_values[k].begin(), + mapping_data); + + // then copy over to target: + for (unsigned int q=0; q > shape_grads1 (n_q_points); + std::vector > shape_grads2 (n_q_points); + + Assert (data.shape_gradients.n_rows() == this->dofs_per_cell * dim, + ExcInternalError()); + Assert (data.shape_gradients.n_cols() == n_q_points, + ExcInternalError()); + + // loop over all shape + // functions, and treat the + // gradients of each shape + // function at all quadrature + // points + for (unsigned int k=0; kdofs_per_cell; ++k) + { + // treat the gradients of + // this particular shape + // function at all + // q-points. if Dv is the + // gradient of the shape + // function on the unit + // cell, then + // (J^-T)Dv(J^-1) is the + // value we want to have on + // the real cell. so, we + // will have to apply a + // covariant transformation + // to Dv twice. since the + // interface only allows + // multiplication with + // (J^-1) from the right, + // we have to trick a + // little in between + Assert (fe_data.shape_gradients[k].size() == n_q_points, + ExcInternalError()); + // do first transformation + mapping.transform_covariant(&*shape_grads1.begin(), + &*shape_grads1.end(), + fe_data.shape_gradients[k].begin(), + mapping_data); + // transpose matrix + for (unsigned int q=0; qcompute_2nd (mapping, cell, 0, mapping_data, fe_data, data); +} + + + +template +void +FE_RaviartThomas::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; + + // get the flags indicating the + // fields that have to be filled + const UpdateFlags flags(fe_data.current_update_flags()); + + const unsigned int n_q_points = quadrature.n_quadrature_points; + + // fill shape function + // values. these are vector-valued, + // so we have to transform + // them. since the output format + // (in data.shape_values) is a + // sequence of doubles (one for + // each non-zero shape function + // value, and for each quadrature + // point, rather than a sequence of + // small vectors, we have to use a + // number of conversions + if (flags & update_values) + { + Assert (fe_data.shape_values.n_cols() == + GeometryInfo::faces_per_cell * n_q_points, + ExcInternalError()); + + std::vector > shape_values (n_q_points); + + Assert (data.shape_values.n_rows() == this->dofs_per_cell * dim, + ExcInternalError()); + Assert (data.shape_values.n_cols() == n_q_points, + ExcInternalError()); + + for (unsigned int k=0; kdofs_per_cell; ++k) + { + // first transform shape + // values... + mapping.transform_covariant(&*shape_values.begin(), + &*shape_values.end(), + fe_data.shape_values[k].begin()+offset, + mapping_data); + + // then copy over to target: + for (unsigned int q=0; q::faces_per_cell * n_q_points, + ExcInternalError()); + + std::vector > shape_grads1 (n_q_points); + std::vector > shape_grads2 (n_q_points); + + Assert (data.shape_gradients.n_rows() == this->dofs_per_cell * dim, + ExcInternalError()); + Assert (data.shape_gradients.n_cols() == n_q_points, + ExcInternalError()); + + // loop over all shape + // functions, and treat the + // gradients of each shape + // function at all quadrature + // points + for (unsigned int k=0; kdofs_per_cell; ++k) + { + // treat the gradients of + // this particular shape + // function at all + // q-points. if Dv is the + // gradient of the shape + // function on the unit + // cell, then + // (J^-T)Dv(J^-1) is the + // value we want to have on + // the real cell. so, we + // will have to apply a + // covariant transformation + // to Dv twice. since the + // interface only allows + // multiplication with + // (J^-1) from the right, + // we have to trick a + // little in between + // + // do first transformation + mapping.transform_covariant(&*shape_grads1.begin(), + &*shape_grads1.end(), + fe_data.shape_gradients[k].begin()+offset, + mapping_data); + // transpose matrix + for (unsigned int q=0; qcompute_2nd (mapping, cell, offset, mapping_data, fe_data, data); +} + + + +template +void +FE_RaviartThomas::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 + // faces are stored contiguously) + const unsigned int offset = ((face * GeometryInfo::subfaces_per_face + subface) + * quadrature.n_quadrature_points); + + // get the flags indicating the + // fields that have to be filled + const UpdateFlags flags(fe_data.current_update_flags()); + + const unsigned int n_q_points = quadrature.n_quadrature_points; + + // fill shape function + // values. these are vector-valued, + // so we have to transform + // them. since the output format + // (in data.shape_values) is a + // sequence of doubles (one for + // each non-zero shape function + // value, and for each quadrature + // point, rather than a sequence of + // small vectors, we have to use a + // number of conversions + if (flags & update_values) + { + Assert (fe_data.shape_values.n_cols() == + GeometryInfo::faces_per_cell * + GeometryInfo::subfaces_per_face * + n_q_points, + ExcInternalError()); + + std::vector > shape_values (n_q_points); + + Assert (data.shape_values.n_rows() == this->dofs_per_cell * dim, + ExcInternalError()); + Assert (data.shape_values.n_cols() == n_q_points, + ExcInternalError()); + + for (unsigned int k=0; kdofs_per_cell; ++k) + { + // first transform shape + // values... + mapping.transform_covariant(&*shape_values.begin(), + &*shape_values.end(), + fe_data.shape_values[k].begin()+offset, + mapping_data); + + // then copy over to target: + for (unsigned int q=0; q::faces_per_cell * + GeometryInfo::subfaces_per_face * + n_q_points, + ExcInternalError()); + + std::vector > shape_grads1 (n_q_points); + std::vector > shape_grads2 (n_q_points); + + Assert (data.shape_gradients.n_rows() == this->dofs_per_cell * dim, + ExcInternalError()); + Assert (data.shape_gradients.n_cols() == n_q_points, + ExcInternalError()); + + // loop over all shape + // functions, and treat the + // gradients of each shape + // function at all quadrature + // points + for (unsigned int k=0; kdofs_per_cell; ++k) + { + // treat the gradients of + // this particular shape + // function at all + // q-points. if Dv is the + // gradient of the shape + // function on the unit + // cell, then + // (J^-T)Dv(J^-1) is the + // value we want to have on + // the real cell. so, we + // will have to apply a + // covariant transformation + // to Dv twice. since the + // interface only allows + // multiplication with + // (J^-1) from the right, + // we have to trick a + // little in between + // + // do first transformation + mapping.transform_covariant(&*shape_grads1.begin(), + &*shape_grads1.end(), + fe_data.shape_gradients[k].begin()+offset, + mapping_data); + // transpose matrix + for (unsigned int q=0; qcompute_2nd (mapping, cell, offset, mapping_data, fe_data, data); +} + + + +template +unsigned int +FE_RaviartThomas::n_base_elements () const +{ + return 1; +} + + + +template +const FiniteElement & +FE_RaviartThomas::base_element (const unsigned int index) const +{ + Assert (index==0, ExcIndexRange(index, 0, 1)); + return *this; +} + + + +template +unsigned int +FE_RaviartThomas::element_multiplicity (const unsigned int index) const +{ + Assert (index==0, ExcIndexRange(index, 0, 1)); + return 1; +} + + + +template +bool +FE_RaviartThomas::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)); + + switch (degree) + { + case 0: + { + switch (dim) + { + case 2: + { + // only on the one + // non-adjacent face + // are the values + // actually zero. list + // these in a table + const unsigned int + opposite_faces[GeometryInfo<2>::faces_per_cell] + = { 2, 3, 0, 1}; + + return (face_index != opposite_faces[shape_index]); + }; + + default: Assert (false, ExcNotImplemented()); + }; + }; + + default: // other degree + Assert (false, ExcNotImplemented()); + }; + + return true; +} + + + +template +unsigned int +FE_RaviartThomas::memory_consumption () const +{ + Assert (false, ExcNotImplemented ()); + return 0; +} + + + +template +unsigned int +FE_RaviartThomas::get_degree () const +{ + return degree; +} + + + +template class FE_RaviartThomas; diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index 53aca37dc6..911398c66a 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -21,6 +21,11 @@ #include #include +#ifdef HAVE_STD_STRINGSTREAM +# include +#else +# include +#endif /* ----------------------- FESystem::InternalData ------------------- */ @@ -222,6 +227,35 @@ FESystem::~FESystem () +template +std::string +FESystem::get_name () const +{ +#ifdef HAVE_STD_STRINGSTREAM + std::ostringstream namebuf; +#else + std::ostrstream namebuf; +#endif + + namebuf << "FESystem<" << dim << ">["; + for (unsigned int i=0; i FiniteElement* FESystem::clone() const @@ -407,6 +441,96 @@ FESystem::shape_grad_grad_component (const unsigned int i, +template +void +FESystem:: +get_interpolation_matrix (const FiniteElementBase &x_source_fe, + FullMatrix &interpolation_matrix) const +{ + Assert (interpolation_matrix.m() == this->dofs_per_cell, + ExcDimensionMismatch (interpolation_matrix.m(), + this->dofs_per_cell)); + Assert (interpolation_matrix.n() == x_source_fe.dofs_per_cell, + ExcDimensionMismatch (interpolation_matrix.m(), + x_source_fe.dofs_per_cell)); + + // there are certain conditions + // that the two elements have to + // satisfy so that this can work. + // + // condition 1: the other element + // must also be a system element + AssertThrow ((x_source_fe.get_name().find ("FESystem<") == 0) + || + (dynamic_cast*>(&x_source_fe) != 0), + typename FiniteElementBase:: + ExcInterpolationNotImplemented()); + + // ok, source is a system element, + // so we may be able to do the work + const FESystem &source_fe + = dynamic_cast&>(x_source_fe); + + // condition 2: same number of + // basis elements + AssertThrow (n_base_elements() == source_fe.n_base_elements(), + typename FiniteElementBase:: + ExcInterpolationNotImplemented()); + + // condition 3: same number of + // basis elements + for (unsigned int i=0; i:: + ExcInterpolationNotImplemented()); + + // ok, so let's try whether it + // works: + + // first let's see whether all the + // basis elements actually generate + // their interpolation matrices. if + // we get past the following loop, + // then apparently none of the + // called base elements threw an + // exception, so we're fine + // continuing and assembling the + // one big matrix from the small + // ones of the base elements + std::vector > base_matrices (n_base_elements()); + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int j=0; jsystem_to_base_table[i].first == + source_fe.system_to_base_table[j].first) + interpolation_matrix(i,j) + = (base_matrices[this->system_to_base_table[i].first.first] + (this->system_to_base_table[i].second, + source_fe.system_to_base_table[j].second)); +} + + + +//---------------------------------------------------------------------- +// Data field initialization +//---------------------------------------------------------------------- + + + template UpdateFlags FESystem::update_once (const UpdateFlags flags) const diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index 5746b9be20..815e2f5816 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -28,15 +28,50 @@ #include + +namespace +{ + // forwarder function for + // FE::get_interpolation_matrix. we + // will want to call that function + // for arbitrary FullMatrix + // types, but it only accepts + // double arguments. since it is a + // virtual function, this can also + // not be changed. so have a + // forwarder function that calls + // that function directly if + // T==double, and otherwise uses a + // temporary + template + void gim_forwarder (const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix) + { + fe2.get_interpolation_matrix (fe1, interpolation_matrix); + } + + + template + void gim_forwarder (const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix) + { + FullMatrix tmp (interpolation_matrix.m(), + interpolation_matrix.n()); + fe2.get_interpolation_matrix (fe1, tmp); + interpolation_matrix = tmp; + } +} + + template -void FETools::get_interpolation_matrix(const FiniteElement &fe1, - const FiniteElement &fe2, - FullMatrix &interpolation_matrix) +void FETools::get_interpolation_matrix (const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix) { Assert (fe1.n_components() == fe2.n_components(), ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); - Assert (fe1.is_primitive() == true, ExcFEMustBePrimitive()); - Assert (fe2.is_primitive() == true, ExcFEMustBePrimitive()); Assert(interpolation_matrix.m()==fe2.dofs_per_cell && interpolation_matrix.n()==fe1.dofs_per_cell, ExcMatrixDimensionMismatch(interpolation_matrix.m(), @@ -44,6 +79,30 @@ void FETools::get_interpolation_matrix(const FiniteElement &fe1, fe2.dofs_per_cell, fe1.dofs_per_cell)); + // first try the easy way: maybe + // the FE wants to implement things + // itself: + bool fe_implements_interpolation = true; + try + { + gim_forwarder (fe1, fe2, interpolation_matrix); + } + catch (typename FiniteElementBase::ExcInterpolationNotImplemented &) + { + // too bad.... + fe_implements_interpolation = false; + } + if (fe_implements_interpolation == true) + return; + + // uh, so this was not the + // case. hm. then do it the hard + // way. note that this will only + // work if the element is + // primitive, so check this first + Assert (fe1.is_primitive() == true, ExcFEMustBePrimitive()); + Assert (fe2.is_primitive() == true, ExcFEMustBePrimitive()); + // Initialize FEValues for fe1 at // the unit support points of the // fe2 element. @@ -76,8 +135,6 @@ void FETools::get_back_interpolation_matrix(const FiniteElement &fe1, { Assert (fe1.n_components() == fe2.n_components(), ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); - Assert (fe1.is_primitive() == true, ExcFEMustBePrimitive()); - Assert (fe2.is_primitive() == true, ExcFEMustBePrimitive()); Assert(interpolation_matrix.m()==fe1.dofs_per_cell && interpolation_matrix.n()==fe1.dofs_per_cell, ExcMatrixDimensionMismatch(interpolation_matrix.m(), @@ -104,8 +161,6 @@ void FETools::get_interpolation_difference_matrix(const FiniteElement &fe1, { Assert (fe1.n_components() == fe2.n_components(), ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); - Assert (fe1.is_primitive() == true, ExcFEMustBePrimitive()); - Assert (fe2.is_primitive() == true, ExcFEMustBePrimitive()); Assert(difference_matrix.m()==fe1.dofs_per_cell && difference_matrix.n()==fe1.dofs_per_cell, ExcMatrixDimensionMismatch(difference_matrix.m(), -- 2.39.5