#include <fe/fe_update_flags.h>
#include <fe/mapping.h>
+#include <string>
+
template<int dim> class FESystem;
const std::vector<bool> &restriction_is_additive_flags,
const std::vector<std::vector<bool> > &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
* 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<dim>::ExcInterpolationNotImplemented}.
+ */
+ virtual void
+ get_interpolation_matrix (const FiniteElementBase<dim> &source,
+ FullMatrix<double> &matrix) const;
+
/**
* Comparison operator. We also
* check for equality of the
<< "The interface matrix has a size of " << arg1
<< "x" << arg2
<< ", which is not reasonable in the present dimension.");
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInterpolationNotImplemented);
protected:
/**
*/
FE_DGP (const unsigned int k);
+ /**
+ * Return a string that uniquely
+ * identifies a finite
+ * element. This class returns
+ * @p{FE_DGP<dim>(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
*/
FE_DGPNonparametric (const unsigned int k);
+ /**
+ * Return a string that uniquely
+ * identifies a finite
+ * element. This class returns
+ * @p{FE_DGPNonparametric<dim>(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
*/
FE_DGQ (const unsigned int k);
+ /**
+ * Return a string that uniquely
+ * identifies a finite
+ * element. This class returns
+ * @p{FE_DGQ<dim>(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
* 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<dim>::ExcInterpolationNotImplemented}
+ * is thrown.
+ */
+ virtual void
+ get_interpolation_matrix (const FiniteElementBase<dim> &source,
+ FullMatrix<double> &matrix) const;
/**
* Number of base elements in a
* implemented there.
*
*
- * @author Wolfgang Bangerth, Anna Schneebeli, 2002
+ * @author Wolfgang Bangerth, Anna Schneebeli, 2002, 2003
*/
template <int dim>
class FE_Nedelec : public FiniteElement<dim>
*/
FE_Nedelec (const unsigned int p);
+ /**
+ * Return a string that uniquely
+ * identifies a finite
+ * element. This class returns
+ * @p{FE_Nedelec<dim>(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
*/
FE_Q (const unsigned int p);
+ /**
+ * Return a string that uniquely
+ * identifies a finite
+ * element. This class returns
+ * @p{FE_Q<dim>(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
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<dim>::ExcInterpolationNotImplemented}
+ * is thrown.
+ */
+ virtual void
+ get_interpolation_matrix (const FiniteElementBase<dim> &source,
+ FullMatrix<double> &matrix) const;
+
+ /**
* Number of base elements in a
* mixed discretization. Since
* this is a scalar element,
*/
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<dim>::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}
std::vector<unsigned int>
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[];
*/
FE_Q_Hierarchical (const unsigned int p);
+ /**
+ * Return a string that uniquely
+ * identifies a finite
+ * element. This class returns
+ * @p{FE_Q_Hierarchical<dim>(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
--- /dev/null
+//---------------------------------------------------------------
+// $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 <base/config.h>
+#include <base/tensor_product_polynomials.h>
+#include <grid/geometry_info.h>
+#include <fe/fe.h>
+
+template <int dim> 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 <int dim>
+class FE_RaviartThomas : public FiniteElement<dim>
+{
+ 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<dim>(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<dim> &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<dim> &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<dim> &p,
+ const unsigned int component) const;
+
+ /**
+ * Return the polynomial degree
+ * of this finite element,
+ * i.e. the value passed to the
+ * constructor.
+ */
+ unsigned int get_degree () const;
+
+ /**
+ * 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<dim>::ExcInterpolationNotImplemented}
+ * is thrown.
+ */
+ virtual void
+ get_interpolation_matrix (const FiniteElementBase<dim> &source,
+ FullMatrix<double> &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<dim> &
+ 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<dim> * clone() const;
+
+ /**
+ * Prepare internal data
+ * structures and fill in values
+ * independent of the cell.
+ */
+ virtual
+ typename Mapping<dim>::InternalDataBase *
+ get_data (const UpdateFlags,
+ const Mapping<dim>& mapping,
+ const Quadrature<dim>& quadrature) const ;
+
+ /**
+ * Implementation of the same
+ * function in
+ * @ref{FiniteElement}.
+ */
+ virtual void
+ fill_fe_values (const Mapping<dim> &mapping,
+ const typename DoFHandler<dim>::cell_iterator &cell,
+ const Quadrature<dim> &quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_internal,
+ typename Mapping<dim>::InternalDataBase &fe_internal,
+ FEValuesData<dim>& data) const;
+
+ /**
+ * Implementation of the same
+ * function in
+ * @ref{FiniteElement}.
+ */
+ virtual void
+ fill_fe_face_values (const Mapping<dim> &mapping,
+ const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const Quadrature<dim-1> &quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_internal,
+ typename Mapping<dim>::InternalDataBase &fe_internal,
+ FEValuesData<dim>& data) const;
+
+ /**
+ * Implementation of the same
+ * function in
+ * @ref{FiniteElement}.
+ */
+ virtual void
+ fill_fe_subface_values (const Mapping<dim> &mapping,
+ const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int sub_no,
+ const Quadrature<dim-1> &quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_internal,
+ typename Mapping<dim>::InternalDataBase &fe_internal,
+ FEValuesData<dim>& data) const;
+
+ private:
+ /**
+ * 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<AnisotropicPolynomials<dim> > 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<std::pair<unsigned int, unsigned int> > renumber;
+
+
+ /**
+ * Generate the polynomial spaces
+ * for the @ref{polynomials}
+ * member.
+ */
+ static std::vector<AnisotropicPolynomials<dim> >
+ 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<unsigned int>
+ 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<bool>
+ get_ria_vector (const unsigned int degree);
+
+ /**
+ * Compute the values of the
+ * @p{renumber} field.
+ */
+ static std::vector<std::pair<unsigned int, unsigned int> >
+ 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<dim>::InternalDataBase
+ {
+ public:
+ /**
+ * Array with shape function
+ * values in quadrature
+ * points. There is one row
+ * for each shape function,
+ * containing values for each
+ * quadrature point. 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 <int dim1> friend class FE_RaviartThomas;
+};
+
+
+/* -------------- declaration of explicit specializations ------------- */
+
+template <>
+void FE_RaviartThomas<1>::initialize_unit_face_support_points ();
+
+template <>
+std::vector<unsigned int> FE_RaviartThomas<1>::get_dpo_vector (const unsigned int);
+
+template <>
+std::vector<AnisotropicPolynomials<1> >
+FE_RaviartThomas<1>::create_polynomials (const unsigned int);
+
+template <>
+std::vector<AnisotropicPolynomials<2> >
+FE_RaviartThomas<2>::create_polynomials (const unsigned int);
+
+template <>
+std::vector<AnisotropicPolynomials<3> >
+FE_RaviartThomas<3>::create_polynomials (const unsigned int);
+
+template <>
+std::vector<std::pair<unsigned int, unsigned int> >
+FE_RaviartThomas<1>::compute_renumber (const unsigned int);
+
+template <>
+std::vector<std::pair<unsigned int, unsigned int> >
+FE_RaviartThomas<2>::compute_renumber (const unsigned int);
+
+template <>
+std::vector<std::pair<unsigned int, unsigned int> >
+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
* 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 <int dim>
class FESystem : public FiniteElement<dim>
*/
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<dim>[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
* the computation of these
* values to the base elements.
*/
- virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const;
+ virtual
+ Tensor<2,dim>
+ shape_grad_grad_component (const unsigned int i,
+ const Point<dim> &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<dim>::ExcInterpolationNotImplemented}
+ * is thrown.
+ */
+ virtual void
+ get_interpolation_matrix (const FiniteElementBase<dim> &source,
+ FullMatrix<double> &matrix) const;
- /**
+ /**
* Number of different base
* elements of this object.
*
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.");
};
static_cast<unsigned int>(-1));
}
+
+
+template <int dim>
+void
+FiniteElementBase<dim>::
+get_interpolation_matrix (const FiniteElementBase<dim> &,
+ FullMatrix<double> &) const
+{
+ // by default, no interpolation
+ // implemented. so throw exception,
+ // as documentation says
+ AssertThrow (false,
+ typename FiniteElementBase<dim>::
+ ExcInterpolationNotImplemented());
+}
+
template <int dim>
-bool FiniteElementBase<dim>::operator == (const FiniteElementBase<dim> &f) const
+bool
+FiniteElementBase<dim>::operator == (const FiniteElementBase<dim> &f) const
{
return ((static_cast<const FiniteElementData<dim>&>(*this) ==
static_cast<const FiniteElementData<dim>&>(f)) &&
#include <fe/fe_dgp.h>
#include <fe/fe_values.h>
+#ifdef HAVE_STD_STRINGSTREAM
+# include <sstream>
+#else
+# include <strstream>
+#endif
template <int dim>
+template <int dim>
+std::string
+FE_DGP<dim>::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 <int dim>
FiniteElement<dim> *
FE_DGP<dim>::clone() const
#include <fe/fe_dgp_nonparametric.h>
#include <fe/fe_values.h>
+#ifdef HAVE_STD_STRINGSTREAM
+# include <sstream>
+#else
+# include <strstream>
+#endif
template <int dim>
+template <int dim>
+std::string
+FE_DGPNonparametric<dim>::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 <int dim>
FiniteElement<dim> *
FE_DGPNonparametric<dim>::clone() const
#include <fe/fe_dgq.h>
#include <fe/fe_values.h>
+#ifdef HAVE_STD_STRINGSTREAM
+# include <sstream>
+#else
+# include <strstream>
+#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 <int dim> 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<unsigned int>(-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<unsigned int>(-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<N, ExcInternalError());
+ if (N==1)
+ return Point<1> (.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<N, ExcInternalError());
+
+ if (N==1)
+ return Point<2> (.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<N, ExcInternalError());
+ if (N==1)
+ return Point<3> (.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 <int dim>
+template <int dim>
+std::string
+FE_DGQ<dim>::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 <int dim>
FiniteElement<dim> *
FE_DGQ<dim>::clone() const
+template <int dim>
+void
+FE_DGQ<dim>::
+get_interpolation_matrix (const FiniteElementBase<dim> &x_source_fe,
+ FullMatrix<double> &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<const FE_DGQ<dim>*>(&x_source_fe) != 0),
+ typename FiniteElementBase<dim>::
+ ExcInterpolationNotImplemented());
+
+ // ok, source is a Q element, so
+ // we will be able to do the work
+ const FE_DGQ<dim> &source_fe
+ = dynamic_cast<const FE_DGQ<dim>&>(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<double> cell_interpolation (this->dofs_per_cell,
+ this->dofs_per_cell);
+ FullMatrix<double> source_interpolation (this->dofs_per_cell,
+ source_fe.dofs_per_cell);
+ FullMatrix<double> tmp (this->dofs_per_cell,
+ source_fe.dofs_per_cell);
+ for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+ {
+ // generate a point on this
+ // cell and evaluate the
+ // shape functions there
+ const Point<dim> p = generate_unit_point (j, this->dofs_per_cell,
+ int2type<dim>());
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ cell_interpolation(j,i)
+ = polynomial_space.compute_value (i, p);
+
+ for (unsigned int i=0; i<source_fe.dofs_per_cell; ++i)
+ source_interpolation(j,i)
+ = source_fe.polynomial_space.compute_value (i, p);
+ }
+
+ // then compute the
+ // interpolation matrix matrix
+ // for this coordinate
+ // direction
+ cell_interpolation.gauss_jordan ();
+ cell_interpolation.mmult (interpolation_matrix,
+ source_interpolation);
+
+ // cut off very small values
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
+ if (std::fabs(interpolation_matrix(i,j)) < 1e-15)
+ interpolation_matrix(i,j) = 0.;
+
+ // make sure that the row sum of
+ // each of the matrices is 1 at
+ // this point. this must be so
+ // since the shape functions sum up
+ // to 1
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ {
+ double sum = 0.;
+ for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
+ sum += interpolation_matrix(i,j);
+
+ Assert (std::fabs(sum-1) < 5e-14*std::max(degree,1U)*dim,
+ ExcInternalError());
+ }
+}
+
+
+
+
//----------------------------------------------------------------------
// Data field initialization
#include <fe/fe_nedelec.h>
#include <fe/fe_values.h>
+#ifdef HAVE_STD_STRINGSTREAM
+# include <sstream>
+#else
+# include <strstream>
+#endif
+
template <int dim>
FE_Nedelec<dim>::FE_Nedelec (const unsigned int degree)
+template <int dim>
+std::string
+FE_Nedelec<dim>::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 <int dim>
FiniteElement<dim> *
FE_Nedelec<dim>::clone() const
if (flags & update_values)
{
Assert (fe_data.shape_values.n_cols() ==
- GeometryInfo<dim>::faces_per_cell * n_q_points,
- ExcInternalError());
+ GeometryInfo<dim>::subfaces_per_face *
+ GeometryInfo<dim>::faces_per_cell *
+ n_q_points,
+ ExcInternalError());
std::vector<Tensor<1,dim> > shape_values (n_q_points);
if (flags & update_gradients)
{
Assert (fe_data.shape_gradients.n_cols() ==
- GeometryInfo<dim>::faces_per_cell * n_q_points,
+ GeometryInfo<dim>::faces_per_cell *
+ GeometryInfo<dim>::subfaces_per_face *
+ n_q_points,
ExcInternalError());
std::vector<Tensor<2,dim> > shape_grads1 (n_q_points);
#include <fe/fe_q.h>
#include <fe/fe_values.h>
+#ifdef HAVE_STD_STRINGSTREAM
+# include <sstream>
+#else
+# include <strstream>
+#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 <int dim> struct int2type {};
+
+ // given a permutation array,
+ // compute and return the inverse
+ // permutation
#ifdef DEAL_II_ANON_NAMESPACE_BUG
static
#endif
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<unsigned int>(-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<unsigned int>(-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<N, ExcInternalError());
+ 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<N, ExcInternalError());
+ 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<N, ExcInternalError());
+ 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);
+ }
+
+
+
+ // 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<N, ExcInternalError());
+ const double h = 1./(N+1);
+ return Point<1>((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<N, ExcInternalError());
+
+ const unsigned int N1d = int_sqrt(N+4);
+ const double h = 1./(N1d+1);
+
+ // i gives the index in the list
+ // of points excluding the four
+ // vertices. convert this into an
+ // index for all N1d**2 points
+ //
+ // we do so by
+ // - adding one if the point is
+ // beyond the lower left vertex
+ // (actually, all points are)
+
+ // - adding one if the point is
+ // beyond the lower right one
+ // - adding one if it is beyond
+ // the upper left one
+ // - not adding one for the upper
+ // right vertex, since no point
+ // can be beyond that one anyway
+ // :-)
+ const unsigned int true_i = (1
+ +
+ (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<unsigned int>(-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<unsigned int>(-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;
+ }
}
+template <int dim>
+std::string
+FE_Q<dim>::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 <int dim>
FiniteElement<dim> *
FE_Q<dim>::clone() const
}
+
+template <int dim>
+void
+FE_Q<dim>::
+get_interpolation_matrix (const FiniteElementBase<dim> &x_source_fe,
+ FullMatrix<double> &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<const FE_Q<dim>*>(&x_source_fe) != 0),
+ typename FiniteElementBase<dim>::
+ ExcInterpolationNotImplemented());
+
+ // ok, source is a Q element, so
+ // we will be able to do the work
+ const FE_Q<dim> &source_fe
+ = dynamic_cast<const FE_Q<dim>&>(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<double> cell_interpolation (this->dofs_per_cell,
+ this->dofs_per_cell);
+ FullMatrix<double> source_interpolation (this->dofs_per_cell,
+ source_fe.dofs_per_cell);
+ FullMatrix<double> tmp (this->dofs_per_cell,
+ source_fe.dofs_per_cell);
+ for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+ {
+ // generate a point on this
+ // cell and evaluate the
+ // shape functions there
+ const Point<dim> p = generate_unit_point (j, this->dofs_per_cell,
+ int2type<dim>());
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ cell_interpolation(renumber[j],renumber[i])
+ = polynomial_space.compute_value (i, p);
+
+ for (unsigned int i=0; i<source_fe.dofs_per_cell; ++i)
+ source_interpolation(renumber[j],source_fe.renumber[i])
+ = source_fe.polynomial_space.compute_value (i, p);
+ }
+
+ // then compute the
+ // interpolation matrix matrix
+ // for this coordinate
+ // direction
+ cell_interpolation.gauss_jordan ();
+ cell_interpolation.mmult (interpolation_matrix,
+ source_interpolation);
+
+ // cut off very small values
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
+ if (std::fabs(interpolation_matrix(i,j)) < 1e-15)
+ interpolation_matrix(i,j) = 0.;
+
+ // make sure that the row sum of
+ // each of the matrices is 1 at
+ // this point. this must be so
+ // since the shape functions sum up
+ // to 1
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ {
+ double sum = 0.;
+ for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
+ sum += interpolation_matrix(i,j);
+
+ Assert (std::fabs(sum-1) < 2e-14*degree*dim,
+ ExcInternalError());
+ }
+}
+
+
//----------------------------------------------------------------------
// Auxiliary functions
//----------------------------------------------------------------------
const unsigned int n = degree+1;
-
if (degree == 0)
{
Assert ((fe_data.dofs_per_vertex == 0) &&
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<renumber.size(), ExcInternalError());
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)
{
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<degree ;++jx)
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<degree; jy++)
}
-#if (deal_II_dimension == 1)
+#if deal_II_dimension == 1
template <>
std::vector<unsigned int>
#endif
+#if deal_II_dimension == 1
+template <>
+void
+FE_Q<1>::initialize_constraints ()
+{
+ // no constraints in 1d
+}
-template <int dim>
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
void
-FE_Q<dim>::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<dim-1>
+ 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<double> face_interpolation (n_small_functions, this->dofs_per_face);
+ FullMatrix<double> subface_interpolation (n_small_functions, n_small_functions);
+
+ const std::vector<unsigned int>
+ face_renumber_inverse (invert_numbering(face_renumber));
+
+ for (unsigned int i=0; i<n_small_functions; ++i)
{
- this->interface_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<dim-1> p_face = generate_face_unit_point (i, n_small_functions,
+ int2type<dim-1>());
+
+ // 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; j<this->dofs_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<n_small_functions; ++j)
+ {
+ // first thing is to check
+ // which face the present
+ // point is on,
+ // i.e. whether it is left
+ // or right of the middle
+ // vertex in 2d, or
+ // something more complex
+ // in 3d (note that we
+ // might actually be
+ // sitting on top of the
+ // center vertex, or on on
+ // interface between
+ // children, but that
+ // doesn't really bother
+ // us: the shape functions
+ // associated with that
+ // have the same value
+ // whether we consider the
+ // left or the right
+ // subface, and all other
+ // shape functions should
+ // be zero there as well,
+ // so it doesn't really
+ // matter whether we
+ // account for this fact or
+ // not...)
+ const unsigned int subface
+ = GeometryInfo<dim-1>::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<dim-1> p_subface
+ = GeometryInfo<dim-1>::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<double> x;
+ x.TableBase<2,double>::reinit (this->interface_constraints_size());
+ x.fill (Matrices::constraint_matrices[degree-1]);
+
+ for (unsigned int i=0; i<x.m(); ++i)
+ for (unsigned int j=0; j<x.n(); ++j)
+ Assert (std::fabs (x(i,j) - this->interface_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 <int dim>
void
FE_Q<dim>::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<GeometryInfo<dim>::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; row<this->dofs_per_cell; ++row)
- {
- double sum = 0;
- for (unsigned int col=0; col<this->dofs_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<double> cell_interpolation (this->dofs_per_cell,
+ this->dofs_per_cell);
+ FullMatrix<double> subcell_interpolation (this->dofs_per_cell,
+ this->dofs_per_cell);
+ for (unsigned int child=0; child<GeometryInfo<dim>::children_per_cell; ++child)
+ this->prolongation[child].reinit (this->dofs_per_cell,
+ this->dofs_per_cell);
+ for (unsigned int child=0; child<GeometryInfo<dim>::children_per_cell; ++child)
+ {
+ for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+ {
+ // generate a point on
+ // the child cell and
+ // evaluate the shape
+ // functions there
+ const Point<dim> p_subcell = generate_unit_point (j, this->dofs_per_cell,
+ int2type<dim>());
+ const Point<dim> p_cell =
+ GeometryInfo<dim>::child_to_cell_coordinates (p_subcell, child);
+
+ for (unsigned int i=0; i<this->dofs_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; i<this->dofs_per_cell; ++i)
+ for (unsigned int j=0; j<this->dofs_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; row<this->dofs_per_cell; ++row)
+ {
+ double sum = 0;
+ for (unsigned int col=0; col<this->dofs_per_cell; ++col)
+ sum += this->prolongation[child](row,col);
+ Assert (std::fabs(sum-1.) < 2e-14*degree*dim,
+ ExcInternalError());
+ };
+ }
}
void
FE_Q<dim>::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<GeometryInfo<dim>::children_per_cell; ++c)
this->restriction[c].reinit (this->dofs_per_cell, this->dofs_per_cell);
- switch (dim)
+ for (unsigned int i=0; i<this->dofs_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<GeometryInfo<dim>::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<GeometryInfo<dim>::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<dim> p_cell = generate_unit_point (i, this->dofs_per_cell,
+ int2type<dim>());
+ unsigned int mother_dof = 0;
+ for (; mother_dof<this->dofs_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; j<this->dofs_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<GeometryInfo<dim>::children_per_cell;
+ ++child)
+ {
+ // first initialize this
+ // column of the matrix
+ for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+ this->restriction[child](mother_dof, j) = 0.;
+
+ // then check whether this
+ // interpolation point is
+ // inside this child cell
+ const Point<dim> p_subcell
+ = GeometryInfo<dim>::cell_to_child_coordinates (p_cell, child);
+ if (GeometryInfo<dim>::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<GeometryInfo<dim>::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_dof<this->dofs_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; j<this->dofs_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.;
+ }
+ }
}
}
#include <fe/fe_q.h>
-// 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 <>
#include <fe/fe_q.h>
-// 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;
#include <fe/fe_q.h>
-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
#include <fe/fe_values.h>
#include <cmath>
+#ifdef HAVE_STD_STRINGSTREAM
+# include <sstream>
+#else
+# include <strstream>
+#endif
+
namespace
{
+template <int dim>
+std::string
+FE_Q_Hierarchical<dim>::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 <int dim>
FiniteElement<dim> *
FE_Q_Hierarchical<dim>::clone() const
--- /dev/null
+//----------------------------------------------------------------
+// $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 <base/quadrature.h>
+#include <base/table.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe.h>
+#include <fe/mapping.h>
+#include <fe/fe_raviart_thomas.h>
+#include <fe/fe_values.h>
+
+#ifdef HAVE_STD_STRINGSTREAM
+# include <sstream>
+#else
+# include <strstream>
+#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 <int dim> 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<unsigned int>(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<unsigned int>(std::pow(1.*N, 1./3.));
+ const unsigned int N2 = static_cast<unsigned int>(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 <int dim>
+FE_RaviartThomas<dim>::FE_RaviartThomas (const unsigned int degree)
+ :
+ FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree),
+ dim),
+ get_ria_vector (degree),
+ std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),dim).dofs_per_cell,
+ std::vector<bool>(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<std::pair<unsigned,unsigned> > tmp1, tmp2;
+ this->system_to_component_table.swap (tmp1);
+ this->face_system_to_component_table.swap (tmp2);
+}
+
+
+
+template <int dim>
+std::string
+FE_RaviartThomas<dim>::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 <int dim>
+FiniteElement<dim> *
+FE_RaviartThomas<dim>::clone() const
+{
+ return new FE_RaviartThomas<dim>(degree);
+}
+
+
+template <int dim>
+double
+FE_RaviartThomas<dim>::shape_value_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
+ Assert (component < 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 <int dim>
+Tensor<1,dim>
+FE_RaviartThomas<dim>::shape_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
+ Assert (component < 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 <int dim>
+Tensor<2,dim>
+FE_RaviartThomas<dim>::shape_grad_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
+ Assert (component < 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<double> &) const
+{
+ Assert (false, ExcNotUsefulInThisDimension());
+}
+
+#endif
+
+
+template <int dim>
+void
+FE_RaviartThomas<dim>::
+get_interpolation_matrix (const FiniteElementBase<dim> &x_source_fe,
+ FullMatrix<double> &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<const FE_RaviartThomas<dim>*>(&x_source_fe) != 0),
+ typename FiniteElementBase<dim>::
+ ExcInterpolationNotImplemented());
+
+ // ok, source is a RT element, so
+ // we will be able to do the work
+ const FE_RaviartThomas<dim> &source_fe
+ = dynamic_cast<const FE_RaviartThomas<dim>&>(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<dim; ++d)
+ Assert (polynomials[d].n() == dofs_per_coordinate, ExcInternalError());
+
+ const unsigned int source_dofs_per_coordinate = source_fe.dofs_per_cell/dim;
+ Assert (source_dofs_per_coordinate*dim == source_fe.dofs_per_cell,
+ ExcInternalError());
+ for (unsigned int d=0; d<dim; ++d)
+ Assert (source_fe.polynomials[d].n() == source_dofs_per_coordinate, ExcInternalError());
+
+ FullMatrix<double> cell_interpolation (dofs_per_coordinate,
+ dofs_per_coordinate);
+ FullMatrix<double> source_interpolation (dofs_per_coordinate,
+ source_dofs_per_coordinate);
+ FullMatrix<double> tmp (dofs_per_coordinate,
+ source_dofs_per_coordinate);
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ for (unsigned int j=0; j<dofs_per_coordinate; ++j)
+ {
+ // generate a point on this
+ // cell and evaluate the
+ // shape functions there
+ //
+ // see the comment for that
+ // function to see why the
+ // third parameter is
+ // necessary
+ const Point<dim> p = generate_unit_point (j, dofs_per_coordinate,
+ d, int2type<dim>());
+ for (unsigned int i=0; i<dofs_per_coordinate; ++i)
+ cell_interpolation(j,i) = polynomials[d].compute_value (i, p);
+
+ for (unsigned int i=0; i<source_dofs_per_coordinate; ++i)
+ source_interpolation(j,i) = source_fe.polynomials[d].compute_value (i, p);
+ }
+
+ // then compute the
+ // interpolation matrix matrix
+ // for this coordinate
+ // direction
+ cell_interpolation.gauss_jordan ();
+ cell_interpolation.mmult (tmp, source_interpolation);
+
+ // finally transfer the
+ // results for this
+ // coordinate into the matrix
+ // corresponding to the
+ // entire space on this
+ // cell. cut off very small
+ // values here
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ if (renumber[i].first == d)
+ for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
+ if (source_fe.renumber[j].first == d)
+ if (std::fabs(tmp(renumber[i].second,
+ source_fe.renumber[j].second)) > 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; i<this->dofs_per_cell; ++i)
+ for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
+ Assert ((interpolation_matrix(i,j) == 0.) ||
+ (renumber[i].first == source_fe.renumber[j].first),
+ ExcInternalError());
+}
+
+
+
+//----------------------------------------------------------------------
+// Auxiliary and internal functions
+//----------------------------------------------------------------------
+
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+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<Polynomials::Polynomial<double> >
+ face_polynomials (Polynomials::Hierarchical::
+ generate_complete_basis (degree));
+ Assert (face_polynomials.size() == this->dofs_per_face, ExcInternalError());
+
+ FullMatrix<double> face_interpolation (2*this->dofs_per_face, this->dofs_per_face);
+ FullMatrix<double> 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<GeometryInfo<dim>::subfaces_per_face; ++subface)
+ for (unsigned int i=0; i<this->dofs_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; j<this->dofs_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; i<this->interface_constraints.m(); ++i)
+ for (unsigned int j=0; j<this->interface_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 <int dim>
+void
+FE_RaviartThomas<dim>::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<dim; ++d)
+ Assert (polynomials[d].n() == dofs_per_coordinate, ExcInternalError());
+
+ FullMatrix<double> cell_interpolation (dofs_per_coordinate,
+ dofs_per_coordinate);
+ FullMatrix<double> subcell_interpolation (dofs_per_coordinate,
+ dofs_per_coordinate);
+ FullMatrix<double> tmp (dofs_per_coordinate,
+ dofs_per_coordinate);
+ for (unsigned int child=0; child<GeometryInfo<dim>::children_per_cell; ++child)
+ this->prolongation[child].reinit (this->dofs_per_cell,
+ this->dofs_per_cell);
+ for (unsigned int child=0; child<GeometryInfo<dim>::children_per_cell; ++child)
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ for (unsigned int j=0; j<dofs_per_coordinate; ++j)
+ {
+ // generate a point on
+ // the child cell and
+ // evaluate the shape
+ // functions there
+ //
+ // see the comment for
+ // that function to see
+ // why the third
+ // parameter is necessary
+ const Point<dim> p_subcell = generate_unit_point (j, dofs_per_coordinate,
+ d, int2type<dim>());
+ const Point<dim> p_cell =
+ GeometryInfo<dim>::child_to_cell_coordinates (p_subcell, child);
+
+ for (unsigned int i=0; i<dofs_per_coordinate; ++i)
+ {
+ cell_interpolation(j,i) = polynomials[d].compute_value (i, p_cell);
+ subcell_interpolation(j,i) = polynomials[d].compute_value (i, p_subcell);
+ }
+ }
+
+ // then compute the embedding
+ // matrix for this child and
+ // this coordinate direction
+ subcell_interpolation.gauss_jordan ();
+ subcell_interpolation.mmult (tmp, cell_interpolation);
+
+ // finally transfer the
+ // results for this
+ // coordinate into the matrix
+ // corresponding to the
+ // entire space on this
+ // cell. cut off very small
+ // values here
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ if (renumber[i].first == d)
+ for (unsigned int j=0; j<this->dofs_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<GeometryInfo<dim>::children_per_cell; ++child)
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ for (unsigned int j=0; j<this->dofs_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<GeometryInfo<dim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::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 <int dim>
+void
+FE_RaviartThomas<dim>::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; i<this->dofs_per_face; ++i)
+ this->unit_support_points[0*this->dofs_per_face+i] = Point<dim>(.5, .0);
+ for (unsigned int i=0; i<this->dofs_per_face; ++i)
+ this->unit_support_points[1*this->dofs_per_face+i] = Point<dim>(1., .5);
+ for (unsigned int i=0; i<this->dofs_per_face; ++i)
+ this->unit_support_points[2*this->dofs_per_face+i] = Point<dim>(.5, 1.);
+ for (unsigned int i=0; i<this->dofs_per_face; ++i)
+ this->unit_support_points[3*this->dofs_per_face+i] = Point<dim>(.0, .5);
+
+ // associate the rest with
+ // the cell center
+ for (unsigned int i=4*this->dofs_per_face; i<this->dofs_per_cell; ++i)
+ this->unit_support_points[i] = Point<dim>(.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; i<this->dofs_per_face; ++i)
+ this->unit_support_points[0*this->dofs_per_face+i] = Point<dim>(.5, .0, .5);
+ for (unsigned int i=0; i<this->dofs_per_face; ++i)
+ this->unit_support_points[1*this->dofs_per_face+i] = Point<dim>(.5, 1., .5);
+ for (unsigned int i=0; i<this->dofs_per_face; ++i)
+ this->unit_support_points[2*this->dofs_per_face+i] = Point<dim>(.5, .5, 0.);
+ for (unsigned int i=0; i<this->dofs_per_face; ++i)
+ this->unit_support_points[3*this->dofs_per_face+i] = Point<dim>(1., .5, .5);
+ for (unsigned int i=0; i<this->dofs_per_face; ++i)
+ this->unit_support_points[4*this->dofs_per_face+i] = Point<dim>(.5, .5, 1.);
+ for (unsigned int i=0; i<this->dofs_per_face; ++i)
+ this->unit_support_points[5*this->dofs_per_face+i] = Point<dim>(.0, .5, .5);
+
+ // associate the rest with
+ // the cell center
+ for (unsigned int i=6*this->dofs_per_face; i<this->dofs_per_cell; ++i)
+ this->unit_support_points[i] = Point<dim>(.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 <int dim>
+void FE_RaviartThomas<dim>::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; i<this->dofs_per_face; ++i)
+ this->unit_face_support_points[i] = (dim == 2 ?
+ Point<dim-1>(.5) :
+ Point<dim-1>(.5,.5));
+}
+
+
+#if deal_II_dimension == 1
+
+template <>
+std::vector<unsigned int>
+FE_RaviartThomas<1>::get_dpo_vector (const unsigned int)
+{
+ Assert (false, ExcNotUsefulInThisDimension());
+ return std::vector<unsigned int>();
+}
+
+#endif
+
+
+template <int dim>
+std::vector<unsigned int>
+FE_RaviartThomas<dim>::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<dim-1; ++d)
+ dofs_per_face *= degree+1;
+
+ // and then there are interior dofs
+ const unsigned int
+ interior_dofs = dim*degree*dofs_per_face;
+
+ std::vector<unsigned int> dpo(dim+1);
+ dpo[dim-1] = dofs_per_face;
+ dpo[dim] = interior_dofs;
+
+ return dpo;
+}
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+std::vector<bool>
+FE_RaviartThomas<1>::get_ria_vector (const unsigned int)
+{
+ Assert (false, ExcNotUsefulInThisDimension());
+ return std::vector<bool>();
+}
+
+#endif
+
+
+template <int dim>
+std::vector<bool>
+FE_RaviartThomas<dim>::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<dim>(get_dpo_vector(degree),dim).dofs_per_cell ==
+ dofs_per_cell,
+ ExcInternalError());
+ Assert (FiniteElementData<dim>(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<bool> ret_val(dofs_per_cell,false);
+ for (unsigned int i=GeometryInfo<dim>::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<AnisotropicPolynomials<1> >
+FE_RaviartThomas<1>::create_polynomials (const unsigned int)
+{
+ Assert (false, ExcNotUsefulInThisDimension());
+ return std::vector<AnisotropicPolynomials<1> > ();
+}
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
+std::vector<AnisotropicPolynomials<2> >
+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<Polynomials::Polynomial<double> > 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<std::vector<Polynomials::Polynomial<double> > >
+ 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<dim> anisotropic[dim]
+ = { AnisotropicPolynomials<dim> (pols_vector_1),
+ AnisotropicPolynomials<dim> (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<AnisotropicPolynomials<dim> >
+ ret_val (&anisotropic[0], &anisotropic[dim]);
+ return ret_val;
+}
+
+#endif
+
+
+#if deal_II_dimension == 3
+
+template <>
+std::vector<AnisotropicPolynomials<3> >
+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<Polynomials::Polynomial<double> > 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<std::vector<Polynomials::Polynomial<double> > >
+ 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<dim> anisotropic[dim]
+ = { AnisotropicPolynomials<dim> (pols_vector_1),
+ AnisotropicPolynomials<dim> (pols_vector_2),
+ AnisotropicPolynomials<dim> (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<AnisotropicPolynomials<dim> >
+ ret_val (&anisotropic[0], &anisotropic[dim]);
+ return ret_val;
+}
+
+#endif
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+std::vector<std::pair<unsigned int, unsigned int> >
+FE_RaviartThomas<1>::compute_renumber (const unsigned int)
+{
+ Assert (false, ExcNotUsefulInThisDimension());
+ return std::vector<std::pair<unsigned int, unsigned int> > ();
+}
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
+std::vector<std::pair<unsigned int, unsigned int> >
+FE_RaviartThomas<2>::compute_renumber (const unsigned int degree)
+{
+ const unsigned int dim = 2;
+
+ std::vector<std::pair<unsigned int, unsigned int> > 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<degree+1; ++i)
+ ret_val.push_back (std::make_pair (1, i));
+
+ // face 1
+ for (unsigned int i=0; i<degree+1; ++i)
+ ret_val.push_back (std::make_pair (0, (degree+2)*i+1));
+
+ // face 2
+ for (unsigned int i=0; i<degree+1; ++i)
+ ret_val.push_back (std::make_pair (1, (degree+1)+i));
+
+ // face 3
+ for (unsigned int i=0; i<degree+1; ++i)
+ ret_val.push_back (std::make_pair (0, (degree+2)*i));
+
+ // then go on with interior bubble
+ // functions, first for the
+ // x-direction, then for the
+ // y-direction
+ for (unsigned int x=0; x<degree; ++x)
+ for (unsigned int y=0; y<degree+1; ++y)
+ {
+ const unsigned int index_in_component = (x+2) + y*(degree+2);
+ Assert (index_in_component < (degree+1)*(degree+2),
+ ExcInternalError());
+ ret_val.push_back (std::make_pair(0, index_in_component));
+ }
+ for (unsigned int x=0; x<degree+1; ++x)
+ for (unsigned int y=0; y<degree; ++y)
+ {
+ const unsigned int index_in_component = 2*(degree+1) + y + x*degree;
+ Assert (index_in_component < (degree+1)*(degree+2),
+ ExcInternalError());
+ ret_val.push_back (std::make_pair(1, index_in_component));
+ }
+
+#ifdef DEBUG
+ // make sure we have actually used
+ // up all elements of the tensor
+ // product polynomial
+ Assert (ret_val.size() == 2*(degree+1)*(degree+2),
+ ExcInternalError());
+ std::vector<bool> test[dim] = { std::vector<bool>(ret_val.size()/dim, false),
+ std::vector<bool>(ret_val.size()/dim, false) };
+ for (unsigned int i=0; i<ret_val.size(); ++i)
+ {
+ Assert (ret_val[i].first < dim, ExcInternalError());
+ Assert (ret_val[i].second < test[dim].size(), ExcInternalError());
+ Assert (test[ret_val[i].first][ret_val[i].second] == false,
+ ExcInternalError());
+
+ test[ret_val[i].first][ret_val[i].second] = true;
+ }
+ for (unsigned int d=0; d<dim; ++d)
+ Assert (std::find (test[d].begin(), test[d].end(), false) == test[d].end(),
+ ExcInternalError());
+#endif
+
+ return ret_val;
+}
+
+#endif
+
+
+#if deal_II_dimension == 3
+
+template <>
+std::vector<std::pair<unsigned int, unsigned int> >
+FE_RaviartThomas<3>::compute_renumber (const unsigned int /*degree*/)
+{
+ Assert (false, ExcNotImplemented());
+ return std::vector<std::pair<unsigned int, unsigned int> > ();
+}
+
+#endif
+
+
+
+
+template <int dim>
+UpdateFlags
+FE_RaviartThomas<dim>::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 <int dim>
+UpdateFlags
+FE_RaviartThomas<dim>::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 <int dim>
+typename Mapping<dim>::InternalDataBase *
+FE_RaviartThomas<dim>::get_data (const UpdateFlags update_flags,
+ const Mapping<dim> &mapping,
+ const Quadrature<dim> &quadrature) const
+{
+ // generate a new data object and
+ // initialize some fields
+ InternalData* data = new InternalData;
+
+ // check what needs to be
+ // initialized only once and what
+ // on every cell/face/subface we
+ // visit
+ data->update_once = update_once(update_flags);
+ data->update_each = update_each(update_flags);
+ data->update_flags = data->update_once | data->update_each;
+
+ const UpdateFlags flags(data->update_flags);
+ const unsigned int n_q_points = quadrature.n_quadrature_points;
+
+ // 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; i<this->dofs_per_cell; ++i)
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ if (flags & update_values)
+ for (unsigned int c=0; c<dim; ++c)
+ data->shape_values[i][q][c]
+ = shape_value_component(i,quadrature.point(q),c);
+
+ if (flags & update_gradients)
+ for (unsigned int c=0; c<dim; ++c)
+ data->shape_gradients[i][q][c]
+ = shape_grad_component(i,quadrature.point(q),c);
+ }
+
+ return data;
+}
+
+
+
+
+//----------------------------------------------------------------------
+// Fill data of FEValues
+//----------------------------------------------------------------------
+
+template <int dim>
+void
+FE_RaviartThomas<dim>::fill_fe_values (const Mapping<dim> &mapping,
+ const typename DoFHandler<dim>::cell_iterator &cell,
+ const Quadrature<dim> &quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_data,
+ typename Mapping<dim>::InternalDataBase &fedata,
+ FEValuesData<dim> &data) const
+{
+ // convert data object to internal
+ // data for this class. fails with
+ // an exception if that is not
+ // possible
+ InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
+
+ // 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<Tensor<1,dim> > 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; k<this->dofs_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<n_q_points; ++q)
+ for (unsigned int d=0; d<dim; ++d)
+ data.shape_values[k*dim+d][q] = shape_values[q][d];
+ };
+ };
+
+
+ if (flags & update_gradients)
+ {
+ std::vector<Tensor<2,dim> > shape_grads1 (n_q_points);
+ std::vector<Tensor<2,dim> > 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; k<this->dofs_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; q<n_q_points; ++q)
+ shape_grads2[q] = transpose(shape_grads1[q]);
+ // do second transformation
+ mapping.transform_covariant(&*shape_grads1.begin(),
+ &*shape_grads1.end(),
+ &*shape_grads2.begin(),
+ mapping_data);
+ // transpose back
+ for (unsigned int q=0; q<n_q_points; ++q)
+ shape_grads2[q] = transpose(shape_grads1[q]);
+
+ // then copy over to target:
+ for (unsigned int q=0; q<n_q_points; ++q)
+ for (unsigned int d=0; d<dim; ++d)
+ data.shape_gradients[k*dim+d][q] = shape_grads2[q][d];
+ };
+ }
+
+ if (flags & update_second_derivatives)
+ this->compute_2nd (mapping, cell, 0, mapping_data, fe_data, data);
+}
+
+
+
+template <int dim>
+void
+FE_RaviartThomas<dim>::fill_fe_face_values (const Mapping<dim> &mapping,
+ const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face,
+ const Quadrature<dim-1> &quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_data,
+ typename Mapping<dim>::InternalDataBase &fedata,
+ FEValuesData<dim> &data) const
+{
+ // convert data object to internal
+ // data for this class. fails with
+ // an exception if that is not
+ // possible
+ InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
+
+ // offset determines which data set
+ // to take (all data sets for all
+ // faces are stored contiguously)
+ const unsigned int offset = face * quadrature.n_quadrature_points;
+
+ // 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<dim>::faces_per_cell * n_q_points,
+ ExcInternalError());
+
+ std::vector<Tensor<1,dim> > 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; k<this->dofs_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<n_q_points; ++q)
+ for (unsigned int d=0; d<dim; ++d)
+ data.shape_values[k*dim+d][q] = shape_values[q][d];
+ };
+ };
+
+
+ if (flags & update_gradients)
+ {
+ Assert (fe_data.shape_gradients.n_cols() ==
+ GeometryInfo<dim>::faces_per_cell * n_q_points,
+ ExcInternalError());
+
+ std::vector<Tensor<2,dim> > shape_grads1 (n_q_points);
+ std::vector<Tensor<2,dim> > 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; k<this->dofs_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; q<n_q_points; ++q)
+ shape_grads2[q] = transpose(shape_grads1[q]);
+ // do second transformation
+ mapping.transform_covariant(&*shape_grads1.begin(),
+ &*shape_grads1.end(),
+ &*shape_grads2.begin(),
+ mapping_data);
+ // transpose back
+ for (unsigned int q=0; q<n_q_points; ++q)
+ shape_grads2[q] = transpose(shape_grads1[q]);
+
+ // then copy over to target:
+ for (unsigned int q=0; q<n_q_points; ++q)
+ for (unsigned int d=0; d<dim; ++d)
+ data.shape_gradients[k*dim+d][q] = shape_grads2[q][d];
+ };
+ }
+
+ if (flags & update_second_derivatives)
+ this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+}
+
+
+
+template <int dim>
+void
+FE_RaviartThomas<dim>::fill_fe_subface_values (const Mapping<dim> &mapping,
+ const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face,
+ const unsigned int subface,
+ const Quadrature<dim-1> &quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_data,
+ typename Mapping<dim>::InternalDataBase &fedata,
+ FEValuesData<dim> &data) const
+{
+ // convert data object to internal
+ // data for this class. fails with
+ // an exception if that is not
+ // possible
+ InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
+
+ // offset determines which data set
+ // to take (all data sets for all
+ // faces are stored contiguously)
+ const unsigned int offset = ((face * GeometryInfo<dim>::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<dim>::faces_per_cell *
+ GeometryInfo<dim>::subfaces_per_face *
+ n_q_points,
+ ExcInternalError());
+
+ std::vector<Tensor<1,dim> > 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; k<this->dofs_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<n_q_points; ++q)
+ for (unsigned int d=0; d<dim; ++d)
+ data.shape_values[k*dim+d][q] = shape_values[q][d];
+ };
+ };
+
+
+ if (flags & update_gradients)
+ {
+ Assert (fe_data.shape_gradients.n_cols() ==
+ GeometryInfo<dim>::faces_per_cell *
+ GeometryInfo<dim>::subfaces_per_face *
+ n_q_points,
+ ExcInternalError());
+
+ std::vector<Tensor<2,dim> > shape_grads1 (n_q_points);
+ std::vector<Tensor<2,dim> > 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; k<this->dofs_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; q<n_q_points; ++q)
+ shape_grads2[q] = transpose(shape_grads1[q]);
+ // do second transformation
+ mapping.transform_covariant(&*shape_grads1.begin(),
+ &*shape_grads1.end(),
+ &*shape_grads2.begin(),
+ mapping_data);
+ // transpose back
+ for (unsigned int q=0; q<n_q_points; ++q)
+ shape_grads2[q] = transpose(shape_grads1[q]);
+
+ // then copy over to target:
+ for (unsigned int q=0; q<n_q_points; ++q)
+ for (unsigned int d=0; d<dim; ++d)
+ data.shape_gradients[k*dim+d][q] = shape_grads2[q][d];
+ };
+ }
+
+ if (flags & update_second_derivatives)
+ this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+}
+
+
+
+template <int dim>
+unsigned int
+FE_RaviartThomas<dim>::n_base_elements () const
+{
+ return 1;
+}
+
+
+
+template <int dim>
+const FiniteElement<dim> &
+FE_RaviartThomas<dim>::base_element (const unsigned int index) const
+{
+ Assert (index==0, ExcIndexRange(index, 0, 1));
+ return *this;
+}
+
+
+
+template <int dim>
+unsigned int
+FE_RaviartThomas<dim>::element_multiplicity (const unsigned int index) const
+{
+ Assert (index==0, ExcIndexRange(index, 0, 1));
+ return 1;
+}
+
+
+
+template <int dim>
+bool
+FE_RaviartThomas<dim>::has_support_on_face (const unsigned int shape_index,
+ const unsigned int face_index) const
+{
+ Assert (shape_index < this->dofs_per_cell,
+ ExcIndexRange (shape_index, 0, this->dofs_per_cell));
+ Assert (face_index < GeometryInfo<dim>::faces_per_cell,
+ ExcIndexRange (face_index, 0, GeometryInfo<dim>::faces_per_cell));
+
+ 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 <int dim>
+unsigned int
+FE_RaviartThomas<dim>::memory_consumption () const
+{
+ Assert (false, ExcNotImplemented ());
+ return 0;
+}
+
+
+
+template <int dim>
+unsigned int
+FE_RaviartThomas<dim>::get_degree () const
+{
+ return degree;
+}
+
+
+
+template class FE_RaviartThomas<deal_II_dimension>;
#include <fe/fe_system.h>
#include <fe/fe_values.h>
+#ifdef HAVE_STD_STRINGSTREAM
+# include <sstream>
+#else
+# include <strstream>
+#endif
/* ----------------------- FESystem::InternalData ------------------- */
+template <int dim>
+std::string
+FESystem<dim>::get_name () const
+{
+#ifdef HAVE_STD_STRINGSTREAM
+ std::ostringstream namebuf;
+#else
+ std::ostrstream namebuf;
+#endif
+
+ namebuf << "FESystem<" << dim << ">[";
+ for (unsigned int i=0; i<n_base_elements(); ++i)
+ {
+ namebuf << base_element(i).get_name();
+ if (element_multiplicity(i) != 1)
+ namebuf << '^' << element_multiplicity(i);
+ if (i != n_base_elements()-1)
+ namebuf << '-';
+ }
+ namebuf << ']';
+
+#ifndef HAVE_STD_STRINGSTREAM
+ namebuf << std::ends;
+#endif
+ return namebuf.str();
+}
+
+
+
template <int dim>
FiniteElement<dim>*
FESystem<dim>::clone() const
+template <int dim>
+void
+FESystem<dim>::
+get_interpolation_matrix (const FiniteElementBase<dim> &x_source_fe,
+ FullMatrix<double> &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<const FESystem<dim>*>(&x_source_fe) != 0),
+ typename FiniteElementBase<dim>::
+ ExcInterpolationNotImplemented());
+
+ // ok, source is a system element,
+ // so we may be able to do the work
+ const FESystem<dim> &source_fe
+ = dynamic_cast<const FESystem<dim>&>(x_source_fe);
+
+ // condition 2: same number of
+ // basis elements
+ AssertThrow (n_base_elements() == source_fe.n_base_elements(),
+ typename FiniteElementBase<dim>::
+ ExcInterpolationNotImplemented());
+
+ // condition 3: same number of
+ // basis elements
+ for (unsigned int i=0; i<n_base_elements(); ++i)
+ AssertThrow (element_multiplicity(i) ==
+ source_fe.element_multiplicity(i),
+ typename FiniteElementBase<dim>::
+ 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<FullMatrix<double> > base_matrices (n_base_elements());
+ for (unsigned int i=0; i<n_base_elements(); ++i)
+ {
+ base_matrices[i].reinit (base_element(i).dofs_per_cell,
+ source_fe.base_element(i).dofs_per_cell);
+ base_element(i).get_interpolation_matrix (source_fe.base_element(i),
+ base_matrices[i]);
+ }
+
+ // first clear big matrix, to make
+ // sure that entries that would
+ // couple different bases (or
+ // multiplicity indices) are really
+ // zero. then assign entries
+ interpolation_matrix.clear ();
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
+ if (this->system_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 <int dim>
UpdateFlags
FESystem<dim>::update_once (const UpdateFlags flags) const
#include <dofs/dof_accessor.h>
+
+namespace
+{
+ // forwarder function for
+ // FE::get_interpolation_matrix. we
+ // will want to call that function
+ // for arbitrary FullMatrix<T>
+ // 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 <int dim>
+ void gim_forwarder (const FiniteElement<dim> &fe1,
+ const FiniteElement<dim> &fe2,
+ FullMatrix<double> &interpolation_matrix)
+ {
+ fe2.get_interpolation_matrix (fe1, interpolation_matrix);
+ }
+
+
+ template <int dim, typename number>
+ void gim_forwarder (const FiniteElement<dim> &fe1,
+ const FiniteElement<dim> &fe2,
+ FullMatrix<number> &interpolation_matrix)
+ {
+ FullMatrix<double> tmp (interpolation_matrix.m(),
+ interpolation_matrix.n());
+ fe2.get_interpolation_matrix (fe1, tmp);
+ interpolation_matrix = tmp;
+ }
+}
+
+
template <int dim, typename number>
-void FETools::get_interpolation_matrix(const FiniteElement<dim> &fe1,
- const FiniteElement<dim> &fe2,
- FullMatrix<number> &interpolation_matrix)
+void FETools::get_interpolation_matrix (const FiniteElement<dim> &fe1,
+ const FiniteElement<dim> &fe2,
+ FullMatrix<number> &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(),
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<dim>::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.
{
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(),
{
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(),