#define __deal2__fe_dgp_h
#include <base/config.h>
-#include <base/polynomial.h>
#include <base/polynomial_space.h>
-#include <grid/geometry_info.h>
-#include <fe/fe.h>
+#include <fe/fe_poly.h>
-template <int dim> class PolynomialSpace;
template <int dim> class MappingQ;
-
/*!@addtogroup fe */
/*@{*/
* Discontinuous finite elements based on Legendre polynomials.
*
* This finite element implements complete polynomial spaces, that is,
- * $d$-dimensional polynomials of order $k$. The underlying
+ * dim-dimensional polynomials of degree p. The underlying
* polynomials form a Legendre basis on the unit square. Thus, the
* mass matrix is diagonal, if the grid cells are parallelograms.
*
- * @author Guido Kanschat, 2001, 2002
+ * @author Guido Kanschat, 2001, 2002, Ralf Hartmann 2004
*/
template <int dim>
-class FE_DGP : public FiniteElement<dim>
+class FE_DGP : public FE_Poly<PolynomialSpace<dim>,dim>
{
public:
/**
* Constructor for tensor product
- * polynomials of degree @p{k}.
+ * polynomials of degree @p p.
*/
- FE_DGP (const unsigned int k);
+ FE_DGP (const unsigned int p);
/**
* Return a string that uniquely
*/
virtual std::string get_name () const;
- /**
- * Return the value 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 (const unsigned int i,
- const Point<dim> &p) const;
-
- /**
- * Return the value of the
- * @p{component}th vector
- * component of the @p{i}th shape
- * function at the point
- * @p{p}. See the
- * @ref{FiniteElementBase} base
- * class for more information
- * about the semantics of this
- * function.
- *
- * Since this element is scalar,
- * the returned value is the same
- * as if the function without the
- * @p{_component} suffix were
- * called, provided that the
- * specified component is zero.
- */
- virtual double shape_value_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const;
-
- /**
- * Return the gradient of the
- * @p{i}th shape function at the
- * point @p{p}. See the
- * @ref{FiniteElementBase} base
- * class for more information
- * about the semantics of this
- * function.
- */
- virtual Tensor<1,dim> shape_grad (const unsigned int i,
- const Point<dim> &p) const;
-
- /**
- * Return the gradient of the
- * @p{component}th vector
- * component of the @p{i}th shape
- * function at the point
- * @p{p}. See the
- * @ref{FiniteElementBase} base
- * class for more information
- * about the semantics of this
- * function.
- *
- * Since this element is scalar,
- * the returned value is the same
- * as if the function without the
- * @p{_component} suffix were
- * called, provided that the
- * specified component is zero.
- */
- virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const;
-
- /**
- * Return the tensor of second
- * derivatives of the @p{i}th
- * shape function at point @p{p}
- * on the unit cell. See the
- * @ref{FiniteElementBase} base
- * class for more information
- * about the semantics of this
- * function.
- */
- virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
- const Point<dim> &p) const;
-
- /**
- * Return the second derivative
- * of the @p{component}th vector
- * component of the @p{i}th shape
- * function at the point
- * @p{p}. See the
- * @ref{FiniteElementBase} base
- * class for more information
- * about the semantics of this
- * function.
- *
- * Since this element is scalar,
- * the returned value is the same
- * as if the function without the
- * @p{_component} suffix were
- * called, provided that the
- * specified component is zero.
- */
- virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const;
-
/**
* Return the polynomial degree
* of this finite element,
* constructor.
*/
unsigned int get_degree () const;
-
- /**
- * Number of base elements in a
- * mixed discretization. Since
- * this is a scalar element,
- * return one.
- */
- virtual unsigned int n_base_elements () const;
-
- /**
- * Access to base element
- * objects. Since this element is
- * scalar, @p{base_element(0)} is
- * @p{this}, and all other
- * indices throw an error.
- */
- virtual const FiniteElement<dim> &
- base_element (const unsigned int index) const;
-
- /**
- * Multiplicity of base element
- * @p{index}. Since this is a
- * scalar element,
- * @p{element_multiplicity(0)}
- * returns one, and all other
- * indices will throw an error.
- */
- virtual unsigned int element_multiplicity (const unsigned int index) const;
/**
* Check for non-zero values on a face.
* 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:
* @p{FiniteElementData}.
*/
static std::vector<unsigned int> get_dpo_vector(unsigned int degree);
-
- /**
- * 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;
/**
* Degree of the polynomials.
*/
const unsigned int degree;
-
- /**
- * Pointer to an object
- * representing the polynomial
- * space used here.
- */
- const PolynomialSpace<dim> polynomial_space;
-
- /**
- * Fields of cell-independent data.
- *
- * For information about the
- * general purpose of this class,
- * see the documentation of the
- * base class.
- */
- class InternalData : public FiniteElementBase<dim>::InternalDataBase
- {
- public:
- /**
- * Array with shape function
- * values in quadrature
- * points. There is one
- * row for each shape
- * function, containing
- * values for each quadrature
- * point.
- *
- * In this array, we store
- * the values of the shape
- * function in the quadrature
- * points on the unit
- * cell. Since these values
- * do not change under
- * transformation to the real
- * cell, we only need to copy
- * them over when visiting a
- * concrete cell.
- */
- Table<2,double> shape_values;
-
- /**
- * Array with shape function
- * gradients in quadrature
- * points. There is one
- * row for each shape
- * function, containing
- * values for each quadrature
- * point.
- *
- * We store the gradients in
- * the quadrature points on
- * the unit cell. We then
- * only have to apply the
- * transformation (which is a
- * matrix-vector
- * multiplication) when
- * visiting an actual cell.
- */
- Table<2,Tensor<1,dim> > shape_gradients;
- };
/**
* Allow access from other dimensions.
#include <base/config.h>
#include <base/polynomials_p.h>
-#include <grid/geometry_info.h>
-#include <fe/fe.h>
+#include <fe/fe_poly.h>
template <int dim> class MappingQ;
/**
* Discontinuous finite elements based on monomials of degree up to
- * <tt>k</tt>
+ * <tt>p</tt>.
*
- * This finite element makes use of the @ref{PolynomialsP} class which
- * implements $d$-dimensional polynomials of degree $k$ based the
- * @ref{Polynomials::Monomial} and the @ref{PolynomialSpace} classes.
+ * This finite element makes use of the PolynomialsP class which
+ * implements <tt>dim</tt>-dimensional polynomials of degree
+ * <tt>p</tt> based the Polynomials::Polynomial and the
+ * PolynomialSpace classes.
*
* @author Ralf Hartmann, 2004
*/
template <int dim>
-class FE_DGPMonomial : public FiniteElement<dim>
+class FE_DGPMonomial : public FE_Poly<PolynomialsP<dim>,dim>
{
public:
/**
* Constructor for the polynomial
- * space of degree @p{k}.
+ * space of degree <tt>p</tt>.
*/
- FE_DGPMonomial (const unsigned int k);
+ FE_DGPMonomial (const unsigned int p);
/**
* Return a string that uniquely
* identifies a finite
* element. This class returns
- * @p{FE_DGPMonomial<dim>(degree)}, with
- * @p{dim} and @p{degree}
- * replaced by appropriate
- * values.
+ * <tt>FE_DGPMonomial<dim>(degree)</tt>,
+ * with <tt>dim</tt> and
+ * <tt>p</tt> replaced by
+ * appropriate values.
*/
virtual std::string get_name () const;
- /**
- * Return the value 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 (const unsigned int i,
- const Point<dim> &p) const;
-
- /**
- * Return the value of the
- * @p{component}th vector
- * component of the @p{i}th shape
- * function at the point
- * @p{p}. See the
- * @ref{FiniteElementBase} base
- * class for more information
- * about the semantics of this
- * function.
- *
- * Since this element is scalar,
- * the returned value is the same
- * as if the function without the
- * @p{_component} suffix were
- * called, provided that the
- * specified component is zero.
- */
- virtual double shape_value_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const;
-
- /**
- * Return the gradient of the
- * @p{i}th shape function at the
- * point @p{p}. See the
- * @ref{FiniteElementBase} base
- * class for more information
- * about the semantics of this
- * function.
- */
- virtual Tensor<1,dim> shape_grad (const unsigned int i,
- const Point<dim> &p) const;
-
- /**
- * Return the gradient of the
- * @p{component}th vector
- * component of the @p{i}th shape
- * function at the point
- * @p{p}. See the
- * @ref{FiniteElementBase} base
- * class for more information
- * about the semantics of this
- * function.
- *
- * Since this element is scalar,
- * the returned value is the same
- * as if the function without the
- * @p{_component} suffix were
- * called, provided that the
- * specified component is zero.
- */
- virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const;
-
- /**
- * Return the tensor of second
- * derivatives of the @p{i}th
- * shape function at point @p{p}
- * on the unit cell. See the
- * @ref{FiniteElementBase} base
- * class for more information
- * about the semantics of this
- * function.
- */
- virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
- const Point<dim> &p) const;
-
- /**
- * Return the second derivative
- * of the @p{component}th vector
- * component of the @p{i}th shape
- * function at the point
- * @p{p}. See the
- * @ref{FiniteElementBase} base
- * class for more information
- * about the semantics of this
- * function.
- *
- * Since this element is scalar,
- * the returned value is the same
- * as if the function without the
- * @p{_component} suffix were
- * called, provided that the
- * specified component is zero.
- */
- virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const;
-
/**
* Return the polynomial degree
* of this finite element,
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,
- * return one.
- */
- virtual unsigned int n_base_elements () const;
-
- /**
- * Access to base element
- * objects. Since this element is
- * scalar, @p{base_element(0)} is
- * @p{this}, and all other
- * indices throw an error.
- */
- virtual const FiniteElement<dim> &
- base_element (const unsigned int index) const;
-
- /**
- * Multiplicity of base element
- * @p{index}. Since this is a
- * scalar element,
- * @p{element_multiplicity(0)}
- * returns one, and all other
- * indices will throw an error.
- */
- virtual unsigned int element_multiplicity (const unsigned int index) const;
/**
* Check for non-zero values on a face.
* 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:
* constructor.
*/
void initialize_restriction ();
-
- /**
- * 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;
-
- /**
- * Pointer to an object
- * representing the polynomial
- * space <tt>P_k</tt>.
- */
- const PolynomialsP<dim> polynomial_space;
-
- /**
- * Fields of cell-independent data.
- *
- * For information about the
- * general purpose of this class,
- * see the documentation of the
- * base class.
- */
- class InternalData : public FiniteElementBase<dim>::InternalDataBase
- {
- public:
- /**
- * Array with shape function
- * values in quadrature
- * points. There is one
- * row for each shape
- * function, containing
- * values for each quadrature
- * point.
- *
- * In this array, we store
- * the values of the shape
- * function in the quadrature
- * points on the unit
- * cell. Since these values
- * do not change under
- * transformation to the real
- * cell, we only need to copy
- * them over when visiting a
- * concrete cell.
- */
- Table<2,double> shape_values;
-
- /**
- * Array with shape function
- * gradients in quadrature
- * points. There is one
- * row for each shape
- * function, containing
- * values for each quadrature
- * point.
- *
- * We store the gradients in
- * the quadrature points on
- * the unit cell. We then
- * only have to apply the
- * transformation (which is a
- * matrix-vector
- * multiplication) when
- * visiting an actual cell.
- */
- Table<2,Tensor<1,dim> > shape_gradients;
- };
/**
* Allows @p{MappingQ} class
inline unsigned int
FE_DGPMonomial<dim>::get_degree () const
{
- return polynomial_space.degree();
+ return this->poly_space.degree();
}
// $Id$
// Version: $Name$
//
-// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 by the deal.II authors
+// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
template <int dim>
FE_DGP<dim>::FE_DGP (const unsigned int degree)
:
- FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree), 1, degree),
- std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(degree), 1, degree).dofs_per_cell,true),
- std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree), 1, degree).dofs_per_cell,
- std::vector<bool>(1,true))),
- degree(degree),
- polynomial_space (Polynomials::Legendre::generate_complete_basis(degree))
+ FE_Poly<PolynomialSpace<dim>, dim> (
+ PolynomialSpace<dim>(Polynomials::Legendre::generate_complete_basis(degree)),
+ FiniteElementData<dim>(get_dpo_vector(degree), 1, degree),
+ std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(degree), 1, degree).dofs_per_cell,true),
+ std::vector<std::vector<bool> >(FiniteElementData<dim>(
+ get_dpo_vector(degree), 1, degree).dofs_per_cell, std::vector<bool>(1,true))),
+ degree(degree)
{
// if defined, copy over matrices
// from precomputed arrays
-template <int dim>
-double
-FE_DGP<dim>::shape_value (const unsigned int i,
- const Point<dim> &p) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
- return polynomial_space.compute_value(i, p);
-}
-
-
-
-template <int dim>
-double
-FE_DGP<dim>::shape_value_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
- Assert (component == 0, ExcIndexRange (component, 0, 1));
- return polynomial_space.compute_value(i, p);
-}
-
-
-
-template <int dim>
-Tensor<1,dim>
-FE_DGP<dim>::shape_grad (const unsigned int i,
- const Point<dim> &p) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
- return polynomial_space.compute_grad(i, p);
-}
-
-
-template <int dim>
-Tensor<1,dim>
-FE_DGP<dim>::shape_grad_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
- Assert (component == 0, ExcIndexRange (component, 0, 1));
- return polynomial_space.compute_grad(i, p);
-}
-
-
-
-template <int dim>
-Tensor<2,dim>
-FE_DGP<dim>::shape_grad_grad (const unsigned int i,
- const Point<dim> &p) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
- return polynomial_space.compute_grad_grad(i, p);
-}
-
-
-
-template <int dim>
-Tensor<2,dim>
-FE_DGP<dim>::shape_grad_grad_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
- Assert (component == 0, ExcIndexRange (component, 0, 1));
- return polynomial_space.compute_grad_grad(i, p);
-}
-
-
//----------------------------------------------------------------------
// Auxiliary functions
//----------------------------------------------------------------------
}
-template <int dim>
-UpdateFlags
-FE_DGP<dim>::update_once (const UpdateFlags flags) const
-{
- // for this kind of elements, only
- // the values can be precomputed
- // once and for all. set this flag
- // if the values are requested at
- // all
- return (update_default | (flags & update_values));
-}
-
-
-template <int dim>
-UpdateFlags
-FE_DGP<dim>::update_each (const UpdateFlags flags) const
-{
- UpdateFlags out = update_default;
-
- if (flags & update_gradients)
- out |= update_gradients | update_covariant_transformation;
-
- if (flags & update_second_derivatives)
- out |= update_second_derivatives | update_covariant_transformation;
-
- return out;
-}
-
-
-//----------------------------------------------------------------------
-// Data field initialization
-//----------------------------------------------------------------------
-
-template <int dim>
-typename Mapping<dim>::InternalDataBase *
-FE_DGP<dim>::get_data (const UpdateFlags update_flags,
- const Mapping<dim> &mapping,
- const Quadrature<dim> &quadrature) const
-{
- // generate a new data object
- 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;
-
- // have some scratch arrays
- std::vector<double> values(0);
- std::vector<Tensor<1,dim> > grads(0);
- std::vector<Tensor<2,dim> > grad_grads(0);
-
- // initialize fields only if really
- // necessary. otherwise, don't
- // allocate memory
- if (flags & update_values)
- {
- values.resize (this->dofs_per_cell);
- data->shape_values.reinit (this->dofs_per_cell,
- n_q_points);
- }
-
- if (flags & update_gradients)
- {
- grads.resize (this->dofs_per_cell);
- data->shape_gradients.reinit (this->dofs_per_cell,
- n_q_points);
- }
-
- // if second derivatives through
- // finite differencing is required,
- // then initialize some objects for
- // that
- if (flags & update_second_derivatives)
- data->initialize_2nd (this, mapping, quadrature);
-
- // next already fill those fields
- // of which we have information by
- // now. note that the shape
- // gradients are only those on the
- // unit cell, and need to be
- // transformed when visiting an
- // actual cell
- if (flags & (update_values | update_gradients))
- for (unsigned int i=0; i<n_q_points; ++i)
- {
- polynomial_space.compute(quadrature.point(i),
- values, grads, grad_grads);
- for (unsigned int k=0; k<this->dofs_per_cell; ++k)
- {
- if (flags & update_values)
- data->shape_values[k][i] = values[k];
- if (flags & update_gradients)
- data->shape_gradients[k][i] = grads[k];
- }
- }
- return data;
-}
-
-
-
-//----------------------------------------------------------------------
-// Fill data of FEValues
-//----------------------------------------------------------------------
-
-template <int dim>
-void
-FE_DGP<dim>::fill_fe_values (const Mapping<dim> &mapping,
- const typename DoFHandler<dim>::cell_iterator &cell,
- const Quadrature<dim> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename Mapping<dim>::InternalDataBase &fedata,
- FEValuesData<dim> &data) const
-{
- // convert data object to internal
- // data for this class. fails with
- // an exception if that is not
- // possible
- InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-
- const UpdateFlags flags(fe_data.current_update_flags());
-
- for (unsigned int k=0; k<this->dofs_per_cell; ++k)
- {
- if (flags & update_values)
- for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
- data.shape_values(k,i) = fe_data.shape_values[k][i];
-
- if (flags & update_gradients)
- {
- Assert (data.shape_gradients[k].size() <=
- fe_data.shape_gradients[k].size(),
- ExcInternalError());
- mapping.transform_covariant(data.shape_gradients[k].begin(),
- data.shape_gradients[k].end(),
- fe_data.shape_gradients[k].begin(),
- mapping_data);
- }
- }
-
- if (flags & update_second_derivatives)
- this->compute_2nd (mapping, cell,
- QProjector<dim>::DataSetDescriptor::cell(),
- mapping_data, fe_data, data);
-}
-
-
-
-template <int dim>
-void
-FE_DGP<dim>::fill_fe_face_values (const Mapping<dim> &mapping,
- const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face,
- const Quadrature<dim-1> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename Mapping<dim>::InternalDataBase &fedata,
- FEValuesData<dim> &data) const
-{
- // convert data object to internal
- // data for this class. fails with
- // an exception if that is not
- // possible
- InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-
- // offset determines which data set
- // to take (all data sets for all
- // faces are stored contiguously)
- const typename QProjector<dim>::DataSetDescriptor offset
- = (QProjector<dim>::DataSetDescriptor::
- face (face, cell->face_orientation(face),
- quadrature.n_quadrature_points));
-
- const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
-
- for (unsigned int k=0; k<this->dofs_per_cell; ++k)
- {
- for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
- if (flags & update_values)
- data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
-
- if (flags & update_gradients)
- {
- Assert (data.shape_gradients[k].size() + offset <=
- fe_data.shape_gradients[k].size(),
- ExcInternalError());
- mapping.transform_covariant(data.shape_gradients[k].begin(),
- data.shape_gradients[k].end(),
- fe_data.shape_gradients[k].begin()+offset,
- mapping_data);
- };
- }
-
- if (flags & update_second_derivatives)
- this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
-}
-
-
-
-template <int dim>
-void
-FE_DGP<dim>::fill_fe_subface_values (const Mapping<dim> &mapping,
- const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face,
- const unsigned int subface,
- const Quadrature<dim-1> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename Mapping<dim>::InternalDataBase &fedata,
- FEValuesData<dim> &data) const
-{
- // convert data object to internal
- // data for this class. fails with
- // an exception if that is not
- // possible
- InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-
- // offset determines which data set
- // to take (all data sets for all
- // sub-faces are stored contiguously)
- const typename QProjector<dim>::DataSetDescriptor offset
- = (QProjector<dim>::DataSetDescriptor::
- sub_face (face, subface, cell->face_orientation(face),
- quadrature.n_quadrature_points));
-
- const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
-
- for (unsigned int k=0; k<this->dofs_per_cell; ++k)
- {
- for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
- if (flags & update_values)
- data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
-
- if (flags & update_gradients)
- {
- Assert (data.shape_gradients[k].size() + offset <=
- fe_data.shape_gradients[k].size(),
- ExcInternalError());
- mapping.transform_covariant(data.shape_gradients[k].begin(),
- data.shape_gradients[k].end(),
- fe_data.shape_gradients[k].begin()+offset,
- mapping_data);
- };
- }
-
- if (flags & update_second_derivatives)
- this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
-}
-
-
-
-template <int dim>
-unsigned int
-FE_DGP<dim>::n_base_elements () const
-{
- return 1;
-}
-
-
-
-template <int dim>
-const FiniteElement<dim> &
-FE_DGP<dim>::base_element (const unsigned int index) const
-{
- Assert (index==0, ExcIndexRange(index, 0, 1));
- return *this;
-}
-
-
-
-template <int dim>
-unsigned int
-FE_DGP<dim>::element_multiplicity (const unsigned int index) const
-{
- Assert (index==0, ExcIndexRange(index, 0, 1));
- return 1;
-}
-
-
-
template <int dim>
bool
FE_DGP<dim>::has_support_on_face (const unsigned int,
template <int dim>
FE_DGPMonomial<dim>::FE_DGPMonomial (const unsigned int degree)
:
- FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree), 1, degree),
- std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(degree), 1, degree).dofs_per_cell,true),
- std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree), 1, degree).dofs_per_cell,
- std::vector<bool>(1,true))),
- polynomial_space(degree)
+ FE_Poly<PolynomialsP<dim>, dim> (
+ PolynomialsP<dim>(degree),
+ FiniteElementData<dim>(get_dpo_vector(degree), 1, degree),
+ std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(degree), 1, degree).dofs_per_cell,true),
+ std::vector<std::vector<bool> >(FiniteElementData<dim>(
+ get_dpo_vector(degree), 1, degree).dofs_per_cell, std::vector<bool>(1,true)))
{
- Assert(polynomial_space.n()==this->dofs_per_cell, ExcInternalError());
+ Assert(this->poly_space.n()==this->dofs_per_cell, ExcInternalError());
// DG doesn't have constraints, so
// leave them empty
-template <int dim>
-double
-FE_DGPMonomial<dim>::shape_value (const unsigned int i,
- const Point<dim> &p) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
- return polynomial_space.compute_value(i, p);
-}
-
-
-
-template <int dim>
-double
-FE_DGPMonomial<dim>::shape_value_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
- Assert (component == 0, ExcIndexRange (component, 0, 1));
- return polynomial_space.compute_value(i, p);
-}
-
-
-
-template <int dim>
-Tensor<1,dim>
-FE_DGPMonomial<dim>::shape_grad (const unsigned int i,
- const Point<dim> &p) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
- return polynomial_space.compute_grad(i, p);
-}
-
-
-template <int dim>
-Tensor<1,dim>
-FE_DGPMonomial<dim>::shape_grad_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
- Assert (component == 0, ExcIndexRange (component, 0, 1));
- return polynomial_space.compute_grad(i, p);
-}
-
-
-
-template <int dim>
-Tensor<2,dim>
-FE_DGPMonomial<dim>::shape_grad_grad (const unsigned int i,
- const Point<dim> &p) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
- return polynomial_space.compute_grad_grad(i, p);
-}
-
-
-
-template <int dim>
-Tensor<2,dim>
-FE_DGPMonomial<dim>::shape_grad_grad_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
- Assert (component == 0, ExcIndexRange (component, 0, 1));
- return polynomial_space.compute_grad_grad(i, p);
-}
-
-
-
template <int dim>
void
FE_DGPMonomial<dim>::
for (unsigned int i=0; i<this->dofs_per_cell; ++i)
{
- cell_interpolation(j,i) = polynomial_space.compute_value (i, p_cell);
- subcell_interpolation(j,i) = polynomial_space.compute_value (i, p_subcell);
+ cell_interpolation(j,i) = this->poly_space.compute_value (i, p_cell);
+ subcell_interpolation(j,i) = this->poly_space.compute_value (i, p_subcell);
}
}
}
-template <int dim>
-UpdateFlags
-FE_DGPMonomial<dim>::update_once (const UpdateFlags flags) const
-{
- // for this kind of elements, only
- // the values can be precomputed
- // once and for all. set this flag
- // if the values are requested at
- // all
- return (update_default | (flags & update_values));
-}
-
-
-template <int dim>
-UpdateFlags
-FE_DGPMonomial<dim>::update_each (const UpdateFlags flags) const
-{
- UpdateFlags out = update_default;
-
- if (flags & update_gradients)
- out |= update_gradients | update_covariant_transformation;
-
- if (flags & update_second_derivatives)
- out |= update_second_derivatives | update_covariant_transformation;
-
- return out;
-}
-
-
-//----------------------------------------------------------------------
-// Data field initialization
-//----------------------------------------------------------------------
-
-template <int dim>
-typename Mapping<dim>::InternalDataBase *
-FE_DGPMonomial<dim>::get_data (const UpdateFlags update_flags,
- const Mapping<dim> &mapping,
- const Quadrature<dim> &quadrature) const
-{
- // generate a new data object
- 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;
-
- // have some scratch arrays
- std::vector<double> values(0);
- std::vector<Tensor<1,dim> > grads(0);
- std::vector<Tensor<2,dim> > grad_grads(0);
-
- // initialize fields only if really
- // necessary. otherwise, don't
- // allocate memory
- if (flags & update_values)
- {
- values.resize (this->dofs_per_cell);
- data->shape_values.reinit (this->dofs_per_cell,
- n_q_points);
- }
-
- if (flags & update_gradients)
- {
- grads.resize (this->dofs_per_cell);
- data->shape_gradients.reinit (this->dofs_per_cell,
- n_q_points);
- }
-
- // if second derivatives through
- // finite differencing is required,
- // then initialize some objects for
- // that
- if (flags & update_second_derivatives)
- data->initialize_2nd (this, mapping, quadrature);
-
- // next already fill those fields
- // of which we have information by
- // now. note that the shape
- // gradients are only those on the
- // unit cell, and need to be
- // transformed when visiting an
- // actual cell
- if (flags & (update_values | update_gradients))
- for (unsigned int i=0; i<n_q_points; ++i)
- {
- polynomial_space.compute(quadrature.point(i),
- values, grads, grad_grads);
- for (unsigned int k=0; k<this->dofs_per_cell; ++k)
- {
- if (flags & update_values)
- data->shape_values[k][i] = values[k];
- if (flags & update_gradients)
- data->shape_gradients[k][i] = grads[k];
- }
- }
- return data;
-}
-
-
-
-//----------------------------------------------------------------------
-// Fill data of FEValues
-//----------------------------------------------------------------------
-
-template <int dim>
-void
-FE_DGPMonomial<dim>::fill_fe_values (const Mapping<dim> &mapping,
- const typename DoFHandler<dim>::cell_iterator &cell,
- const Quadrature<dim> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename Mapping<dim>::InternalDataBase &fedata,
- FEValuesData<dim> &data) const
-{
- // convert data object to internal
- // data for this class. fails with
- // an exception if that is not
- // possible
- InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-
- const UpdateFlags flags(fe_data.current_update_flags());
-
- for (unsigned int k=0; k<this->dofs_per_cell; ++k)
- {
- if (flags & update_values)
- for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
- data.shape_values(k,i) = fe_data.shape_values[k][i];
-
- if (flags & update_gradients)
- {
- Assert (data.shape_gradients[k].size() <=
- fe_data.shape_gradients[k].size(),
- ExcInternalError());
- mapping.transform_covariant(data.shape_gradients[k].begin(),
- data.shape_gradients[k].end(),
- fe_data.shape_gradients[k].begin(),
- mapping_data);
- }
- }
-
- if (flags & update_second_derivatives)
- this->compute_2nd (mapping, cell,
- QProjector<dim>::DataSetDescriptor::cell(),
- mapping_data, fe_data, data);
-}
-
-
-
-template <int dim>
-void
-FE_DGPMonomial<dim>::fill_fe_face_values (const Mapping<dim> &mapping,
- const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face,
- const Quadrature<dim-1> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename Mapping<dim>::InternalDataBase &fedata,
- FEValuesData<dim> &data) const
-{
- // convert data object to internal
- // data for this class. fails with
- // an exception if that is not
- // possible
- InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-
- // offset determines which data set
- // to take (all data sets for all
- // faces are stored contiguously)
- const typename QProjector<dim>::DataSetDescriptor offset
- = (QProjector<dim>::DataSetDescriptor::
- face (face, cell->face_orientation(face),
- quadrature.n_quadrature_points));
-
- const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
-
- for (unsigned int k=0; k<this->dofs_per_cell; ++k)
- {
- for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
- if (flags & update_values)
- data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
-
- if (flags & update_gradients)
- {
- Assert (data.shape_gradients[k].size() + offset <=
- fe_data.shape_gradients[k].size(),
- ExcInternalError());
- mapping.transform_covariant(data.shape_gradients[k].begin(),
- data.shape_gradients[k].end(),
- fe_data.shape_gradients[k].begin()+offset,
- mapping_data);
- };
- }
-
- if (flags & update_second_derivatives)
- this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
-}
-
-
-
-template <int dim>
-void
-FE_DGPMonomial<dim>::fill_fe_subface_values (const Mapping<dim> &mapping,
- const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face,
- const unsigned int subface,
- const Quadrature<dim-1> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename Mapping<dim>::InternalDataBase &fedata,
- FEValuesData<dim> &data) const
-{
- // convert data object to internal
- // data for this class. fails with
- // an exception if that is not
- // possible
- InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-
- // offset determines which data set
- // to take (all data sets for all
- // sub-faces are stored contiguously)
- const typename QProjector<dim>::DataSetDescriptor offset
- = (QProjector<dim>::DataSetDescriptor::
- sub_face (face, subface, cell->face_orientation(face),
- quadrature.n_quadrature_points));
-
- const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
-
- for (unsigned int k=0; k<this->dofs_per_cell; ++k)
- {
- for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
- if (flags & update_values)
- data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
-
- if (flags & update_gradients)
- {
- Assert (data.shape_gradients[k].size() + offset <=
- fe_data.shape_gradients[k].size(),
- ExcInternalError());
- mapping.transform_covariant(data.shape_gradients[k].begin(),
- data.shape_gradients[k].end(),
- fe_data.shape_gradients[k].begin()+offset,
- mapping_data);
- };
- }
-
- if (flags & update_second_derivatives)
- this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
-}
-
-
-
-template <int dim>
-unsigned int
-FE_DGPMonomial<dim>::n_base_elements () const
-{
- return 1;
-}
-
-
-
-template <int dim>
-const FiniteElement<dim> &
-FE_DGPMonomial<dim>::base_element (const unsigned int index) const
-{
- Assert (index==0, ExcIndexRange(index, 0, 1));
- return *this;
-}
-
-
-
-template <int dim>
-unsigned int
-FE_DGPMonomial<dim>::element_multiplicity (const unsigned int index) const
-{
- Assert (index==0, ExcIndexRange(index, 0, 1));
- return 1;
-}
-
-
#if deal_II_dimension == 1
template <>
else
{
unsigned int degrees[2];
- polynomial_space.directional_degrees(shape_index, degrees);
+ this->poly_space.directional_degrees(shape_index, degrees);
if ((face_index==0 && degrees[1]==0) ||
(face_index==3 && degrees[0]==0))
support_on_face=true;
else
{
unsigned int degrees[3];
- polynomial_space.directional_degrees(shape_index, degrees);
+ this->poly_space.directional_degrees(shape_index, degrees);
if ((face_index==0 && degrees[1]==0) ||
(face_index==2 && degrees[2]==0) ||
(face_index==5 && degrees[0]==0))