// $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
#define __deal2__fe_dgq_h
#include <base/config.h>
-#include <base/polynomial.h>
#include <base/tensor_product_polynomials.h>
-#include <fe/fe.h>
+#include <fe/fe_poly.h>
-template <int dim> class TensorProductPolynomials;
template <int dim> class MappingQ;
-
/*!@addtogroup fe */
/*@{*/
/**
- * Discontinuous tensor product elements based on equidistant support points.
+ * Discontinuous tensor product elements based on equidistant support
+ * points.
*
- * This is a discontinuous finite element using interpolating tensor
- * product polynomials. The shape functions are Lagrangian
+ * This is a discontinuous finite element based on tensor products of
+ * Lagrangian polynomials. The shape functions are Lagrangian
* interpolants of an equidistant grid of points on the unit cell. The
- * points are numbered in lexicographical order, @p{x} running fastest.
+ * points are numbered in lexicographical order, with @p x running
+ * fastest.
*
- * @author Guido Kanschat, Ralf Hartmann, 2001
+ * @author Ralf Hartmann, 2001, 2004, Guido Kanschat 2001
*/
template <int dim>
-class FE_DGQ : public FiniteElement<dim>
+class FE_DGQ : public FE_Poly<TensorProductPolynomials<dim>,dim>
{
public:
/**
* Constructor for tensor product
- * polynomials of degree @p{k}.
+ * polynomials of degree @p p.
*/
- FE_DGQ (const unsigned int k);
+ FE_DGQ (const unsigned int p);
/**
* Return a string that uniquely
* identifies a finite
* element. This class returns
- * @p{FE_DGQ<dim>(degree)}, with
- * @p{dim} and @p{degree}
+ * @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
- * 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:
/**
* @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;
/**
* Compute renumbering for rotation
* Degree of the polynomials.
*/
const unsigned int degree;
-
- /**
- * Pointer to the tensor
- * product polynomials.
- */
- const TensorProductPolynomials<dim> polynomial_space;
-
- /**
- * Fields of cell-independent data.
- *
- * For information about the
- * general purpose of this class,
- * see the documentation of the
- * base class.
- */
- class InternalData : public FiniteElementBase<dim>::InternalDataBase
- {
- public:
- /**
- * Array with shape function
- * values in quadrature
- * points. There is one
- * row for each shape
- * function, containing
- * values for each quadrature
- * point.
- *
- * In this array, we store
- * the values of the shape
- * function in the quadrature
- * points on the unit
- * cell. Since these values
- * do not change under
- * transformation to the real
- * cell, we only need to copy
- * them over when visiting a
- * concrete cell.
- */
- Table<2,double> shape_values;
-
- /**
- * Array with shape function
- * gradients in quadrature
- * points. There is one
- * row for each shape
- * function, containing
- * values for each quadrature
- * point.
- *
- * We store the gradients in
- * the quadrature points on
- * the unit cell. We then
- * only have to apply the
- * transformation (which is a
- * matrix-vector
- * multiplication) when
- * visiting an actual cell.
- */
- Table<2,Tensor<1,dim> > shape_gradients;
- };
/**
* Allow access from other dimensions.
// $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_DGQ<dim>::FE_DGQ (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::LagrangeEquidistant::generate_complete_basis(degree))
+ FE_Poly<TensorProductPolynomials<dim>, dim> (
+ TensorProductPolynomials<dim>(Polynomials::LagrangeEquidistant::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)
{
// generate permutation/rotation
// index sets to generate some
}
-
-template <int dim>
-double
-FE_DGQ<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_DGQ<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_DGQ<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_DGQ<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_DGQ<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_DGQ<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_DGQ<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_DGQ<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;
-}
-
-
template <int dim>
void
int2type<dim>());
for (unsigned int i=0; i<this->dofs_per_cell; ++i)
cell_interpolation(j,i)
- = polynomial_space.compute_value (i, p);
+ = this->poly_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);
+ = source_fe.poly_space.compute_value (i, p);
}
// then compute the
-
-
-//----------------------------------------------------------------------
-// Data field initialization
-//----------------------------------------------------------------------
-
-template <int dim>
-typename Mapping<dim>::InternalDataBase *
-FE_DGQ<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_DGQ<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_DGQ<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_DGQ<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_DGQ<dim>::n_base_elements () const
-{
- return 1;
-}
-
-
-
-template <int dim>
-const FiniteElement<dim> &
-FE_DGQ<dim>::base_element (const unsigned int index) const
-{
- Assert (index==0, ExcIndexRange(index, 0, 1));
- return *this;
-}
-
-
-
-template <int dim>
-unsigned int
-FE_DGQ<dim>::element_multiplicity (const unsigned int index) const
-{
- Assert (index==0, ExcIndexRange(index, 0, 1));
- return 1;
-}
-
-
template <int dim>
bool